regridStates.F90

Go to the documentation of this file.

 program regridstates

 use mpi_f08
 use esmf

 use grids_io, only     : setup_grid, &
                          write_from_fields, &
                          read_into_fields, &
                          n_tiles, &
                          grid_setup_type

 use utilities, only    : error_handler

 implicit none

 integer, parameter             :: max_vars = 10

 ! namelist inputs
 character(len=15)              :: variable_list(max_vars)
 integer                        :: n_vars, n_tims, extrap_levs
 integer                        :: time_list(10)
 logical                        :: add_time_dim
 real(esmf_kind_r8)             :: missing_value ! value given to unmapped cells in the output grid

 type(grid_setup_type)          :: grid_setup_in, grid_setup_out

 integer                        :: ierr, localpet, npets
 integer                        :: v, t, srcterm

 character(100)                 :: fname_time

 type(esmf_vm)                  :: vm
 type(esmf_grid), allocatable   :: grid_in(:)
 type(esmf_grid)                :: grid_out
 type(esmf_field), allocatable  :: fields_in(:,:)
 type(esmf_field), allocatable  :: fields_out(:,:)
 type(esmf_routehandle)         :: regrid_route
 real(esmf_kind_r8), pointer    :: ptr_out(:,:)

 integer :: ut

 real :: t1, t2, t3, t4
 character(len=3)               :: tstr

 ! see README for details of namelist variables.
 namelist /config/ n_vars, variable_list, missing_value, extrap_levs, time_list, add_time_dim

! INITIALIZE
!-------------------------------------------------------------------------

 call cpu_time(t1)

 call mpi_init(ierr)

 call esmf_initialize(rc=ierr)
 if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("INITIALIZING ESMF", ierr)

 call esmf_vmgetglobal(vm, rc=ierr)
 if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN VMGetGlobal", ierr)

 call esmf_vmget(vm, localpet=localpet, petcount=npets, rc=ierr)
 if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN VMGet", ierr)

!-------------------------------------------------------------------------
! RUN
!-------------------------------------------------------------------------

 print*,'** pets: local, total: ',localpet, npets

 ! checks

 if (mod(npets,n_tiles) /= 0) then
   call error_handler("must run with a task count that is a multiple of 6", 1)
 endif

!------------------------
! read in namelist

 ! defaults
 missing_value=-999.
 extrap_levs=2
 time_list=-1

 open(newunit=ut, file='regrid.nml', iostat=ierr)
 if (ierr /= 0) call error_handler("OPENING regrid NAMELIST.", ierr)
 read(ut, nml=config, iostat=ierr)
 if (ierr /= 0) call error_handler("OPENING config NAMELIST.", ierr)
 call readin_setup(ut,"input",grid_setup_in)
 call readin_setup(ut,"output",grid_setup_out)
 close (ut)

 
 n_tims = 0
 do t=1,10
   if (time_list(t) .lt. 0) exit
   n_tims = n_tims + 1
 enddo
 if (n_tims < 1) then
   call error_handler("n_tims < 1. must have at least one valid increment hour in time_list", 1)
 endif

!------------------------
! Create esmf grid objects for input and output grids, and add land masks

! TO DO - can we make the number of tasks more flexible for fv3

 if (localpet==0) print*,'** Setting up grids'
 allocate(grid_in(n_tims))
 do t = 1, n_tims
   if (grid_setup_in%mask_from_input) then
     call setup_grid(localpet, npets, grid_setup_in, grid_in(t), time_list(t) )
   else
     call setup_grid(localpet, npets, grid_setup_in, grid_in(t))
   endif
 enddo
 call setup_grid(localpet, npets, grid_setup_out, grid_out )

!------------------------
! Create input and output fields

 if (localpet==0) print*,'** Creating/Reading fields'

