search_util.F90

Go to the documentation of this file.

 
 
 module search_util
 
 private
 
 public :: search
 
 contains
 
 subroutine search (field, mask, idim, jdim, tile, field_num, latitude, terrain_land, soilt_climo)
 
 use mpi_f08
 use esmf
 
 implicit none
 
 
 integer, intent(in)               :: idim, jdim, tile, field_num
 integer(esmf_kind_i8), intent(in) :: mask(idim,jdim)
 
 real(esmf_kind_r8), intent(in), optional :: latitude(idim,jdim)
 real(esmf_kind_r8), intent(in), optional :: terrain_land(idim,jdim)
 real(esmf_kind_r8), intent(in), optional :: soilt_climo(idim,jdim)
 
 real(esmf_kind_r8), intent(inout) :: field(idim,jdim)
 
 integer                           :: i, j, krad, ii, jj
 integer                           :: istart, iend
 integer                           :: jstart, jend
 integer                           :: ierr
 
 real                              :: default_value
 real(esmf_kind_r8)                :: field_save(idim,jdim)
 integer                           :: repl_nearby, repl_default
 
!-----------------------------------------------------------------------
! Set default value.
!-----------------------------------------------------------------------
 
 select case (field_num)
   case (0) ! most nst fields
     default_value = 0.0
   case (1) ! ifd
     default_value = 1.0
   case (7) ! terrain height, flag value to turn off terrain adjustment
            ! of soil temperatures.
     default_value = -99999.9
   case (11) ! water temperature will use latitude-dependent value
     default_value = -999.0
   case (21) ! ice temperature
     default_value = 265.0
   case (30) ! xz
     default_value = 30.0
   case (65) ! snow liq equivalent
     default_value = 0.0
   case (66) ! snow depth
     default_value = 0.0
   case (83) ! z0 (cm)
     default_value = 0.01
   case (85) ! soil temp
     default_value = 280.0
   case (86) ! soil moisture (volumetric)
     default_value = 0.18
   case (91) ! sea ice fraction
     default_value = 0.5
   case (92) ! sea ice depth (meters)
     default_value = 1.0
   case (223) ! canopy moist
     default_value = 0.0
   case (224) ! soil type, flag value to turn off soil moisture rescaling.
     default_value = -99999.9
   case (225) ! vegetation type, flag value to be replaced
     default_value = -99999.9
   case (226) ! vegetation fraction, flag value to be replaced
     default_value = 0.5
   case (227) ! max vegetation fraction, flag value to be replaced
     default_value = 0.5
   case (228) ! min vegetation fraction, flag value to be replaced
     default_value = 0.5
   case (229) ! lai, flag value to be replaced
     default_value = 1.0
   case (230) ! soil type on the input grid
     default_value = 11.0
   case (231) ! soil type at land ice points
     default_value = 16.0
   case default
     print*,'- FATAL ERROR.  UNIDENTIFIED FIELD NUMBER : ', field_num
     call mpi_abort(mpi_comm_world, 77, ierr)
 end select
 
!-----------------------------------------------------------------------
! Perform search and replace.
!-----------------------------------------------------------------------
 
 field_save = field
 repl_nearby = 0
 repl_default = 0
!$OMP PARALLEL DO DEFAULT(NONE), &
!$OMP SHARED(IDIM,JDIM,MASK,FIELD_SAVE,FIELD,TILE,LATITUDE,DEFAULT_VALUE,FIELD_NUM,REPL_NEARBY,REPL_DEFAULT,SOILT_CLIMO,TERRAIN_LAND), &
!$OMP PRIVATE(I,J,KRAD,ISTART,IEND,JSTART,JEND,II,JJ)
 
 j_loop : do j = 1, jdim
   i_loop : do i = 1, idim
 
     if (mask(i,j) == 1 .and. field_save(i,j) < -9999.0) then
 
       krad_loop : do krad = 1, 100
 
         istart = i - krad
         iend   = i + krad
         jstart = j - krad
         jend   = j + krad
 
         jj_loop : do jj = jstart, jend
         ii_loop : do ii = istart, iend
 
