interp.F90

Go to the documentation of this file.

 subroutine interp(localpet, method, input_file)

 use esmf
 use netcdf
 use model_grid
 use source_grid
 use utils
 use mpi

 implicit none

 character(len=*), intent(in)       :: input_file

 integer                            :: rc, localpet
 integer                            :: i, j, ij, tile, n, ncid, status
 integer                            :: l(1), u(1), t
 integer                            :: clb_mdl(2), cub_mdl(2)
 integer                            :: varid, record
 integer                            :: tile_num, pt_loc_this_tile
 integer                            :: isrctermprocessing

 integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:)
 integer(esmf_kind_i4), pointer     :: unmapped_ptr(:)

 real(esmf_kind_r4), pointer        :: data_mdl_ptr(:,:)
 real(esmf_kind_r4), pointer        :: data_src_ptr(:,:)
 real(esmf_kind_r4), allocatable    :: data_src_global(:,:)
 real(esmf_kind_r4), allocatable    :: data_mdl_one_tile(:,:)
 real(esmf_kind_r4), allocatable    :: vegt_mdl_one_tile(:,:)
 real(esmf_kind_r4), allocatable    :: lat_mdl_one_tile(:,:)
 real(esmf_kind_r4), allocatable    :: lon_mdl_one_tile(:,:)

 type(esmf_regridmethod_flag),intent(in) :: method
 type(esmf_field)                        :: data_field_src
 type(esmf_routehandle)                  :: regrid_data
 type(esmf_polemethod_flag)              :: pole
 
 print*,"- CALL FieldCreate FOR SOURCE GRID DATA."
 data_field_src = esmf_fieldcreate(grid_src, &
                                  typekind=esmf_typekind_r4, &
                                  indexflag=esmf_index_global, &
                                  staggerloc=esmf_staggerloc_center, &
                                  name="source data", &
                                  rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldCreate", rc)

 print*,"- CALL FieldGet FOR SOURCE GRID DATA."
 nullify(data_src_ptr)
 call esmf_fieldget(data_field_src, &
                    farrayptr=data_src_ptr,  &
                    rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldGet", rc)

 print*,"- CALL FieldGet FOR MODEL GRID DATA."
 nullify(data_mdl_ptr)
 call esmf_fieldget(data_field_mdl, &
                    farrayptr=data_mdl_ptr,  &
                    computationallbound=clb_mdl, &
                    computationalubound=cub_mdl, &
                    rc=rc)
 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN FieldGet", rc)

 if (localpet == 0) then
   allocate(data_src_global(i_src,j_src))
 else
   allocate(data_src_global(0,0))
 endif

 print*,'- OPEN SOURCE FILE ', trim(input_file)
 status = nf90_open(input_file, nf90_nowrite, ncid)
 call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE")

 if (localpet == 0) then
   allocate(data_mdl_one_tile(i_mdl,j_mdl))
   allocate(mask_mdl_one_tile(i_mdl,j_mdl))
   allocate(lat_mdl_one_tile(i_mdl,j_mdl))
   allocate(lon_mdl_one_tile(i_mdl,j_mdl))
 else
   allocate(data_mdl_one_tile(0,0))
   allocate(mask_mdl_one_tile(0,0))
   allocate(lat_mdl_one_tile(0,0))
   allocate(lon_mdl_one_tile(0,0))
 endif

 record = 0

 time_loop : do t = 1, num_time_recs ! loop over each time period

 field_loop : do n = 1, num_fields ! loop over each surface field.

   record = record + 1

   if (localpet == 0) then
     status = nf90_inq_varid(ncid, field_names(n), varid)
     call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID")
     status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,t/), count=(/i_src,j_src,1/))
     call netcdf_err(status, "IN ROUTINE INTERP READING FIELD")
   endif

   print*,"- CALL FieldScatter FOR SOURCE GRID DATA."
   call esmf_fieldscatter(data_field_src, data_src_global, rootpet=0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldScatter.", rc)

   isrctermprocessing = 1

   if (record == 1) then

     if (method == esmf_regridmethod_bilinear) then
       pole = esmf_polemethod_allavg
     else
       pole = esmf_polemethod_none
     endif

     print*,"- CALL FieldRegridStore."
     nullify(unmapped_ptr)
     call esmf_fieldregridstore(data_field_src, &
                            data_field_mdl, &
                            srcmaskvalues=(/0/), &
                            dstmaskvalues=(/0/), &
                            polemethod=pole, &                     
                            unmappedaction=esmf_unmappedaction_ignore, &
                            normtype=esmf_normtype_fracarea, &
                            srctermprocessing=isrctermprocessing, &
                            routehandle=regrid_data, &
                            regridmethod=method, &
                            unmappeddstlist=unmapped_ptr,  &
                            rc=rc)
     if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
        call error_handler("IN FieldRegridStore.", rc)

   endif

   print*,"- CALL Field_Regrid."
   call esmf_fieldregrid(data_field_src, &
                         data_field_mdl, &
                         routehandle=regrid_data, &
                         termorderflag=esmf_termorder_srcseq, &
                         rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
      call error_handler("IN FieldRegrid.", rc)
 
!-----------------------------------------------------------------------
! Unmapped points are stored in "unmapped_ptr".  The pointer contains
! "ij" global indices as follows:
!
! tile 1: 1 thru (itile*jtile)
! tile n: (n-1)*(itile*jtile) thru n*(itile*jtile)
! 
! This "ij" index is converted to the tile number and i/j index of the
! field object.  This logic assumes the model grid object was
! created using "GLOBAL" indices. 
!
! Unmapped data points are given the flag value of -9999.9, which
! will be replaced in routine "search".
!-----------------------------------------------------------------------

   l = lbound(unmapped_ptr)
   u = ubound(unmapped_ptr)
   do ij = l(1), u(1)

     tile_num = ((unmapped_ptr(ij)-1) / (i_mdl*j_mdl)) ! tile number minus 1
     pt_loc_this_tile = unmapped_ptr(ij) - (tile_num * i_mdl * j_mdl) 
                                                  ! "ij" location of point within tile.

     j = (pt_loc_this_tile - 1) / i_mdl + 1
     i = mod(pt_loc_this_tile, i_mdl)
     if (i==0) i = i_mdl
     data_mdl_ptr(i,j) = -9999.9 
                                 
   enddo

! These fields are adjusted at landice.

   select case (trim(field_names(n)))
     case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type')
     if (localpet == 0) then
       allocate(vegt_mdl_one_tile(i_mdl,j_mdl))
     else
       allocate(vegt_mdl_one_tile(0,0))
     endif
   end select

   output_loop : do tile = 1, num_tiles

     print*,"- CALL FieldGather FOR MODEL LATITUDE."
     call esmf_fieldgather(latitude_field_mdl, lat_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
     if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
        call error_handler("IN FieldGather.", rc)

     print*,"- CALL FieldGather FOR MODEL LONGITUDE."
     call esmf_fieldgather(longitude_field_mdl, lon_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
     if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
        call error_handler("IN FieldGather.", rc)

     print*,"- CALL FieldGather FOR MODEL GRID DATA."
     call esmf_fieldgather(data_field_mdl, data_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
     if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
        call error_handler("IN FieldGather.", rc)

     print*,"- CALL FieldGather FOR MODEL GRID LAND MASK."
     call esmf_fieldgather(mask_field_mdl, mask_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
     if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
        call error_handler("IN FieldGather.", rc)

     select case (trim(field_names(n)))
       case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type')
         print*,"- CALL FieldGather FOR MODEL GRID VEG TYPE."
         call esmf_fieldgather(vegt_field_mdl, vegt_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
         if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
            call error_handler("IN FieldGather.", rc)
     end select

     if (localpet == 0) then
       print*,'- CALL SEARCH FOR TILE ',tile
       call search (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, tile, field_names(n))
       select case (field_names(n))
         case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type')
           call adjust_for_landice (data_mdl_one_tile, vegt_mdl_one_tile, i_mdl, j_mdl, field_names(n))
       end select
       where(mask_mdl_one_tile == 0) data_mdl_one_tile = missing
       call output (data_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, tile, record, t, n)
     endif

     if (field_names(n) == 'vegetation_type') then
       print*,"- CALL FieldScatter FOR MODEL GRID VEGETATION TYPE."
       call esmf_fieldscatter(vegt_field_mdl, data_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
       if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
          call error_handler("IN FieldScatter.", rc)
     endif

   enddo output_loop

   if (allocated(vegt_mdl_one_tile)) deallocate(vegt_mdl_one_tile)

 enddo field_loop
 enddo time_loop

 status=nf90_close(ncid)

 deallocate(data_mdl_one_tile, mask_mdl_one_tile)
 deallocate(data_src_global, lat_mdl_one_tile, lon_mdl_one_tile)

 print*,"- CALL FieldRegridRelease."
 call esmf_fieldregridrelease(routehandle=regrid_data, rc=rc)

 print*,"- CALL FieldDestroy."
 call esmf_fielddestroy(data_field_src, rc=rc)

 end subroutine interp

 subroutine adjust_for_landice(field, vegt, idim, jdim, field_ch)

 use esmf
 use mpi

 implicit none

 character(len=*), intent(in)      :: field_ch

 integer, intent(in)               :: idim, jdim

 real(esmf_kind_i4), intent(in)    :: vegt(idim,jdim)
 real(esmf_kind_r4), intent(inout) :: field(idim,jdim)

 integer, parameter                :: landice=15

 integer                           :: i, j, ierr

 real                              :: landice_value

 select case (field_ch)
   case ('substrate_temperature') ! soil substrate temp
     landice_value = 273.15
     do j = 1, jdim
     do i = 1, idim
       if (nint(vegt(i,j)) == landice) then
         field(i,j) = min(field(i,j), landice_value)
       endif
     enddo
     enddo
   case ('vegetation_greenness') ! vegetation greenness
     landice_value = 0.01 ! 1.0% is bare ground
     do j = 1, jdim
     do i = 1, idim
       if (nint(vegt(i,j)) == landice) then
         field(i,j) = landice_value
       endif
     enddo
     enddo
   case ('leaf_area_index') ! leaf area index
     landice_value = 0.0 ! bare ground
     do j = 1, jdim
     do i = 1, idim
       if (nint(vegt(i,j)) == landice) then
         field(i,j) = landice_value
       endif
     enddo
     enddo
   case ('slope_type') ! slope type
     landice_value = 9.0
     do j = 1, jdim
     do i = 1, idim
       if (nint(vegt(i,j)) == landice) then
         field(i,j) = landice_value
       else
         if (nint(field(i,j)) == nint(landice_value)) field(i,j) = 2.0
       endif
     enddo
     enddo
   case ('soil_type') ! soil type
     landice_value = 16.0
     do j = 1, jdim
     do i = 1, idim
       if (nint(vegt(i,j)) == landice) then
         field(i,j) = landice_value
       else
         if (nint(field(i,j)) == nint(landice_value)) field(i,j) = 6.0
       endif
     enddo
     enddo
   case default
     print*,'- FATAL ERROR IN ROUTINE ADJUST_FOR_LANDICE.  UNIDENTIFIED FIELD : ', field_ch
     call mpi_abort(mpi_comm_world, 57, ierr)
 end select

 end subroutine adjust_for_landice