EnKF/EnKF-MPI-TOPAZ/Tools/p_fixhycom_eco.F90

! File:          p_fixhycom.F90
!
! Created:       ???
!
! Last modified: 29/06/2010
!
! Purpose:       Fixes EnKF output.
!
! Description:   
!
! Modifications:
!                1/03/2011 Ehouarn:
!                  -modification of the fixhycom subroutine: interpolation of 
!                    biogeochemical tracers on the analysis grid according to
!                    hycom remapping in order to insure conservation.  
!                29/06/2010 PS:
!                  - set the maximum ICEC to 0.995 to match the model
!                ?/?/? KAL:
!                  - Modification of the "fixhycom" subroutine, into separate
!                    program, working on a file-by-file basis
!                Prior history:
!                  Not documented.

! KAL -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
! KAL -- Input arguments:
! KAL --     template restart file
! KAL --     ensemble member
! KAL -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
! KAL -- TODO: Fix for ice fields

program fixhycom_eco




   use mod_raw_io
   use m_parse_blkdat
   use m_put_mod_fld
   use m_get_mod_fld
   use m_get_mod_grid
#if defined (ECO)
   use m_fixhycom_eco_metno
#endif   
   implicit none

   integer*4, external :: iargc
   real, parameter :: onem=9806.
   real, parameter :: PSTARB0 = 1000;

   integer imem                  ! ensemble member
   character(len=80) :: restart,icerestart ! restart template

   character(len=80) :: afile,newfile, char80
   integer          :: fnd, rstind, tmpindx, iafile
   logical          :: ex, allok, nomatch
   character(len=8) :: cfld, ctmp
   character(len=3) :: cproc,cmem
   integer          :: tlevel, vlevel, nproc
   real             :: bmin, bmax, rdummy
   integer          :: idm,jdm,kdm
   real, allocatable:: fld(:,:)
   real*8, allocatable, dimension(:,:) :: &
      ficem,hicem,hsnwm,ticem,tsrfm
   real, allocatable, dimension(:,:) :: depths,modlon,modlat,saln
   real*4, allocatable:: fldr4(:,:)
   real*4 :: spval,amin,amax
   real, allocatable :: press(:)

   real, allocatable, dimension(:,:)   :: dpsum
   real, allocatable, dimension(:,:,:)   :: dp, dpold

   integer,parameter :: numfields=2
   integer :: ios,ios2, reclICE,ifld
   integer :: i,j,k,ktr

   real :: mindx,meandx
   
#if defined (ECO)   
    real,dimension(:,:,:), allocatable::dpfor,cfi 
    character(len=80) :: restfor
    real,dimension(:,:,:,:), allocatable::tracerf
    real, dimension(:,:), allocatable::trcraij
    real, dimension(:), allocatable::prsf,dpf
    integer::ntracr,ktrcr
    real::dpthin
    character(2)::ctrcr
    logical, dimension(:), allocatable::lcm
    
    real, dimension(:,:,:), allocatable::temp,sal
    integer::kisop
#endif

   icerestart=''
   if (iargc()==2 .or. iargc()==3) then
      call getarg(1,restart)
      call getarg(2,ctmp)
      read(ctmp,*) imem
      write(cmem,'(i3.3)') imem
      if (iargc()==3) call getarg(3,icerestart)
   else
      print *,'"fixhycom" -- A crude routine to correct restart files for obvious errors'
      print *
      print *,'usage: '
      print *,'   fixhycom restart_file ensemble_member <ice_file>'
      print *,'   "restart_file"    the restart file you want to fix (.a-file)"'
      print *,'   "ensemble_member" is the ensemble member - should corr. to that of restart file'
      print *,'   "ice_file"        is optional - it is the restart file for ice fields'
      call exit(1)
   endif

   ! Get dimensions from blkdat
   call parse_blkdat('idm   ','integer',rdummy,idm)
   call parse_blkdat('jdm   ','integer',rdummy,jdm)
   call parse_blkdat('kdm   ','integer',rdummy,kdm)
#if defined (ECO)
   call parse_blkdat('ntracr','integer',rdummy,ntracr) 	  
#endif
      
   if (idm>0 .and. idm < 1e4 .and. jdm>0 .and. jdm<1e4) then
      allocate(fld  (idm,jdm))
      allocate(fldr4(idm,jdm))
      allocate(saln (idm,jdm))
      allocate(ficem(idm,jdm))
      allocate(hicem(idm,jdm))
      allocate(hsnwm(idm,jdm))
      allocate(ticem(idm,jdm))
      allocate(tsrfm(idm,jdm))
      allocate(dpsum(idm,jdm))
      allocate(depths(idm,jdm))
      allocate(modlon(idm,jdm))
      allocate(modlat(idm,jdm))
      allocate(dpold(idm,jdm,kdm))
      allocate(dp   (idm,jdm,kdm))
      allocate(press(kdm+1))
