program_setup.F90

Go to the documentation of this file.

 
 
 module program_setup
 
 use esmf
 use utilities, only                    : error_handler, to_lower
 
 implicit none
 
 private
 
 character(len=500), public      :: varmap_file = "NULL"
 character(len=500), public      :: atm_files_input_grid(6) = "NULL"
 character(len=500), public      :: atm_core_files_input_grid(7) = "NULL"
 character(len=500), public      :: atm_tracer_files_input_grid(6) = "NULL"
 character(len=500), public      :: data_dir_input_grid = "NULL"
 character(len=500), public      :: fix_dir_target_grid = "NULL"
 character(len=500), public      :: mosaic_file_input_grid = "NULL"
 character(len=500), public      :: mosaic_file_target_grid = "NULL"
 character(len=500), public      :: nst_files_input_grid = "NULL"
 character(len=500), public      :: grib2_file_input_grid = "NULL"
 character(len=500), public      :: geogrid_file_input_grid = "NULL"
 character(len=500), public      :: orog_dir_input_grid = "NULL"
 character(len=500), public      :: orog_files_input_grid(6) = "NULL"
 character(len=500), public      :: orog_dir_target_grid = "NULL"
 character(len=500), public      :: orog_files_target_grid(6) = "NULL"
 character(len=500), public      :: sfc_files_input_grid(6) = "NULL"
 character(len=500), public      :: vcoord_file_target_grid = "NULL"
 character(len=500), public      :: thomp_mp_climo_file= "NULL"
 character(len=15),  public      :: cres_target_grid = "NULL"
 character(len=500), public      :: atm_weight_file="NULL"
 character(len=25),  public      :: input_type="restart"
 character(len=20),  public      :: external_model="GFS"
 
 integer, parameter, public      :: max_tracers=100
 integer, public                 :: num_tracers
 integer, public                 :: num_tracers_input
 
 logical, allocatable, public    :: read_from_input(:)
 
 character(len=20), public       :: tracers(max_tracers)="NULL"
 character(len=20), public       :: tracers_input(max_tracers)="NULL"
 character(len=20), allocatable, public      :: missing_var_methods(:)
 character(len=20), allocatable, public      :: chgres_var_names(:)
 character(len=20), allocatable, public      :: field_var_names(:)
 
 character(len=500), public                  :: wam_parm_file="msis21.parm"
 
 integer, public                 :: cycle_year = -999
 integer, public                 :: cycle_mon = -999
 integer, public                 :: cycle_day = -999
 integer, public                 :: cycle_hour = -999
 integer, public                 :: regional = 0
 integer, public                 :: halo_bndy = 0
 integer, public                 :: halo_blend = 0
 integer, public                 :: nsoill_out = 4
 
 logical, public                 :: convert_atm = .false. 
 logical, public                 :: convert_nst = .false. 
 logical, public                 :: convert_sfc = .false. 
 logical, public                 :: wam_cold_start = .false. 
 
 ! Options for replacing vegetation/soil type, veg fraction, and lai with data from the grib2 file
 ! Default is to use climatology instead
 logical, public                 :: vgtyp_from_climo = .true. 
 logical, public                 :: sotyp_from_climo = .true. 
 logical, public                 :: vgfrc_from_climo = .true. 
 
 logical, public                 :: minmax_vgfrc_from_climo = .true. 
 logical, public                 :: lai_from_climo = .true. 
 logical, public                 :: tg3_from_soil = .false. 
 logical, public                 :: use_thomp_mp_climo=.false. 
 
 real, allocatable, public       :: drysmc_input(:)
 real, allocatable, public       :: drysmc_target(:)
 real, allocatable, public       :: maxsmc_input(:)
 real, allocatable, public       :: maxsmc_target(:)
 real, allocatable, public       :: refsmc_input(:)
 real, allocatable, public       :: refsmc_target(:)
 real, allocatable, public       :: wltsmc_input(:)
 real, allocatable, public       :: wltsmc_target(:)
 real, allocatable, public       :: bb_target(:)
 real, allocatable, public       :: satpsi_target(:)
 real(kind=esmf_kind_r4), allocatable, public  :: missing_var_values(:) 
 
 public :: read_setup_namelist
 public :: calc_soil_params_driver
 public :: read_varmap
 public :: get_var_cond
 
 contains
 
 subroutine read_setup_namelist(filename)
 
 use mpi_f08
 
 implicit none
 
 character(len=*), intent(in), optional :: filename
 character(len=250), allocatable :: filename_to_use
 
 integer                     :: i, is, ie, ierr, myrank
 
 namelist /config/ varmap_file, &
                   mosaic_file_target_grid, &
                   fix_dir_target_grid,     &
                   orog_dir_target_grid,    &
                   orog_files_target_grid,  &
                   mosaic_file_input_grid,  &
                   orog_dir_input_grid,     &
                   orog_files_input_grid,   &
                   nst_files_input_grid,    &
                   sfc_files_input_grid,    &
                   atm_files_input_grid,    &
                   atm_core_files_input_grid,    &
                   atm_tracer_files_input_grid,    &
                   grib2_file_input_grid, &
                   geogrid_file_input_grid, &
                   data_dir_input_grid,     &
                   vcoord_file_target_grid, &
                   cycle_year, cycle_mon, cycle_day,    &
                   cycle_hour, convert_atm, &
                   convert_nst, convert_sfc, &
                   wam_cold_start, &
                   vgtyp_from_climo, &
                   sotyp_from_climo, &
                   vgfrc_from_climo, &
                   minmax_vgfrc_from_climo, &
                   lai_from_climo, tg3_from_soil, &
                   regional, input_type, &
                   external_model, &
                   wam_parm_file, &
                   atm_weight_file, tracers, &
                   tracers_input, &
                   halo_bndy, & 
                   halo_blend, &
                   nsoill_out, &
                   thomp_mp_climo_file
 
 if (present(filename)) then
   filename_to_use = filename
 else
   filename_to_use = "./fort.41"
 endif
 
 call mpi_comm_rank(mpi_comm_world, myrank, ierr)
 
 if (myrank == 0) then
   print*,"- READ SETUP NAMELIST"
   open(41, file=trim(filename_to_use), iostat=ierr)
   if (ierr /= 0) call error_handler("OPENING SETUP NAMELIST.", ierr)
   read(41, nml=config, iostat=ierr)
   if (ierr /= 0) call error_handler("READING SETUP NAMELIST.", ierr)
   close (41)
 endif
 
 call mpi_barrier(mpi_comm_world,ierr)
 call mpi_bcast(varmap_file,len(varmap_file),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(mosaic_file_target_grid,len(mosaic_file_target_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(fix_dir_target_grid,len(fix_dir_target_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(orog_dir_target_grid,len(orog_dir_target_grid),mpi_character,0,mpi_comm_world,ierr)
 do i = 1, 6
   call mpi_bcast(orog_files_target_grid(i),len(orog_files_target_grid),mpi_character,0,mpi_comm_world,ierr)
 enddo
 call mpi_bcast(mosaic_file_input_grid,len(mosaic_file_input_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(orog_dir_input_grid,len(orog_dir_input_grid),mpi_character,0,mpi_comm_world,ierr)
 do i = 1, 6
   call mpi_bcast(orog_files_input_grid(i),len(orog_files_input_grid),mpi_character,0,mpi_comm_world,ierr)
 enddo
 call mpi_bcast(nst_files_input_grid,len(nst_files_input_grid),mpi_character,0,mpi_comm_world,ierr)
 do i = 1, 6
   call mpi_bcast(sfc_files_input_grid(i),len(sfc_files_input_grid),mpi_character,0,mpi_comm_world,ierr)
   call mpi_bcast(atm_files_input_grid(i),len(atm_files_input_grid),mpi_character,0,mpi_comm_world,ierr)
   call mpi_bcast(atm_tracer_files_input_grid(i),len(atm_tracer_files_input_grid),mpi_character,0,mpi_comm_world,ierr)
 enddo
 do i = 1, 7
   call mpi_bcast(atm_core_files_input_grid(i),len(atm_core_files_input_grid),mpi_character,0,mpi_comm_world,ierr)
 enddo
 call mpi_bcast(grib2_file_input_grid,len(grib2_file_input_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(geogrid_file_input_grid,len(geogrid_file_input_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(data_dir_input_grid,len(data_dir_input_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(vcoord_file_target_grid,len(vcoord_file_target_grid),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(cycle_year,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(cycle_mon,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(cycle_day,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(cycle_hour,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(convert_atm,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(convert_nst,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(convert_sfc,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(wam_cold_start,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(vgtyp_from_climo,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(sotyp_from_climo,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(vgfrc_from_climo,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(minmax_vgfrc_from_climo,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(lai_from_climo,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(tg3_from_soil,1,mpi_logical,0,mpi_comm_world,ierr)
 call mpi_bcast(regional,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(input_type,len(input_type),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(external_model,len(external_model),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(wam_parm_file,len(wam_parm_file),mpi_character,0,mpi_comm_world,ierr)
 call mpi_bcast(atm_weight_file,len(atm_weight_file),mpi_character,0,mpi_comm_world,ierr)
 do i = 1, max_tracers
   call mpi_bcast(tracers(i),len(tracers),mpi_character,0,mpi_comm_world,ierr)
   call mpi_bcast(tracers_input(i),len(tracers_input),mpi_character,0,mpi_comm_world,ierr)
 enddo
 call mpi_bcast(halo_bndy,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(halo_blend,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(nsoill_out,1,mpi_integer,0,mpi_comm_world,ierr)
 call mpi_bcast(thomp_mp_climo_file,len(thomp_mp_climo_file),mpi_character,0,mpi_comm_world,ierr)
 
 call to_lower(input_type)
 
 orog_dir_target_grid = trim(orog_dir_target_grid) // '/'
 orog_dir_input_grid = trim(orog_dir_input_grid) // '/'
 
!-------------------------------------------------------------------------
! Determine CRES of target grid from the name of the orography file.
!-------------------------------------------------------------------------
 
 is = 1
 ie = index(orog_files_target_grid(1), "_oro_") - 1
 
 if (ie == 0) then
   call error_handler("CANT DETERMINE CRES FROM OROG FILE.", 1)
 endif
 
 cres_target_grid = orog_files_target_grid(1)(is:ie)
 
 if (.not. convert_sfc .and. .not. convert_atm) then
   call error_handler("MUST CONVERT EITHER AN ATM OR SFC FILE.", 1)
 endif
 
!-------------------------------------------------------------------------
! Flag for processing stand-alone regional grid.  When '1', 
! remove halo from atmospheric and surface data and output
! atmospheric lateral boundary condition file. When '2',
! create lateral boundary file only.  When '0' (the default),
! process normally as a global grid.
!-------------------------------------------------------------------------
 
 if (regional > 0) then
   print*,"- PROCESSING A REGIONAL NEST WITH A BOUNDARY HALO OF ",halo_bndy
   print*,"- PROCESSING A REGIONAL NEST WITH A BLENDING HALO OF ",halo_blend
 else
   halo_bndy = 0
   halo_blend = 0
 endif
 
 num_tracers = 0
 do is = 1, max_tracers
   if (trim(tracers(is)) == "NULL") exit
   num_tracers = num_tracers + 1
   print*,"- TRACER NAME IN OUTPUT FILE ", trim(tracers(is))
 enddo
 
 num_tracers_input = 0
 do is = 1, max_tracers
   if (trim(tracers_input(is)) == "NULL") exit
   num_tracers_input = num_tracers_input + 1
   print*,"- TRACER NAME IN INPUT FILE ", trim(tracers_input(is))
 enddo
 
!-------------------------------------------------------------------------
! Ensure spo, spo2, and spo3 in tracers list if wam ic is on
!-------------------------------------------------------------------------
 
 if( wam_cold_start ) then
    ierr=3
    do is = 1, num_tracers
      if(trim(tracers(is)) == "spo"  ) ierr = ierr - 1
      if(trim(tracers(is)) == "spo2" ) ierr = ierr - 1
      if(trim(tracers(is)) == "spo3" ) ierr = ierr - 1
    enddo
    if (ierr /= 0) then
      print*,"-ERROR: spo, spo2, and spo3 should be in tracers namelist"
      call error_handler("WAM TRACER NAMELIST.", ierr)
    endif
    print*,"- WAM COLDSTART OPTION IS TURNED ON."
 endif  
 
!-------------------------------------------------------------------------
! Ensure program recognizes the input data type.  
!-------------------------------------------------------------------------
 
 select case (trim(input_type))
   case ("restart")
     print*,'- INPUT DATA FROM FV3 TILED RESTART FILES.'
   case ("history")
     print*,'- INPUT DATA FROM FV3 TILED HISTORY FILES.'
#ifdef CHGRES_ALL
   case ("gaussian_nemsio")
     print*,'- INPUT DATA FROM FV3 GAUSSIAN NEMSIO FILE.'
   case ("gfs_gaussian_nemsio")
     print*,'- INPUT DATA FROM SPECTRAL GFS GAUSSIAN NEMSIO FILE.'
   case ("gfs_sigio")
     print*,'- INPUT DATA FROM SPECTRAL GFS SIGIO/SFCIO FILE.'
#endif
   case ("gaussian_netcdf")
     print*,'- INPUT DATA FROM FV3 GAUSSIAN NETCDF FILE.'
   case ("grib2")
     print*,'- INPUT DATA FROM A GRIB2 FILE'
   case default
     call error_handler("UNRECOGNIZED INPUT DATA TYPE.", 1)
 end select
 
!-------------------------------------------------------------------------
! Ensure proper file variable provided for grib2 input  
!-------------------------------------------------------------------------
 
 if (trim(input_type) == "grib2") then
   if (trim(grib2_file_input_grid) == "NULL" .or. trim(grib2_file_input_grid) == "") then
     call error_handler("FOR GRIB2 DATA, PLEASE PROVIDE GRIB2_FILE_INPUT_GRID", 1)
   endif
 endif
 
 !-------------------------------------------------------------------------
! For grib2 input, warn about possibly unsupported external model types
!-------------------------------------------------------------------------
 
 if (trim(input_type) == "grib2") then
   if (.not. any((/character(4)::"GFS","NAM","RAP","HRRR","RRFS"/)==trim(external_model))) then
     call error_handler( "KNOWN SUPPORTED external_model INPUTS ARE GFS, NAM, RAP, HRRR, AND RRFS. " // &
    "IF YOU WISH TO PROCESS GRIB2 DATA FROM ANOTHER MODEL, YOU MAY ATTEMPT TO DO SO AT YOUR OWN RISK. " // &
    "ONE WAY TO DO THIS IS PROVIDE NAM FOR external_model AS IT IS A RELATIVELY STRAIGHT-" // &
    "FORWARD REGIONAL GRIB2 FILE. YOU MAY ALSO COMMENT OUT THIS ERROR MESSAGE IN " // &
    "program_setup.f90 LINE 389. NO GUARANTEE IS PROVIDED THAT THE CODE WILL WORK OR "// &
    "THAT THE RESULTING DATA WILL BE CORRECT OR WORK WITH THE ATMOSPHERIC MODEL.", 1)
   endif
 endif
 
!-------------------------------------------------------------------------
! For grib2 hrrr input without geogrid file input, warn that soil moisture interpolation
! will be less accurate
!-------------------------------------------------------------------------
 
 if (trim(input_type) == "grib2" .and. trim(external_model)=="HRRR") then
   if (trim(geogrid_file_input_grid) == "NULL" .or. trim(grib2_file_input_grid) == "") then
     print*, "HRRR DATA DOES NOT CONTAIN SOIL TYPE INFORMATION. WITHOUT"
     print*, "GEOGRID_FILE_INPUT_GRID SPECIFIED, SOIL MOISTURE INTERPOLATION MAY BE LESS ACCURATE."
   endif
 endif
 
 if (trim(thomp_mp_climo_file) /= "NULL") then
   use_thomp_mp_climo=.true.
   print*,"- WILL PROCESS CLIMO THOMPSON MP TRACERS FROM FILE: ", trim(thomp_mp_climo_file)
 endif
 
 return
 
 end subroutine read_setup_namelist
 
subroutine read_varmap
 
 implicit none
 
 integer                        :: istat, k, nvars, rc, localpet
 integer                        :: idum1(1), idum2(2)
 character(len=500)             :: line
 character(len=20),allocatable  :: var_type(:)
 
 type(esmf_vm)                  :: vm
 
 if (trim(input_type) == "grib2") then 
 
   call esmf_vmgetglobal(vm, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN VMGetGlobal", rc)
 
   call esmf_vmget(vm, localpet=localpet, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
    call error_handler("IN VMGet", rc)
 
   if (localpet == 0) then
 
     print*,"OPEN VARIABLE MAPPING FILE: ", trim(varmap_file)
     open(14, file=trim(varmap_file), form='formatted', iostat=istat)
     if (istat /= 0) then
       call error_handler("OPENING VARIABLE MAPPING FILE", istat)
     endif
 
     nvars = 0
 
!  Loop over lines of file to count the number of variables
     do
       read(14, '(A)', iostat=istat) line
       if (istat/=0) exit
       if ( trim(line) .eq. '' ) cycle
       nvars = nvars+1
     enddo
     if (nvars == 0) call error_handler("VARMAP FILE IS EMPTY.", -1)
 
     idum1(1) = nvars
 
   endif  ! localpet == 0
 
   call esmf_vmbroadcast(vm, idum1, 1, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
   nvars = idum1(1)
 
   allocate(chgres_var_names(nvars))
   allocate(field_var_names(nvars))
   allocate(missing_var_methods(nvars))
   allocate(missing_var_values(nvars))
   allocate(read_from_input(nvars))
 
   read_from_input(:) = .true.
 
   if (localpet == 0) then
 
     num_tracers_input = 0
     allocate(var_type(nvars))
     rewind(14)
 
     do k = 1,nvars
       read(14, *, iostat=istat) chgres_var_names(k), field_var_names(k) , &
                            missing_var_methods(k), missing_var_values(k), var_type(k)
       if (istat /= 0) call error_handler("READING VARIABLE MAPPING FILE", istat)
 
       if(trim(var_type(k))=='T') then
         num_tracers_input = num_tracers_input + 1
         tracers_input(num_tracers_input)=chgres_var_names(k)
         if ((trim(chgres_var_names(k)) == "ice_aero" .or. trim(chgres_var_names(k)) == "liq_aero") .and. &
           trim(thomp_mp_climo_file) .ne. "NULL" .and. trim(input_type) == "grib2") then
           call error_handler("VARMAP TABLE CONTAINS TRACER ENTRIES FOR THOMPSON AEROSOLS liq_aero or "// &
           "ice_aero. REMOVE THESE ENTRIES OR REMOVE THE NAMELIST ENTRY FOR "// &
           "thomp_mp_climo_file AND TRY AGAIN.",1)
         endif
       endif
     enddo
     close(14)
 
     num_tracers = num_tracers_input
     idum2(1) = num_tracers
     idum2(2) = num_tracers_input
 
     deallocate(var_type)
 
   endif  ! localpet = 0
 
   call esmf_vmbroadcast(vm, idum2, 2, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
   num_tracers       = idum2(1)
   num_tracers_input = idum2(2)
 
   call esmf_vmbroadcast(vm, tracers_input, len(tracers_input)*max_tracers, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
   call esmf_vmbroadcast(vm, chgres_var_names, len(chgres_var_names)*nvars, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
   call esmf_vmbroadcast(vm, field_var_names, len(field_var_names)*nvars, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
   call esmf_vmbroadcast(vm, missing_var_methods, len(missing_var_methods)*nvars, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
   call esmf_vmbroadcast(vm, missing_var_values, nvars, 0, rc=rc)
   if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
     call error_handler("IN VMBroadcast", rc)
 
 endif
end subroutine read_varmap
 
subroutine get_var_cond(var_name,this_miss_var_method,this_miss_var_value, &
                            this_field_var_name, loc)
  
  implicit none
  character(len=20), intent(in) :: var_name
  
  character(len=20), optional, intent(out) :: this_miss_var_method, &
                                              this_field_var_name
  real(esmf_kind_r4), optional, intent(out):: this_miss_var_value                                           
  
  integer, optional, intent(out)        :: loc
  
  integer                               :: i, tmp(size(missing_var_methods))
  
  i=0
  
  tmp(:)=0
  where(chgres_var_names==var_name) tmp=1
  
  i = maxloc(merge(1.,0.,chgres_var_names == var_name),dim=1) !findloc(chgres_var_names,var_name)
  print*, i
  if (maxval(tmp).eq.0) then
    print*, "WARNING: NO ENTRY FOR ", trim(var_name), " IN VARMAP TABLE. WILL SKIP " // &
            "VARIABLE IF NOT FOUND IN EXTERNAL MODEL FILE"
            
    if(present(this_miss_var_method)) this_miss_var_method = "skip"
    if(present(this_miss_var_value)) this_miss_var_value = -9999.9_esmf_kind_r4
    if(present(this_field_var_name)) this_field_var_name = "NULL"
    if(present(loc)) loc = 9999
  else
    if(present(this_miss_var_method)) this_miss_var_method = missing_var_methods(i)
    if(present(this_miss_var_value)) this_miss_var_value = missing_var_values(i)
    if(present(this_field_var_name)) this_field_var_name = field_var_names(i)
    if(present(loc)) loc = i
  endif
  
end subroutine get_var_cond
 
 subroutine calc_soil_params_driver(localpet)
 
 implicit none
 
 integer, intent(in)       :: localpet
 
 integer, parameter        :: num_statsgo = 16
 real, parameter           :: smlow_statsgo = 0.5
 real, parameter           :: smhigh_statsgo = 6.0
 
! zobler soil type used by spectral gfs prior to June 2017.
 integer, parameter        :: num_zobler = 9
 real, parameter           :: smlow_zobler = 0.5
 real, parameter           :: smhigh_zobler = 6.0
 
 integer                   :: num_soil_cats
 
 real                      :: bb_statsgo(num_statsgo)
 real                      :: maxsmc_statsgo(num_statsgo)
 real                      :: satdk_statsgo(num_statsgo)
 real                      :: satpsi_statsgo(num_statsgo)
 
 real                      :: bb_zobler(num_zobler)
 real                      :: maxsmc_zobler(num_zobler)
 real                      :: satdk_zobler(num_zobler)
 real                      :: satpsi_zobler(num_zobler)
 
 real, allocatable         :: bb(:)
 real                      :: smlow, smhigh
 real, allocatable         :: satdk(:)
 real, allocatable         :: satpsi(:)
 real, allocatable         :: satdw(:)
 
! using stasgo table
 data bb_statsgo /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
            6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
            5.25, -9.99, 4.05, 4.26/
 
 data maxsmc_statsgo /0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
              0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
              0.464, -9.99, 0.200, 0.421/
 
 data satdk_statsgo /1.7600e-4, 1.4078e-5, 5.2304e-6, 2.8089e-6, 2.8089e-6, &
             3.3770e-6, 4.4518e-6, 2.0348e-6, 2.4464e-6, 7.2199e-6, &
             1.3444e-6, 9.7384e-7, 3.3770e-6,     -9.99, 1.4078e-5, &
             1.4078e-5/
 
 data satpsi_statsgo /0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
              0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
              0.3548, -9.99,  0.0350, 0.0363/
 
 data bb_zobler /4.26,  8.72, 11.55,  4.74, 10.73,  8.17, &
                 6.77,  5.25,  4.26/
 
 data maxsmc_zobler /0.421, 0.464, 0.468, 0.434, 0.406, 0.465, &
                     0.404, 0.439, 0.421/
 
 data satdk_zobler /1.41e-5, 0.20e-5, 0.10e-5, 0.52e-5, 0.72e-5, &
                    0.25e-5, 0.45e-5, 0.34e-5, 1.41e-5/
 
 data satpsi_zobler /0.040, 0.620, 0.470, 0.140, 0.100, 0.260,  &
                     0.140, 0.360, 0.040/
 
!-------------------------------------------------------------------------
! Compute soil parameters for the input grid.
!-------------------------------------------------------------------------
 
 select case (trim(input_type))
   case ("gfs_sigio")
     print*,'- INPUT GRID USED ZOBLER SOIL TYPES.'
     num_soil_cats = num_zobler
   case default
     print*,'- INPUT GRID USED STATSGO SOIL TYPES.'
     num_soil_cats = num_statsgo
 end select
 
 allocate(maxsmc_input(num_soil_cats))
 allocate(wltsmc_input(num_soil_cats))
 allocate(drysmc_input(num_soil_cats))
 allocate(refsmc_input(num_soil_cats))
 allocate(bb(num_soil_cats))
 allocate(satdk(num_soil_cats))
 allocate(satpsi(num_soil_cats))
 allocate(satdw(num_soil_cats))
 
 select case (trim(input_type))
   case ("gfs_sigio")
     smlow  = smlow_zobler
     smhigh = smhigh_zobler
     maxsmc_input = maxsmc_zobler
     bb     = bb_zobler
     satdk  = satdk_zobler
     satpsi = satpsi_zobler
   case default
     smlow  = smlow_statsgo
     smhigh = smhigh_statsgo
     maxsmc_input = maxsmc_statsgo
     bb     = bb_statsgo
     satdk  = satdk_statsgo
     satpsi = satpsi_statsgo
 end select
 
 call calc_soil_params(num_soil_cats, smlow, smhigh, satdk, maxsmc_input, &
                       bb, satpsi, satdw, refsmc_input, drysmc_input, wltsmc_input)
 
 deallocate(bb, satdk, satpsi, satdw)
 
 if (localpet == 0) print*,'maxsmc input grid ',maxsmc_input
 if (localpet == 0) print*,'wltsmc input grid ',wltsmc_input
 
!-------------------------------------------------------------------------
! Compute soil parameters for the target grid.
!-------------------------------------------------------------------------
 
 print*,'- TARGET GRID USEING STATSGO SOIL TYPES.'
 
 num_soil_cats = num_statsgo
 
 allocate(maxsmc_target(num_soil_cats))
 allocate(wltsmc_target(num_soil_cats))
 allocate(drysmc_target(num_soil_cats))
 allocate(refsmc_target(num_soil_cats))
 allocate(bb_target(num_soil_cats))
 allocate(satpsi_target(num_soil_cats))
 allocate(satdk(num_soil_cats))
 allocate(satdw(num_soil_cats))
 
 smlow  = smlow_statsgo
 smhigh = smhigh_statsgo
 maxsmc_target = maxsmc_statsgo
 bb_target     = bb_statsgo
 satdk  = satdk_statsgo
 satpsi_target = satpsi_statsgo
 
 call calc_soil_params(num_soil_cats, smlow, smhigh, satdk, maxsmc_target, &
                       bb_target, satpsi_target, satdw, refsmc_target, drysmc_target, wltsmc_target)
 
 deallocate(satdk, satdw)
 
 if (localpet == 0) print*,'maxsmc target grid ',maxsmc_target
 if (localpet == 0) print*,'wltsmc input grid ',wltsmc_target
 
 end subroutine calc_soil_params_driver
 
 subroutine calc_soil_params(num_soil_cats, smlow, smhigh, satdk,  &
            maxsmc, bb, satpsi, satdw, refsmc, drysmc, wltsmc)
 
 implicit none
 
 integer, intent(in)            :: num_soil_cats
 
 real, intent(in)               :: smlow, smhigh
 real, intent(in)               :: bb(num_soil_cats)
 real, intent(in)               :: maxsmc(num_soil_cats)
 real, intent(in)               :: satdk(num_soil_cats)
 real, intent(in)               :: satpsi(num_soil_cats)
 
 real, intent(out)              :: satdw(num_soil_cats)
 real, intent(out)              :: refsmc(num_soil_cats)
 real, intent(out)              :: drysmc(num_soil_cats)
 real, intent(out)              :: wltsmc(num_soil_cats)
 
 integer                        :: i
 
 real                      :: refsmc1
 real                      :: wltsmc1
 
 satdw = 0.0
 refsmc = 0.0
 wltsmc = 0.0
 drysmc = 0.0
 
 do i = 1, num_soil_cats
 
   if (maxsmc(i) > 0.0) then
 
   satdw(i)  = bb(i)*satdk(i)*(satpsi(i)/maxsmc(i))
   refsmc1 = maxsmc(i)*(5.79e-9/satdk(i)) **(1.0/(2.0*bb(i)+3.0))
   refsmc(i) = refsmc1 + (maxsmc(i)-refsmc1) / smhigh
   wltsmc1 = maxsmc(i) * (200.0/satpsi(i))**(-1.0/bb(i))
   wltsmc(i) = wltsmc1 - smlow * wltsmc1
 
!----------------------------------------------------------------------
!  CURRENT VERSION DRYSMC VALUES THAT EQUATE TO WLTSMC.
!  FUTURE VERSION COULD LET DRYSMC BE INDEPENDENTLY SET VIA NAMELIST.
!----------------------------------------------------------------------
 
   drysmc(i) = wltsmc(i)
 
   end if
 
 END DO
 
 end subroutine calc_soil_params
 
 end module program_setup