! input
 allocate(fields_in(n_tims,n_vars))

 do t = 1, n_tims
     do v = 1, n_vars

        fields_in(t,v)  = esmf_fieldcreate(grid_in(t), &
                            typekind=esmf_typekind_r8, &
                            staggerloc=esmf_staggerloc_center, &
                            name="input for regridding", &
                            rc=ierr)

        if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
           call error_handler("in FieldCreate "//trim(variable_list(v)), ierr)
     end do
 end do

! output
 allocate(fields_out(n_tims,n_vars))

 do t = 1, n_tims
     do v = 1, n_vars

         fields_out(t,v)  = esmf_fieldcreate(grid_out, &
                                           typekind=esmf_typekind_r8, &
                                           staggerloc=esmf_staggerloc_center, &
                                           name="output of regridding", &
                                           rc=ierr)
         if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
         call error_handler("in FieldCreate, field_out", ierr)


         ! set the default output value (for non-mapped cells)
         call esmf_fieldget(fields_out(t,v), &
                            farrayptr=ptr_out, &
                            rc=ierr)
         if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
         call error_handler("IN FieldGet", ierr)

         ptr_out=missing_value

     enddo
 enddo

!------------------------
! read data into input fields

 do t = 1, n_tims
        
        write(tstr,"(I3.3)")time_list(t)
        fname_time = trim(grid_setup_in%fname)//tstr//".nc"
        write(6,*) 'reading into ', trim(fname_time)
        call read_into_fields(localpet, grid_setup_in%ires, grid_setup_in%jres, &
                                 trim(fname_time), trim(grid_setup_in%dir), &
                                 grid_setup_in, n_vars, variable_list(1:n_vars), fields_in(t,:))
 enddo

 call cpu_time(t2)
!------------------------
! regrid the input fields to the output grid

 if (localpet==0) print*,'** Performing regridding'

 srcterm=1
 ! get regriding route for a field (only uses the grid info in the field)
 ! to turn off masking, remove [src/dstMaskVales] argumemnts
 call esmf_fieldregridstore(srcfield=fields_in(1,1), srcmaskvalues=(/0/), &
                            dstfield=fields_out(1,1), dstmaskvalues=(/0/), &
                            ! allow unmapped grid cells, without returning error
                            unmappedaction=esmf_unmappedaction_ignore, &
                            polemethod=esmf_polemethod_allavg, &
                            ! fill un-mapped grid cells with a neighbour
                            extrapmethod=esmf_extrapmethod_creep, &
                            ! number of "levels" of neighbours to search for a value
                            extrapnumlevels=extrap_levs, &
                            ! needed for reproducibility
                            ! (combined with ESMF_TERMORDER_SRCSEQ below)
                            srctermprocessing=srcterm, &
                            routehandle=regrid_route, &
                            ! use bilinear interp (slightly better results than PATCH)
                            regridmethod=esmf_regridmethod_bilinear, rc=ierr)
 if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldRegridStore", ierr)

! do the re-gridding

 call cpu_time(t3)

 do t=1, n_tims
     do v=1, n_vars
         call esmf_fieldregrid(fields_in(t,v), &
                               fields_out(t,v), &
                               routehandle=regrid_route, &
                               zeroregion=esmf_region_select, & ! initialize output with missing_value
                               termorderflag=esmf_termorder_srcseq, rc=ierr)
         if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
          call error_handler("IN FieldRegrid", ierr)
     enddo
 enddo

! TO-DO: terrain-correct temperatures (all layers?)

! write out fields on destination grid. All times into same file.

 if (localpet==0) print*,'** Writing out regridded fields'

 call write_from_fields(localpet, grid_setup_out%ires, grid_setup_out%jres,     &
                          trim(grid_setup_out%fname), trim(grid_setup_out%dir), &
                          n_vars, n_tims, variable_list(1:n_vars), fields_out, add_time_dim)


! clean up

 call esmf_fieldregridrelease(routehandle=regrid_route, rc=ierr)
 if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldRegridRelease", ierr)

 do t = 1, n_tims
     do v = 1, n_vars
         call esmf_fielddestroy(fields_in(t,v),rc=ierr)
         if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
                call error_handler("DESTROYING FIELD", ierr)

         call esmf_fielddestroy(fields_out(t,v),rc=ierr)
         if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
            call error_handler("DESTROYING FIELD", ierr)
     enddo
 
     call esmf_griddestroy(grid_in(t), rc=ierr)
     if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
        call error_handler("DESTROYING GRID", ierr)
 enddo

 call esmf_griddestroy(grid_out,rc=ierr)
 if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("DESTROYING GRID", ierr)

!-------------------------------------------------------------------------
! FINALIZE
!-------------------------------------------------------------------------

 call esmf_finalize(endflag=esmf_end_keepmpi, rc=ierr)
 call mpi_finalize(ierr)

 call cpu_time(t4)
 if (localpet==0) print*, '** time in tile2tile', t4 - t1
 if (localpet==0) print*, '** time in RegridStore', t3 - t2

 print*,"** DONE.", localpet

 end program regridstates