#if defined (ECO)
      allocate(dpfor(idm,jdm,kdm))
      allocate(cfi(kdm,ntracr,2))
      allocate(prsf(kdm+1))
      allocate(tracerf(idm,jdm,kdm,ntracr))
      allocate(trcraij(kdm,ntracr))
      allocate(lcm(kdm))
      allocate(dpf(kdm))
      allocate(sal(idm,jdm,kdm))
      allocate(temp(idm,jdm,kdm))
#endif
   else
      print *,'fld allocate error'
      stop '(EnKF_postprocess)'
   end if


   ! Remove postfix of restart file
   fnd=max(index(restart,'.a'),index(restart,'.b'))


   ! Inquire for existence
   inquire(exist=ex,file=restart(1:fnd-1)//'.b')
   if (.not.ex) then
      write(*,*) 'Can not find '//restart(1:fnd-1)//'.b'
      stop '(EnKF_postprocess)'
   end if

   print *,restart(1:fnd-1)//'.b'
   newfile='fix'//restart(1:fnd-1)


   ! Get model grid
   call get_mod_grid(modlon,modlat,depths,mindx,meandx,idm,jdm)

#if defined (ECO)
   !files where are stored the forecast fields!
   restfor='forecast'//cmem
   dpthin = onem*0.001
#endif

   ! Get layer thickness
   dpsum=0.
   do k=1,kdm
      !call get_mod_fld(dp(:,:,k),1,'dp      ',k,1)
      call get_mod_fld_new(restart(1:fnd-1),dp(:,:,k),imem,'dp      ',k,1,idm,jdm)
      dpsum=dpsum+dp(:,:,k)
#if defined (ECO)
      !reading of forecast fields: tracers, T and S
      call get_mod_fld_new(trim(restfor),dpfor(:,:,k),imem,'dp    ',k,1,idm,jdm)
      do ktrcr=1,ntracr
         write(ctrcr,'(i2.2)') ktrcr
         cfld='tracer'//ctrcr
         call get_mod_fld_new(trim(restfor),tracerf(:,:,k,ktrcr),imem,cfld,k,1,idm,jdm)
      enddo
      call get_mod_fld_new(trim(restfor),temp(:,:,k),imem,'temp    ',k,1,idm,jdm)
      call get_mod_fld_new(trim(restfor),sal(:,:,k),imem,'saln    ',k,1,idm,jdm)
#endif      
   end do
   print *,maxval(dpsum-depths*onem)
   dpold=dp



   ! DP correction
   do j=1,jdm
   do i=1,idm

      !!! Move negative layers to neighbouring layers.
      do k = 1, kdm-1
         dp(i,j,k+1) = dp(i,j,k+1) + min(0.0,dp(i,j,k))
         dp(i,j,k  ) = max(dp(i,j,k),0.0)
         !dp(i,j,k  ) = max(dp(i,j,k),1.e-3*onem)
      end do

      !!! Go backwards to fix lowermost layer.
      do k = kdm, 3, -1
         dp(i,j,k-1) = dp(i,j,k-1) + min(0.0,dp(i,j,k))
         dp(i,j,k)   =   max(dp(i,j,k),0.0)
	 !dp(i,j,k)   =   max(dp(i,j,k),1.e-3*onem)
      end do

      !!! No layers below the sea bed.
      press(  1) = 0.0
#if defined (ECO)
      !computation of the forecast layer interfaces (prsf)!
      prsf(1)=0.
#endif               
      do k = 1, kdm-1
         press(k+1) = press(k) + dp(i,j,k)
         press(k+1) = min(depths(i,j)*onem,press(k+1))
#if defined (ECO)
         prsf(k+1) = prsf(k) + dpfor(i,j,k)
#endif 	 
      end do
      press(kdm+1) = depths(i,j)*onem
#if defined (ECO)
      prsf(kdm+1)=depths(i,j)*onem
#endif     

      do k = 1, kdm
         dp(i,j,k) = press(k+1) - press(k)
#if defined (ECO)
         !definition of the isopycnal layers!
         dpf(k)=max(dpfor(i,j,k),dpthin) 
         kisop=compute_kisop(temp(i,j,:),sal(i,j,:),kdm)
         if     (k.le.max(2,kisop)) then
          lcm(k) = .false.  !fixed layers are never PCM
         else
! ---       thin and isopycnal layers remapped with PCM.
            lcm(k) = dpfor(i,j,k).le.dpthin
         endif 
#endif  	       
      end do

!eho 2/11/11      
      if(depths(i,j)>10000. .or. depths(i,j)<1.)then
        dp(i,j,:)=dpold(i,j,:)
      endif
      
#if defined (ECO)     
!      if(depths(i,j)==0.)then
      if(depths(i,j)>10000. .or. depths(i,j)<1.)then
        cycle
      else	
	call hybgen_weno_coefs(tracerf(i,j,:,:),dpf,lcm,cfi,kdm,ntracr,dpthin)
        call hybgen_weno_remap(tracerf(i,j,:,:),prsf,dpfor(i,j,:),cfi,trcraij,&
                        press,kdm,kdm,ntracr,dpthin)
        tracerf(i,j,1:kdm,1:ntracr)=trcraij(1:kdm,1:ntracr)  
      endif
      
#endif
      
   end do
   end do

   do k = 1, kdm
      print *,maxval(dpold(:,:,k)-dp(:,:,k))/onem
   end do

#if defined (ECO)     
   deallocate(trcraij,dpf,lcm,cfi,prsf,dpfor)
#endif





   ! Loop over restart file
   rstind=1 ! Restart index
   allok=.true.
   do while ( allok)

      ! Get header info from restart
      call rst_header_from_index(restart(1:fnd-1)//'.b', &
            cfld,vlevel,tlevel,rstind,bmin,bmax,.true.)

      allok=tlevel/=-1 ! test to see if read was ok

      print *,cfld



      if (allok) then 

         ! Get actual field  - for now we use the READRAW routine (later on we
         ! should switch to m_get_mod_fld
!         call READRAW(fldr4,amin,amax,idm,jdm,.false.,spval,restart(1:fnd-1)//'.a',rstind)
!         fld=fldr4

         ! Sjekk p at vi har lest rett - samanlign max/min fr filene
!         if     (abs(amin-bmin).gt.abs(bmin)*1.e-4 .or. &
!                 abs(bmax-amax).gt.abs(bmax)*1.e-4     ) then
!            print *,'Inconsistency between .a and .b files'
!            print *,'.a : ',amin,amax
!            print *,'.b : ',bmin,bmax
!            print *,cfld,vlevel,tlevel
!            call exit(1)
!         end if

	 call get_mod_fld_new(restart(1:fnd-1),fld(:,:),imem,cfld,vlevel,1,idm,jdm)

         if (trim(cfld)=='temp') then

            ! need salinity as well
            call get_mod_fld_new(restart(1:fnd-1),saln(:,:),imem,'saln    ',vlevel,1,idm,jdm)

            !if (tlevel==-1) then
            !   print *,'Could not get salinity field'
            !   call exit(1)
            !end if

            ! keep water warmer than freezing point
            do j=1,jdm
            do i=1,idm
               fld(i,j)=max(-.057*saln(i,j),fld(i,j))
               fld(i,j)=min(35.,fld(i,j)) !FC: cut off values that are too warm
            end do
            end do
         else if (trim(cfld)=='saln') then
           do j=1,jdm
           do i=1,idm
              fld(i,j)=max(5.,fld(i,j)) ! LB :no water fresher than 5 psu (Baltic)
              fld(i,j)=min(41.,fld(i,j))! FC: no water saltier than 40 psu
           end do
           end do
         else if (trim(cfld)=='pstarb') then
              fld(:,:) = (sqrt(PSTARB0 ** 2 + fld(:,:) ** 2) + fld(:,:)) / 2.0d0;
         else if (trim(cfld)=='dp') then
            fld = dp(:,:,vlevel) ! NB, one time level 
#if defined (ECO)  
         else if (cfld(1:6)=='tracer') then  	    
	    !updating the file!
	    ktrcr=tracr_get_incr(cfld(7:8))
	    if (ktrcr==-1)then
	      print*,'alert tracer unknow'
	      exit
	    endif
            fld(:,:)= tracerf(:,:,vlevel,ktrcr)
	    
#endif	 
	 end if ! No correction for other fields in the hycom restart file

         call put_mod_fld(trim(newfile),fld,imem,cfld,vlevel,tlevel,rstind,idm,jdm)


            
         rstind=rstind+1
      end if ! read of template was ok
   end do

#if defined (ECO)
   deallocate(tracerf)
#endif


   ! put_mod_fld does not include the header from the original file
   open(10,file=trim(newfile)//'.b',status='old')
   open(20,file=restart(1:fnd-1)//'.b',status='old')
   ! Supports parallel execution with different members
   open(11,file='tmp'//cmem//'.b',status='replace')

   ! Header from original file
   read(20,'(a80)') char80 ; write(11,'(a80)') char80
   read(20,'(a80)') char80 ; write(11,'(a80)') char80
   close(20)

   ! The rest from the newly created file
   ios=0
   do while(ios==0)
      read(10,'(a80)',iostat=ios) char80  
      if (ios==0) write(11,'(a80)') char80
   end do

   close(10)
   close(11)
   close(20)

   ! Move the new .b file to "newfile"
   !SERIAL call system('mv tmp'//cmem//'.b '//trim(newfile)//'.b')



!###################################################################
!####################### FIX   ICE   MODEL #########################
!###################################################################





   if (iargc()==3) then

      inquire(exist=ex,file=trim(icerestart))
      if (.not.ex) then
         print *,icerestart//' does not exist!'
         print *,'(fixhycom)'
         stop
      end if

      inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm  !,iceU,iceV
      open(10,file=icerestart,form='unformatted',access='direct',recl=reclICE)
      read(10,rec=imem)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
      close(10)

      ! PS 25/06/2010 max(ficem) = 0.995 - max the model allows
      ficem=min(max(0.,ficem), 0.995)
      hicem=min(max(0.,hicem),15.)
      hsnwm=min(max(0.,hsnwm), 4.)

      open (10,file='fix'//trim(icerestart),form='unformatted',access='direct',recl=reclICE)
      write(10,rec=imem)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
      close(10)
   end if
      




!KAL
!KAL   ! This is to make it easier for supporting programs
!KAL   open(10,file='file.EnKF_postprocess',status='replace')
!KAL   write(10,'(a)') 'analysis'//cmem
!KAL   close(10)
!KAL
!KAL
!KAL
!KAL   ! ice processing Loop over restart file
!KAL   print *
!KAL   print *
!KAL   print *,'processing ice restart file'
!KAL   rstind=1 ! Restart index
!KAL   allok=.true.
!KAL   call system("cp ensemble_TMP_ICE.uf ensemble_TMP_ICE_final.uf")
!KAL   do ifld=1,numfields
!KAL      cfld=  fieldnames (ifld)
!KAL      vlevel=fieldlevels(ifld)
!KAL      if (trim(cfld)=='icec' .or. trim(cfld)=='hice') then
!KAL         nomatch=.true.
!KAL         do iafile=1,nproc ! List of procs used in analysis
!KAL            write(cproc,'(i3.3)') iafile-1
!KAL
!KAL            ! NB - time level=1
!KAL            ! NB2 - the files dumped in the analysis lack a header (last argument
!KAL            ! is false)
!KAL            call rst_index_from_header(trim(afile)//'.b',cfld,vlevel,1, &
!KAL                                             tmpindx,bmin,bmax,.false.) 
!KAL
!KAL
!KAL            if (tmpindx/=-1) then
!KAL               print '(a8," -- layer:",i4,"  match : record, file",i4," ",a)', cfld, vlevel,tmpindx,trim(afile)
!KAL               nomatch=.false.
!KAL               exit
!KAL            end if
!KAL
!KAL            ! Read field from analysed file 
!KAL            call READRAW(fldr4,amin,amax,idm,jdm,.false.,spval,trim(afile)//'.a',tmpindx)
!KAL            fld=fldr4
!KAL
!KAL           ! Sjekk p at vi har lest rett - samanlign max/min fr filene
!KAL           if     (abs(amin-bmin).gt.abs(bmin)*1.e-4 .or. &
!KAL                   abs(bmax-amax).gt.abs(bmax)*1.e-4     ) then
!KAL              print *,'Inconsistency between .a and .b files'
!KAL              print *,'.a : ',amin,amax
!KAL              print *,'.b : ',bmin,bmax
!KAL              print *,cfld,vlevel,tlevel
!KAL              call exit(1)
!KAL           end if
!KAL
!KAL         end do
!KAL
!KAL
!KAL         ! Check if we got ice concentration or ice thickness
!KAL         if (.not.nomatch) then
!KAL
!KAL
!KAL            inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm  !,iceU,iceV
!KAL            open(10,file=trim(icetemplate),form='unformatted',access='direct',recl=reclICE,status='old')
!KAL            open(11,file='ensemble_TMP_ICE_final.uf',form='unformatted',access='direct',recl=reclICE,status='unknown')
!KAL            read(10,rec=imem,iostat=ios)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
!KAL            if (trim(cfld)=='icec') ficem=fld
!KAL            if (trim(cfld)=='hice') hicem=fld
!KAL            write(11,rec=imem,iostat=ios2)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
!KAL            close(10)
!KAL            close(11)
!KAL
!KAL            if (ios/=0 .or.ios2/=0) then
!KAL               print *,ios
!KAL               print *,ios2
!KAL               print *,'Error when writing to ice ens file'
!KAL               call exit(1)
!KAL            end if
!KAL         end if
!KAL      end if
!KAL   end do
!KAL
!KAL
!KAL   print *,'Normal exit of EnKF_postprocess'
!KAL   print *,'TODO: Process dp inconsistencies'
!KAL   call exit(0)
!KAL
!KAL
!KAL
end program