!-----------------------------------------------------------------------
!          Search only along outer square.
!-----------------------------------------------------------------------
 
           if ((jj == jstart) .or. (jj == jend) .or.   &
               (ii == istart) .or. (ii == iend)) then
 
             if (jj < 1 .or. jj > jdim) cycle jj_loop
             if (ii < 1 .or. ii > idim) cycle ii_loop
 
               if (mask(ii,jj) == 1  .and. field_save(ii,jj) > -9999.0) then
                 field(i,j) = field_save(ii,jj)
                ! write(6,100) field_num,tile,i,j,ii,jj,field(i,j)
                ! When using non-GFS data, there are a lot of these print statements even
                ! when everything is working correctly. Count instead of printing each
                 repl_nearby = repl_nearby + 1
                 cycle i_loop
               endif
 
           endif
 
         enddo ii_loop
         enddo jj_loop
 
       enddo krad_loop
 
       if (field_num == 11) then
         call sst_guess(latitude(i,j), field(i,j))
       elseif (field_num == 91) then  ! sea ice fract
         if (abs(latitude(i,j)) > 55.0) then
           field(i,j) = default_value
           repl_default = repl_default + 1
         else
           field(i,j) = 0.0
           repl_default = repl_default + 1
         endif
       elseif (field_num == 7 .and. PRESENT(terrain_land)) then 
         !  Terrain heights for isolated landice points never get a correct value, so replace
         !  with terrain height from the input grid interpolated to the target grid
         field(i,j) = terrain_land(i,j)
         repl_default = repl_default + 1
       elseif (field_num == 224 .and. PRESENT(soilt_climo)) then
          !  When using input soil type fields instead of climatological data on the
          ! target grid, isolated land locations that exist in the target grid but
          ! not the input grid don't receiving proper soil type information, so replace
          ! with climatological values
         field(i,j) = soilt_climo(i,j)
         repl_default = repl_default + 1
       else
         field(i,j) = default_value  
         repl_default = repl_default + 1
       endif
 
       !write(6,101) field_num,tile,i,j,field(i,j)
 
     endif
   enddo i_loop
 enddo j_loop
!$OMP END PARALLEL DO
 
! 100 format(1x,"- MISSING POINT FIELD ",i4," TILE: ",i2," I/J: ",i5,i5," SET TO VALUE AT: ",i5,i5,". NEW VALUE IS: ",f8.3)
! 101 format(1x,"- MISSING POINT FIELD ",i4," TILE: ",i2," I/J: ",i5,i5," SET TO DEFAULT VALUE OF: ",f8.3)
 print*, "- TOTAL POINTS FOR VAR ", field_num, " REPLACED BY NEARBY VALUES: ", repl_nearby
 print*, "- TOTAL POINTS FOR VAR ", field_num, " REPLACED BY DEFAULT VALUE: ", repl_default
 
 end subroutine search
 
 subroutine sst_guess(latitude, sst)
 
 use esmf
 
 implicit none
 
 real(esmf_kind_r8), parameter :: SST_POLAR_IN_KELVIN = 273.16
 real(esmf_kind_r8), parameter :: SST_TROPICAL_IN_KELVIN = 300.0
 real(esmf_kind_r8), parameter :: POLAR_LATITUDE = 60.0
 real(esmf_kind_r8), parameter :: TROPICAL_LATITUDE = 30.0
 real(esmf_kind_r8), parameter :: DSST_DLAT = -0.8947
 real(esmf_kind_r8), parameter :: SST_Y_INTERCEPT = 326.84
 
 real(esmf_kind_r8), intent(in)  :: latitude
 
 real(esmf_kind_r8), intent(out) :: sst
 
 if (abs(latitude) >= polar_latitude) then
   sst = sst_polar_in_kelvin
 elseif (abs(latitude) <= tropical_latitude) then
   sst = sst_tropical_in_kelvin
 else
   sst = dsst_dlat * abs(latitude) + sst_y_intercept
 endif
 
 end subroutine sst_guess
 
 end module search_util