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

program EnKF_postprocess
! KAL -- The new EnKF is MPI-parallelized, and each thread will dump
! KAL -- to its own file named "analysisXXX_procXXX.[ab]
! KAL --
! KAL -- This routine will gather the
! KAL -- analysis from the separate analyzed files into one complete restart file.
! KAL -- To do this, a "template" restart file must be specified. This file
! KAL -- copies non-existing variables in the analyzed fields from the template,
! KAL -- and into the final analysis. The final files produced by this routine 
! KAL -- are named "analysisXXX.[ab]".
! KAL -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
! KAL -- Input arguments:
! KAL --     template restart file
! KAL --     template ice restart file
! KAL --     ensemble member
! KAL --     number of MPI threads used in analysis
! KAL -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 
! KAL --
! KAL -- NB: This routine will not check or modify the fields. To do that, 
! KAL --     use "consistency" (check) and "fixhycom" (fix)



   use mod_raw_io
   use m_parse_blkdat
   use m_put_mod_fld
   implicit none

   integer*4, external :: iargc

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

   character(len=80) :: afile
   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*4, allocatable:: fldr4(:,:)
   real*4 :: spval,amin,amax

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

   integer,parameter :: numfields=2
   integer :: ios,ios2, reclICE,ifld
   character(len=8) :: fieldnames(numfields)
   integer :: fieldlevels(numfields)


   if (iargc()==4) then
      call getarg(1,template)
      call getarg(2,icetemplate)
      call getarg(3,ctmp)
      read(ctmp,*) imem
      write(cmem,'(i3.3)') imem
      call getarg(4,ctmp)
      read(ctmp,*) nproc
   else
      print *,'usage: EnKF_postprocess restart_template ice_template ensemble_member nproc'
      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 (idm>0 .and. idm < 1e4 .and. jdm>0 .and. jdm<1e4) then
      allocate(fld  (idm,jdm))
      allocate(fldr4(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))
   else
      print *,'fld allocate error'
      stop '(EnKF_postprocess)'
   end if


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


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

   print *,template(1:fnd-1)//'.b'

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

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

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


      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,template(1:fnd-1)//'.a',rstind)
         fld=fldr4


         !print *,cfld,tlevel, vlevel


         ! From the template, we have the stuff we need, now we go looking
         ! in the EnKF analysis files output for the correct input. Note that
         ! all time levels are set equal in the restart. This may introduce 
         ! imbalances, but the alternative is probably worse.
         !
         ! KAL -- need list of analysis files here, for now we hardcode
         nomatch=.true.
         do iafile=1,nproc ! List of procs used in analysis

            write(cproc,'(i3.3)') iafile-1
            
            ! Temporary name, will change
            afile='analysis'//cmem//'_proc'//cproc


            ! NB - time level=1
            ! NB2 - the files dumped in the analysis lack a header (last argument
            ! is false)
            if (trim(cfld)=='ficem') then
              ctmp='icec'
            elseif (trim(cfld)=='hicem') then
              ctmp='hice'
            else
              ctmp=trim(cfld)
            endif
            call rst_index_from_header(trim(afile)//'.b',ctmp,vlevel,1, &
                                          tmpindx,bmin,bmax,.false.) 

            if (tmpindx/=-1) then
               !if (tlevel/=1) print *,'--> replacing time level with 1'
               print '(a8," -- time, layer:",2i4,"  match : record, file",i4," ",a)', cfld,tlevel, vlevel,tmpindx,trim(afile)
               nomatch=.false.

               ! Read field from analysed file 
               call READRAW(fldr4,amin,amax,idm,jdm,.false.,spval,trim(afile)//'.a',tmpindx)
               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



               ! put into final, processed file -- imem is not used, actually
               call put_mod_fld('analysis'//cmem,fld,imem,cfld,vlevel,tlevel,rstind,idm,jdm)
               

               exit
            end if
         end do


         if (nomatch) then
            print '(a8," -- time, layer:",2i4," - no match - replace with template")',cfld,tlevel,vlevel

            ! put template values into final, processed file -- imem is not used, actually
            call put_mod_fld('analysis'//cmem,fld,imem,cfld,vlevel,tlevel,rstind,idm,jdm)
         end if
            
         rstind=rstind+1
      end if ! read of template was ok
   end do

#if ! defined (SINGLE_RESTART)

   ! ice processing Loop over restart file
   print *
   print *
   print *,'processing ice restart file'
   rstind=1 ! Restart index
   allok=.true.

   ! Temporary solution
   fieldnames(1)='icec'
   fieldnames(2)='hice'
   vlevel=0

   ! Copy template record imem to analysysICE.uf
   inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm
   open(10,file=trim(icetemplate),form='unformatted',access='direct',recl=reclICE,status='old')
   read(10,rec=imem,iostat=ios)ficem,hicem,hsnwm,ticem,tsrfm
   close(10)

   open(11,file='analysisICE.uf',form='unformatted',access='direct',recl=reclICE,status='unknown')
   write(11,rec=imem,iostat=ios2)ficem,hicem,hsnwm,ticem,tsrfm
   close(11)

   do ifld=1,numfields
      cfld=  fieldnames (ifld)
      if (trim(cfld)=='icec' .or. trim(cfld)=='hice') then
         nomatch=.true.
         do iafile=1,nproc ! List of procs used in analysis
            write(cproc,'(i3.3)') iafile-1

            ! NB - time level=1
            ! NB2 - the files dumped in the analysis lack a header (last argument
            ! is false)
            call rst_index_from_header(trim(afile)//'.b',cfld,vlevel,1, &
                                             tmpindx,bmin,bmax,.false.) 


            if (tmpindx/=-1) then
               print '(a8," -- layer:",i4,"  match : record, file",i4," ",a)', cfld, vlevel,tmpindx,trim(afile)
               nomatch=.false.

               !print *,'Got match for '//cfld

               ! Read field from analysed file 
               call READRAW(fldr4,amin,amax,idm,jdm,.false.,spval,trim(afile)//'.a',tmpindx)
               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
              exit

            end if

         end do


         ! Check if we got ice concentration or ice thickness
         if (.not.nomatch) then


            inquire(iolength=reclICE)ficem,hicem,hsnwm,ticem,tsrfm  !,iceU,iceV
            open(11,file='analysisICE.uf',form='unformatted',access='direct',recl=reclICE,status='unknown')
            if (trim(cfld)=='icec') ficem=fld
            if (trim(cfld)=='hice') hicem=fld
            write(11,rec=imem,iostat=ios2)ficem,hicem,hsnwm,ticem,tsrfm !,iceU,iceV
            close(11)

            if (ios/=0 .or.ios2/=0) then
               print *,ios
               print *,ios2
               print *,'Error when writing to ice ens file'
               call exit(1)
            end if
         end if
      end if
   end do
#endif

   print *,'Normal exit of EnKF_postprocess'
   print *,'TODO: Process dp inconsistencies'
#if defined(AIX)
   call exit_(0)
#else
   call exit(0)
#endif



end program