phys/module_sf_ctsm.F

   1 module module_sf_ctsm
   2
   3   implicit none
   4   private
   5
   6 #ifdef WRF_USE_CTSM
   7
   8   public :: ctsm_init
   9   public :: ctsm_run
  10
  11   ! FIXME(wjs, 2020-01-01) Introduce a lilac_kinds module, and get this from there
  12   integer, parameter :: r8 = selected_real_kind(12)  ! 8 byte real
  13
  14 contains
  15
  16   !-----------------------------------------------------------------------
  17   subroutine get_num_points(ide, jde, its, ite, jts, jte, &
  18        num_points, ite_limited, jte_limited)
  19     ! Return the number of points owned by this task
  20     !
  21     ! If ite_limited and/or jte_limited are provided, then return those values, too.
  22     ! These are the task end indices on the mass point grid.
  23     integer, intent(in) :: ide  ! domain end index, i
  24     integer, intent(in) :: jde  ! domain end index, j
  25     integer, intent(in) :: its  ! task start index, i
  26     integer, intent(in) :: ite  ! task end index, i
  27     integer, intent(in) :: jts  ! task start index, j
  28     integer, intent(in) :: jte  ! task end index, j
  29     integer, intent(out) :: num_points  ! number of points owned by this task
  30     integer, optional, intent(out) :: ite_limited  ! task end index on the mass point grid, i
  31     integer, optional, intent(out) :: jte_limited  ! task end index on the mass point grid, j
  32
  33     integer :: my_ite_limited  ! task end index on the mass point grid, i
  34     integer :: my_jte_limited  ! task end index on the mass point grid, j
  35
  36     ! The very last index in both row & column space is just used on the momentum grid.
  37     ! Here we are just working with the mass point grid, so we need to ignore that last
  38     ! index.
  39     my_ite_limited = min(ite, ide-1)
  40     my_jte_limited = min(jte, jde-1)
  41
  42     num_points = ((my_ite_limited - its + 1) * (my_jte_limited - jts + 1))
  43
  44     if (present(ite_limited)) then
  45        ite_limited = my_ite_limited
  46     end if
  47     if (present(jte_limited)) then
  48        jte_limited = my_jte_limited
  49     end if
  50
  51   end subroutine get_num_points
  52
  53   subroutine create_gindex(ide, jde, its, ite, jts, jte, gindex)
  54     ! Create a gindex array on each task. This gives the list of global indices owned by
  55     ! each processor, on the mass point grid.
  56     integer, intent(in) :: ide  ! domain end index, i
  57     integer, intent(in) :: jde  ! domain end index, j
  58     integer, intent(in) :: its  ! task start index, i
  59     integer, intent(in) :: ite  ! task end index, i
  60     integer, intent(in) :: jts  ! task start index, j
  61     integer, intent(in) :: jte  ! task end index, j
  62     integer, allocatable, intent(out) :: gindex(:)
  63
  64     integer :: ite_limited  ! task end index on the mass point grid, i
  65     integer :: jte_limited  ! task end index on the mass point grid, j
  66     integer :: num_points
  67     integer :: i, j, n
  68
  69     call get_num_points( &
  70          ide=ide, jde=jde, &
  71          its=its, ite=ite, &
  72          jts=jts, jte=jte, &
  73          num_points = num_points, &
  74          ite_limited = ite_limited, &
  75          jte_limited = jte_limited)
  76
  77     allocate(gindex(num_points))
  78
  79     n = 0
  80     do j = jts, jte_limited
  81        do i = its, ite_limited
  82           n = n + 1
  83           ! In the following, note that we use ide-1 rather than ide for the same reason
  84           ! that we need ite_limited: ide gives the domain end index on the momentum grid,
  85           ! but here we're just dealing with the mass point grid, which has one less point
  86           ! in each direction.
  87           gindex(n) = (j-1)*(ide-1) + i
  88        end do
  89     end do
  90   end subroutine create_gindex
  91   !-----------------------------------------------------------------------
  92
  93   !-----------------------------------------------------------------------
  94   logical function is_land(xland, xice, xice_threshold)
  95     ! Returns true if the given point is land, false if ocean or sea ice
  96     real , intent(in) :: xland          ! land mask (1 for land, 2 for water)
  97     real , intent(in) :: xice           ! fraction of grid that is seaice
  98     real , intent(in) :: xice_threshold ! fraction of grid determining seaice
  99
 100     if (xland >= 1.5) then
 101        ! ocean
 102        is_land = .false.
 103     else if (xice >= xice_threshold) then
 104        ! sea ice
 105        is_land = .false.
 106     else
 107        is_land = .true.
 108     end if
 109   end function is_land
 110   !-----------------------------------------------------------------------
 111
 112   !-----------------------------------------------------------------------
 113   subroutine convert_2d_to_1d (ide, jde, ims, ime, jms, jme, its, ite, jts, jte, var_2d, var_1d)
 114     ! Convert a 2-d native WRF array to a 1-d array appropriate for LILAC
 115     !
 116     ! Allocates var_1d here
 117
 118     use module_wrf_error
 119
 120     ! input/output variables
 121     integer , intent(in)  :: ide ! domain end index, i
 122     integer , intent(in)  :: jde ! domain end index, j
 123     integer , intent(in)  :: ims ! memory start index (includes halo cells), i
 124     integer , intent(in)  :: ime ! memory end index (includes halo cells), i
 125     integer , intent(in)  :: jms ! memory start index (includes halo cells), j
 126     integer , intent(in)  :: jme ! memory end index (includes halo cells), j
 127     integer , intent(in)  :: its ! task start index, i
 128     integer , intent(in)  :: ite ! task end index, i
 129     integer , intent(in)  :: jts ! task start index, j
 130     integer , intent(in)  :: jte ! task end index, j
 131     real    , intent(in)  :: var_2d(ims: , jms: )
 132     real(r8), allocatable, intent(out) :: var_1d(:)
 133
 134     ! local variables
 135     integer :: ite_limited  ! task end index on the mass point grid, i
 136     integer :: jte_limited  ! task end index on the mass point grid, j
 137     integer :: num_points
 138     integer :: i, j, n
 139
 140     if (.not. all(ubound(var_2d) == [ime, jme])) then
 141        call wrf_error_fatal('convert_2d_to_1d: incorrect bounds for var_2d')
 142     end if
 143
 144     call get_num_points( &
 145          ide=ide, jde=jde, &
 146          its=its, ite=ite, &
 147          jts=jts, jte=jte, &
 148          num_points = num_points, &
 149          ite_limited = ite_limited, &
 150          jte_limited = jte_limited)
 151
 152     allocate (var_1d (num_points))
 153
 154     n = 0
 155
 156     do j = jts, jte_limited
 157        do i = its, ite_limited
 158           n = n + 1
 159           var_1d(n) = var_2d(i,j)
 160        end do
 161     end do
 162
 163   end subroutine convert_2d_to_1d
 164   !-----------------------------------------------------------------------
 165
 166   !-----------------------------------------------------------------------
 167   subroutine convert_1d_to_2d (ide, jde, ims, ime, jms, jme, its, ite, jts, jte, &
 168        xland, xice, xice_threshold, &
 169        var_1d, var_2d)
 170     ! Convert a 1-d array from lilac to a 2-d native WRF array
 171     !
 172     ! The output array is only set for land points (as determined by xland, xice, xice_threshold)
 173
 174     use module_wrf_error
 175
 176     ! input/output variables
 177     integer  , intent(in) :: ide                 ! domain end index, i
 178     integer  , intent(in) :: jde                 ! domain end index, j
 179     integer  , intent(in) :: ims                 ! memory start index (includes halo cells), i
 180     integer  , intent(in) :: ime                 ! memory end index (includes halo cells), i
 181     integer  , intent(in) :: jms                 ! memory start index (includes halo cells), j
 182     integer  , intent(in) :: jme                 ! memory end index (includes halo cells), j
 183     integer  , intent(in) :: its                 ! task start index, i
 184     integer  , intent(in) :: ite                 ! task end index, i
 185     integer  , intent(in) :: jts                 ! task start index, j
 186     integer  , intent(in) :: jte                 ! task end index, j
 187     real     , intent(in) :: xland(ims: , jms: ) ! land mask (1 for land, 2 for water)
 188     real     , intent(in) :: xice(ims: , jms: )  ! fraction of grid that is seaice
 189     real     , intent(in) :: xice_threshold      ! fraction of grid determining seaice
 190     real(r8) , intent(in) :: var_1d(:)
 191     real     , intent(inout) :: var_2d(ims: , jms: )
 192
 193     ! local variables
 194     integer :: ite_limited  ! task end index on the mass point grid, i
 195     integer :: jte_limited  ! task end index on the mass point grid, j
 196     integer :: num_points
 197     integer :: i, j, n
 198
 199     if (.not. all(ubound(var_2d) == [ime, jme])) then
 200        call wrf_error_fatal('convert_1d_to_2d: incorrect bounds for var_2d')
 201     end if
 202
 203     call get_num_points( &
 204          ide=ide, jde=jde, &
 205          its=its, ite=ite, &
 206          jts=jts, jte=jte, &
 207          num_points = num_points, &
 208          ite_limited = ite_limited, &
 209          jte_limited = jte_limited)
 210
 211     if (.not. (ubound(var_1d, 1) == num_points)) then
 212        call wrf_error_fatal('convert_1d_to_2d: incorrect size of var_1d')
 213     end if
 214
 215     n = 0
 216
 217     do j = jts, jte_limited
 218        do i = its, ite_limited
 219           n = n + 1
 220           if (is_land(xland(i,j), xice(i,j), xice_threshold)) then
 221              var_2d(i,j) = var_1d(n)
 222           end if
 223        end do
 224     end do
 225
 226   end subroutine convert_1d_to_2d
 227   !-----------------------------------------------------------------------
 228
 229   !-----------------------------------------------------------------------
 230   subroutine ctsm_init( &
 231        ids, ide, jds, jde, &
 232        ims, ime, jms, jme, &
 233        its, ite, jts, jte, &
 234        dt, xlat, xlong,    &
 235        atm_restart)
 236     ! Initialize CTSM via LILAC
 237
 238     use module_wrf_error
 239     use lilac_mod, only : lilac_init1, lilac_init2
 240     use ctsm_LilacCouplingFieldIndices
 241
 242     integer, intent(in) :: ids
 243     integer, intent(in) :: ide
 244     integer, intent(in) :: jds
 245     integer, intent(in) :: jde
 246     integer, intent(in) :: ims
 247     integer, intent(in) :: ime
 248     integer, intent(in) :: jms
 249     integer, intent(in) :: jme
 250     integer, intent(in) :: its
 251     integer, intent(in) :: ite
 252     integer, intent(in) :: jts
 253     integer, intent(in) :: jte
 254
 255     real, intent(in) :: dt       ! time step (seconds)
 256     real, intent(in) :: xlat(ims: , jms: )  ! latitudes (deg)
 257     real, intent(in) :: xlong(ims: , jms: ) ! longitudes (deg)
 258
 259     logical, intent(in) :: atm_restart
 260     integer :: comp_comm
 261     integer , allocatable :: atm_global_index(:)
 262     real                  :: xlong_0360(ims:ime, jms:jme)
 263     real(r8) , allocatable    :: atm_lons(:)
 264     real(r8) , allocatable    :: atm_lats(:)
 265
 266     character(len=128) :: atm_calendar
 267     character(len=128) :: atm_starttype
 268     integer            :: atm_timestep
 269     integer            :: atm_start_year ! (yyyy)
 270     integer            :: atm_start_mon  ! (mm)
 271     integer            :: atm_start_day
 272     integer            :: atm_start_hour
 273     integer            :: atm_start_minute
 274     integer            :: atm_start_second ! seconds after the minute
 275     integer            :: atm_start_secs_since_midnight ! total seconds since midnight
 276
 277     integer            :: atm_global_nx
 278     integer            :: atm_global_ny
 279
 280     character(len=512) :: message
 281
 282     ! TODO(wjs, 2019-12-31) Is this the correct way to get the communicator?
 283     ! (See https://github.com/ESCOMP/CTSM/issues/1078)
 284     call wrf_get_dm_communicator(comp_comm)
 285
 286     call create_gindex(ide, jde, its, ite, jts, jte, atm_global_index)
 287
 288     ! Convert longitude from -180..180 to 0..360
 289     xlong_0360 = xlong
 290     where (xlong_0360 < 0)
 291        xlong_0360 = 360. + xlong_0360
 292     end where
 293
 294     ! reshape lats and lons to 1d for lilac
 295     call convert_2d_to_1d ( &
 296          ide=ide, jde=jde, &
 297          ims=ims, ime=ime, &
 298          jms=jms, jme=jme, &
 299          its=its, ite=ite, &
 300          jts=jts, jte=jte, &
 301          var_2d=xlong_0360, var_1d=atm_lons)
 302     call convert_2d_to_1d ( &
 303          ide=ide, jde=jde, &
 304          ims=ims, ime=ime, &
 305          jms=jms, jme=jme, &
 306          its=its, ite=ite, &
 307          jts=jts, jte=jte, &
 308          var_2d=xlat, var_1d=atm_lats)
 309
 310     atm_global_nx = ide-ids
 311     atm_global_ny = jde-jds
 312
 313     ! calendar stuff
 314
 315     ! TODO(wjs, 2019-12-31) Is this the appropriate way to get the start and end times?
 316     !
 317     ! Some specific questions:
 318     ! - Should we use some passed-in argument rather than getting the namelist values
 319     !   directly here?
 320     ! - nl_get_start_year vs. nl_get_simulation_start_year (and similar for other units):
 321     !   should we be using nl_get_simulation_* here?
 322     ! - Is it correct to form the total seconds from hour, minute & second like this?
 323     ! - I'm using 1 for the first argument; I think this gives the domain id; is this an
 324     !   okay thing to do for now, or can we somehow determine the actual domain id?
 325     !
 326     ! (See https://github.com/ESCOMP/CTSM/issues/1078)
 327     call nl_get_start_year(1, atm_start_year)
 328     call nl_get_start_month(1, atm_start_mon)
 329     call nl_get_start_day(1, atm_start_day)
 330     call nl_get_start_hour(1, atm_start_hour)
 331     call nl_get_start_minute(1, atm_start_minute)
 332     call nl_get_start_second(1, atm_start_second)
 333     atm_start_secs_since_midnight = 60*(60*atm_start_hour + atm_start_minute) + atm_start_second
 334     write(message, '("CTSM start time: ", I4, 1X, I2, 1X, I2, 1X, I5)') &
 335          atm_start_year, atm_start_mon, atm_start_day, atm_start_secs_since_midnight
 336     call wrf_message(message)
 337
 338     atm_calendar      = 'GREGORIAN'
 339
 340     atm_starttype     = 'startup'
 341     if (atm_restart) then
 342         atm_starttype = 'continue'
 343     endif
 344
 345
 346     ! TODO(wjs, 2019-12-31) Is there a way to directly get dt as an integer, rather than
 347     ! relying on converting this real-valued dt to an integer?
 348     ! (See https://github.com/ESCOMP/CTSM/issues/1078)
 349     atm_timestep      = nint(dt)
 350     if (abs(atm_timestep - dt) > (1.e-5 * dt)) then
 351        call wrf_error_fatal('ctsm_init: expect dt representable as integer')
 352     end if
 353
 354     call lilac_init1()
 355     call lilac_init2(  &
 356          mpicom           = comp_comm,        &
 357          atm_global_index = atm_global_index, &
 358          atm_lons         = atm_lons,         &
 359          atm_lats         = atm_lats,         &
 360          atm_global_nx    = atm_global_nx,    &
 361          atm_global_ny    = atm_global_ny,    &
 362          atm_calendar     = atm_calendar,     &
 363          atm_timestep     = atm_timestep,     &
 364          atm_start_year   = atm_start_year,   &
 365          atm_start_mon    = atm_start_mon,    &
 366          atm_start_day    = atm_start_day,    &
 367          atm_start_secs   = atm_start_secs_since_midnight, &
 368          starttype_in     = atm_starttype, &
 369          fields_needed_from_data = [ &
 370          lilac_a2l_Faxa_bcphidry, lilac_a2l_Faxa_bcphodry, lilac_a2l_Faxa_bcphiwet, &
 371          lilac_a2l_Faxa_ocphidry, lilac_a2l_Faxa_ocphodry, lilac_a2l_Faxa_ocphiwet, &
 372          lilac_a2l_Faxa_dstwet1, lilac_a2l_Faxa_dstdry1, &
 373          lilac_a2l_Faxa_dstwet2, lilac_a2l_Faxa_dstdry2, &
 374          lilac_a2l_Faxa_dstwet3, lilac_a2l_Faxa_dstdry3, &
 375          lilac_a2l_Faxa_dstwet4, lilac_a2l_Faxa_dstdry4])
 376
 377   end subroutine ctsm_init
 378   !-----------------------------------------------------------------------
 379
 380   !-----------------------------------------------------------------------
 381   subroutine ctsm_run( &
 382        ! bounds
 383        ids, ide, jds, jde, &
 384        ims, ime, jms, jme, kms, &
 385        its, ite, jts, jte, &
 386
 387        ! restart flag
 388        restart_flag, &
 389
 390        ! general information
 391        dt, xland, xice, xice_threshold, &
 392
 393        ! atm -> lnd variables
 394        dz8w, ht, u_phy, v_phy, p8w, t_phy, th_phy, &
 395        qv_curr, rainbl, sr, &
 396        glw, swvisdir, swvisdif, swnirdir, swnirdif, &
 397
 398        ! lnd -> atm variables
 399        tsk, t2, qsfc, albedo, &
 400        ust, hfx, lh, qfx, emiss, z0, znt)
 401
 402     use lilac_mod, only : lilac_run
 403     use ctsm_LilacCouplingFieldIndices
 404
 405     integer, intent(in) :: ids
 406     integer, intent(in) :: ide
 407     integer, intent(in) :: jds
 408     integer, intent(in) :: jde
 409     integer, intent(in) :: ims
 410     integer, intent(in) :: ime
 411     integer, intent(in) :: jms
 412     integer, intent(in) :: jme
 413     integer, intent(in) :: kms
 414     integer, intent(in) :: its
 415     integer, intent(in) :: ite
 416     integer, intent(in) :: jts
 417     integer, intent(in) :: jte
 418
 419
 420     logical, intent(in) :: restart_flag      ! restart flag sent from WRF to CTSM
 421
 422     real , intent(in) :: dt                  ! timestep [s]
 423     real , intent(in) :: xland(ims: , jms: ) ! land mask (1 for land, 2 for water)
 424     real , intent(in) :: xice(ims: , jms: )  ! fraction of grid that is seaice
 425     real , intent(in) :: xice_threshold      ! fraction of grid determining seaice
 426
 427     real, intent(in) :: dz8w(ims: , kms: , jms: ) ! thickness of atmo layers [m]
 428     real, intent(in) :: ht(ims: , jms: )  ! terrain height [m]
 429     real, intent(in) :: u_phy(ims: , kms: , jms: ) ! 3D U wind component [m/s]
 430     real, intent(in) :: v_phy(ims: , kms: , jms: ) ! 3D V wind component [m/s]
 431     real, intent(in) :: p8w(ims: , kms: , jms: ) ! 3D pressure, valid at interface [Pa]
 432     real, intent(in) :: t_phy(ims: , kms: , jms: ) ! 3D atmospheric temperature valid at mid-levels [K]
 433     real, intent(in) :: th_phy(ims: , kms: , jms: ) ! 3D atmospheric temperature valid at mid-levels [K]
 434     real, intent(in) :: qv_curr(ims: , kms: , jms: ) ! 3D water vapor mixing ratio [kg/kg_dry]
 435     real, intent(in) :: rainbl(ims: , jms: ) ! total input precipitation [mm]
 436     real, intent(in) :: sr(ims: , jms: ) ! frozen precipitation ratio [-]
 437     real, intent(in) :: glw(ims: , jms: ) ! longwave down at surface [W m-2]
 438     real, intent(in) :: swvisdir(ims: , jms: ) ! vis direct beam solar rad onto surface [W m-2]
 439     real, intent(in) :: swvisdif(ims: , jms: ) ! vis diffuse solar rad onto surface [W m-2]
 440     real, intent(in) :: swnirdir(ims: , jms: ) ! nir direct beam solar rad onto surface [W m-2]
 441     real, intent(in) :: swnirdif(ims: , jms: ) ! nir diffuse solar rad onto surface [W m-2]
 442
 443     real, intent(inout) :: tsk(ims: , jms: ) ! surface temperature [K]
 444     real, intent(inout) :: t2(ims: , jms: ) ! diagnostic 2-m temperature [K]
 445     real, intent(inout) :: qsfc(ims: , jms: ) ! bulk surface specific humidity
 446     real, intent(inout) :: albedo(ims: , jms: ) ! total grid albedo [-]
 447     real, intent(inout) :: ust(ims: , jms: ) ! u* in similarity theory [m/s]
 448     real, intent(inout) :: hfx(ims: , jms: ) ! sensible heat flux [W m-2]
 449     real, intent(inout) :: lh(ims: , jms: ) ! latent heat flux [W m-2]
 450     real, intent(inout) :: qfx(ims: , jms: ) ! latent heat flux [kg s-1 m-2]
 451     real, intent(inout) :: emiss(ims: , jms: ) ! surface emissivity [between 0 and 1]
 452     real, intent(inout) :: z0(ims: , jms: ) ! background roughness length [m]
 453     real, intent(inout) :: znt(ims: , jms: ) ! thermal time-varying roughness length [m]
 454
 455     logical, save :: first_call = .true. ! true if and only if this is the first time this subroutine has been called in this run
 456
 457     integer :: i, j
 458
 459     real :: landmask(ims:ime, jms:jme) ! 1 over land, 0 over non-land (ocean, sea ice)
 460     real :: zlvl(ims:ime, jms:jme)   ! mid-point of bottom atm layer [m]
 461     real :: forc_u(ims:ime, jms:jme) ! u wind component, bottom atm layer [m/s]
 462     real :: forc_v(ims:ime, jms:jme) ! v wind component, bottom atm layer [m/s]
 463     real :: pbot(ims:ime, jms:jme)   ! surface pressure [Pa]
 464     real :: forc_th(ims:ime, jms:jme) ! potential temperature [K]
 465     real :: forc_t(ims:ime, jms:jme) ! temperature [K]
 466     real :: forc_q(ims:ime, jms:jme) ! specific humidity [kg/kg]
 467     real :: precip_rate(ims:ime, jms:jme) ! total precipitation rate [mm/s]
 468     real :: rain_convective_rate(ims:ime, jms:jme) ! rate of convective rain [mm/s]
 469     real :: rain_largescale_rate(ims:ime, jms:jme) ! rate of large-scale rain [mm/s]
 470     real :: snow_convective_rate(ims:ime, jms:jme) ! rate of convective snow [mm/s]
 471     real :: snow_largescale_rate(ims:ime, jms:jme) ! rate of large-scale snow [mm/s]
 472
 473     real :: qref(ims:ime, jms:jme) ! 2m surface specific humidity calculated by CTSM [kg/kg]
 474     real :: albedo_visdir(ims:ime, jms:jme) ! albedo calculated by CTSM, visible direct [-]
 475     real :: albedo_visdif(ims:ime, jms:jme) ! albedo calculated by CTSM, visible diffuse [-]
 476     real :: albedo_nirdir(ims:ime, jms:jme) ! albedo calculated by CTSM, near IR direct [-]
 477     real :: albedo_nirdif(ims:ime, jms:jme) ! albedo calculated by CTSM, near IR diffuse [-]
 478
 479     logical  :: no_negative   = .true.      ! flag for converting negative values to zero
 480
 481     ! ------------------------------------------------------------------------
 482     ! Calculate landmask and export it to ctsm via lilac.
 483     !
 484     ! For efficiency, only do this the first time this subroutine is called: we assume
 485     ! landmask stays constant throughout the run (there can be transitions between open
 486     ! ocean and sea ice, but the land itself stays fixed).
 487     ! ------------------------------------------------------------------------
 488
 489     if (first_call) then
 490        do j = jts, jte
 491           do i = its, ite
 492              if (is_land(xland(i,j), xice(i,j), xice_threshold)) then
 493                 landmask(i,j) = 1.
 494              else
 495                 landmask(i,j) = 0.
 496              end if
 497           end do
 498        end do
 499        call export_to_lilac(lilac_a2l_Sa_landfrac, landmask)
 500     end if
 501
 502     ! ------------------------------------------------------------------------
 503     ! Calculate derived variables, and 2-d versions of 3-d fields
 504     ! ------------------------------------------------------------------------
 505
 506     ! dz8w = thickness of full levels; we want the mid-point of the bottom-most level
 507     zlvl(its:ite, jts:jte) = 0.5 * dz8w(its:ite, 1, jts:jte)
 508
 509     forc_u(its:ite, jts:jte) = u_phy(its:ite, 1, jts:jte)
 510     forc_v(its:ite, jts:jte) = v_phy(its:ite, 1, jts:jte)
 511
 512     pbot(its:ite, jts:jte) = p8w(its:ite, 1, jts:jte)
 513
 514     forc_th(its:ite, jts:jte) = th_phy(its:ite, 1, jts:jte)
 515
 516     forc_t(its:ite, jts:jte) = t_phy(its:ite, 1, jts:jte)
 517
 518     ! Convert from mixing ratio to specific humidity
 519     forc_q(its:ite, jts:jte) = qv_curr(its:ite, 1, jts:jte)/(1.0 + qv_curr(its:ite, 1, jts:jte))
 520
 521     ! Separate total precip into rain and snow. Arbitrarily assign all precipitation to
 522     ! convective (CTSM requires separate convective vs. large-scale precipitation, but
 523     ! then just adds the two together).
 524     precip_rate(its:ite, jts:jte) = rainbl(its:ite, jts:jte)/dt
 525     snow_convective_rate(its:ite, jts:jte) = precip_rate(its:ite, jts:jte) * sr(its:ite, jts:jte)
 526     snow_largescale_rate(its:ite, jts:jte) = 0.
 527     rain_convective_rate(its:ite, jts:jte) = precip_rate(its:ite, jts:jte) * (1. - sr(its:ite, jts:jte))
 528     rain_largescale_rate(its:ite, jts:jte) = 0.
 529
 530     ! ------------------------------------------------------------------------
 531     ! Export data to ctsm via lilac
 532     ! ------------------------------------------------------------------------
 533
 534     call export_to_lilac(lilac_a2l_Sa_z, zlvl)
 535     call export_to_lilac(lilac_a2l_Sa_topo, ht)
 536     call export_to_lilac(lilac_a2l_Sa_u, forc_u)
 537     call export_to_lilac(lilac_a2l_Sa_v, forc_v)
 538     call export_to_lilac(lilac_a2l_Sa_ptem, forc_th)
 539     call export_to_lilac(lilac_a2l_Sa_pbot, pbot)
 540     call export_to_lilac(lilac_a2l_Sa_tbot, forc_t)
 541     call export_to_lilac(lilac_a2l_Sa_shum, forc_q)
 542     call export_to_lilac(lilac_a2l_Faxa_lwdn, glw)
 543     call export_to_lilac(lilac_a2l_Faxa_rainc, rain_convective_rate)
 544     call export_to_lilac(lilac_a2l_Faxa_rainl, rain_largescale_rate)
 545     call export_to_lilac(lilac_a2l_Faxa_snowc, snow_convective_rate)
 546     call export_to_lilac(lilac_a2l_Faxa_snowl, snow_largescale_rate)
 547     call export_to_lilac(lilac_a2l_Faxa_swndr, swnirdir, no_negative)
 548     call export_to_lilac(lilac_a2l_Faxa_swvdr, swvisdir, no_negative)
 549     call export_to_lilac(lilac_a2l_Faxa_swndf, swnirdif, no_negative)
 550     call export_to_lilac(lilac_a2l_Faxa_swvdf, swvisdif, no_negative)
 551
 552     ! ------------------------------------------------------------------------
 553     ! Run ctsm via lilac
 554     ! ------------------------------------------------------------------------
 555
 556     ! FIXME(wjs, 2020-01-01) Use correct values for restart and stop alarms
 557     call lilac_run( &
 558          write_restarts_now = restart_flag, &
 559          stop_now = .false.)
 560
 561     ! ------------------------------------------------------------------------
 562     ! Import data from ctsm via lilac
 563     ! ------------------------------------------------------------------------
 564
 565     call import_from_lilac(lilac_l2a_Sl_t, tsk)
 566     call import_from_lilac(lilac_l2a_Sl_tref, t2)
 567     call import_from_lilac(lilac_l2a_Sl_qref, qref)
 568
 569     call import_from_lilac(lilac_l2a_Sl_avsdr, albedo_visdir)
 570     call import_from_lilac(lilac_l2a_Sl_anidr, albedo_nirdir)
 571     call import_from_lilac(lilac_l2a_Sl_avsdf, albedo_visdif)
 572     call import_from_lilac(lilac_l2a_Sl_anidf, albedo_nirdif)
 573     call calculate_total_albedo( &
 574          ims=ims, ime=ime, jms=jms, jme=jme, &
 575          its=its, ite=ite, jts=jts, jte=jte, &
 576          xland          = xland, &
 577          xice           = xice, &
 578          xice_threshold = xice_threshold, &
 579          albedo_visdir  = albedo_visdir, &
 580          albedo_nirdir  = albedo_nirdir, &
 581          albedo_visdif  = albedo_visdif, &
 582          albedo_nirdif  = albedo_nirdif, &
 583          swvisdir       = swvisdir, &
 584          swnirdir       = swnirdir, &
 585          swvisdif       = swvisdif, &
 586          swnirdif       = swnirdif, &
 587          albedo         = albedo)
 588
 589     call import_from_lilac(lilac_l2a_Sl_fv, ust)
 590
 591     call import_from_lilac(lilac_l2a_Fall_sen, hfx)
 592     call import_from_lilac(lilac_l2a_Fall_lat, lh)
 593     call import_from_lilac(lilac_l2a_Fall_evap, qfx)
 594
 595     ! At least for now, use CTSM's Sl_z0m (momentum roughness length) for both background
 596     ! roughness length and thermal time-varying roughness length, even though it's
 597     ! possible that those two should differ in some way.
 598     call import_from_lilac(lilac_l2a_Sl_z0m, z0)
 599
 600     ! ------------------------------------------------------------------------
 601     ! Calculate derived variables
 602     ! ------------------------------------------------------------------------
 603
 604     do j = jts, jte
 605        do i = its, ite
 606           if (is_land(xland(i,j), xice(i,j), xice_threshold)) then
 607
 608              ! Convert from specific humidity to mixing ratio. Note that qref is specific
 609              ! humidity at 2m, whereas qsfc is supposed to be specified right at the
 610              ! surface. So there isn't a perfect correspondence between the two, but
 611              ! given that qsfc is just being used as a diagnostic quantity when using
 612              ! CTSM (for now), we won't worry about this.
 613              qsfc(i,j) = qref(i,j) / (1. - qref(i,j))
 614
 615              ! CTSM assumes an emissivity of 1
 616              emiss(i,j) = 1.
 617
 618              ! At least for now, use CTSM's Sl_z0m (momentum roughness length) for both
 619              ! background roughness length and thermal time-varying roughness length, even
 620              ! though it's possible that those two should differ in some way.
 621              znt(i,j) = z0(i,j)
 622
 623              ! Flip signs of heat and moisture fluxes: CTSM sends as positive downwards;
 624              ! WRF expects them to be positive upwards
 625              hfx(i,j) = -hfx(i,j)
 626              lh(i,j) = -lh(i,j)
 627              qfx(i,j) = -qfx(i,j)
 628           end if
 629        end do
 630     end do
 631
 632     first_call = .false.
 633
 634   contains
 635
 636     subroutine export_to_lilac(field_index, var_2d , no_negative)
 637       ! Reshape var_2d to 1d for LILAC, then set the appropriate atm2lnd variable in LILAC
 638
 639       use ctsm_LilacCouplingFields, only : lilac_atm2lnd
 640
 641       integer, intent(in) :: field_index  ! one of the field indices defined in ctsm_LilacCouplingFieldIndices
 642       real, intent(in) :: var_2d(ims: , jms: )
 643
 644       real(r8), allocatable :: var_1d(:)
 645
 646       logical , optional, intent(in) :: no_negative  ! convert negative values to zero
 647
 648       call convert_2d_to_1d( &
 649            ide=ide, jde=jde, &
 650            ims=ims, ime=ime, &
 651            jms=jms, jme=jme, &
 652            its=its, ite=ite, &
 653            jts=jts, jte=jte, &
 654            var_2d=var_2d, var_1d=var_1d)
 655
 656       ! zero-ing out negative value when we no_negative flag.
 657       if (present (no_negative)) then
 658         if (no_negative) then
 659             var_1d = 0.5* (var_1d +abs(var_1d))
 660         end if
 661       end if
 662
 663       call lilac_atm2lnd( &
 664            field_index = field_index, &
 665            data = var_1d)
 666
 667     end subroutine export_to_lilac
 668
 669     subroutine import_from_lilac(field_index, var_2d)
 670       ! Get the appropriate lnd2atm variable from LILAC, then reshape 1d LILAC variable to
 671       ! fill var_2d
 672
 673       use ctsm_LilacCouplingFields, only : lilac_lnd2atm
 674
 675       integer, intent(in) :: field_index  ! one of the field indices defined in ctsm_LilacCouplingFieldIndices
 676       real, intent(inout) :: var_2d(ims: , jms: )
 677
 678       real(r8), allocatable :: var_1d(:)
 679       integer :: num_points
 680
 681       call get_num_points( &
 682            ide=ide, jde=jde, &
 683            its=its, ite=ite, &
 684            jts=jts, jte=jte, &
 685            num_points = num_points)
 686
 687       allocate(var_1d(num_points))
 688
 689       call lilac_lnd2atm( &
 690            field_index = field_index, &
 691            data = var_1d)
 692
 693       call convert_1d_to_2d( &
 694            ide=ide, jde=jde, &
 695            ims=ims, ime=ime, &
 696            jms=jms, jme=jme, &
 697            its=its, ite=ite, &
 698            jts=jts, jte=jte, &
 699            xland=xland, xice=xice, xice_threshold=xice_threshold, &
 700            var_1d=var_1d, var_2d=var_2d)
 701
 702     end subroutine import_from_lilac
 703
 704   end subroutine ctsm_run
 705   !-----------------------------------------------------------------------
 706
 707   !-----------------------------------------------------------------------
 708   subroutine calculate_total_albedo( &
 709        ims, ime, jms, jme, its, ite, jts, jte, &
 710        xland, xice, xice_threshold, &
 711        albedo_visdir, albedo_nirdir, albedo_visdif, albedo_nirdif, &
 712        swvisdir, swnirdir, swvisdif, swnirdif, &
 713        albedo)
 714     ! Calculate total albedo from its 4 components, based on ratio of incoming SW in the
 715     ! 4 components.
 716
 717     integer, intent(in) :: ims
 718     integer, intent(in) :: ime
 719     integer, intent(in) :: jms
 720     integer, intent(in) :: jme
 721     integer, intent(in) :: its
 722     integer, intent(in) :: ite
 723     integer, intent(in) :: jts
 724     integer, intent(in) :: jte
 725
 726     real, intent(in) :: xland(ims: , jms: ) ! land mask (1 for land, 2 for water)
 727     real, intent(in) :: xice(ims: , jms: )  ! fraction of grid that is seaice
 728     real, intent(in) :: xice_threshold      ! fraction of grid determining seaice
 729
 730     real, intent(in) :: albedo_visdir(ims: , jms: ) ! albedo, visible direct [-]
 731     real, intent(in) :: albedo_nirdir(ims: , jms: ) ! albedo, near IR direct [-]
 732     real, intent(in) :: albedo_visdif(ims: , jms: ) ! albedo, visible diffuse [-]
 733     real, intent(in) :: albedo_nirdif(ims: , jms: ) ! albedo, near IR diffuse [-]
 734     real, intent(in) :: swvisdir(ims: , jms: ) ! solar rad onto surface, visible direct [W m-2]
 735     real, intent(in) :: swnirdir(ims: , jms: ) ! solar rad onto surface, near IR direct [W m-2]
 736     real, intent(in) :: swvisdif(ims: , jms: ) ! solar rad onto surface, visible diffuse [W m-2]
 737     real, intent(in) :: swnirdif(ims: , jms: ) ! solar rad onto surface, near IR diffuse [W m-2]
 738
 739     real, intent(inout) :: albedo(ims: , jms: ) ! total surface albedo [-]
 740
 741     real, parameter :: albedo_default = 0.3
 742
 743     integer :: i, j
 744     real :: sw_tot  ! total solar rad onto surface [W m-2]
 745
 746     do j = jms, jme
 747        do i = ims, ime
 748           if (is_land(xland(i,j), xice(i,j), xice_threshold)) then
 749              sw_tot = swvisdir(i,j) + swnirdir(i,j) + swvisdif(i,j) + swnirdif(i,j)
 750              if (sw_tot > 0.) then
 751                 albedo(i,j) = &
 752                      albedo_visdir(i,j) * (swvisdir(i,j) / sw_tot) + &
 753                      albedo_nirdir(i,j) * (swnirdir(i,j) / sw_tot) + &
 754                      albedo_visdif(i,j) * (swvisdif(i,j) / sw_tot) + &
 755                      albedo_nirdif(i,j) * (swnirdif(i,j) / sw_tot)
 756              else
 757                 ! Night; albedo shouldn't matter; use coefficients from module_sf_clm
 758                 albedo(i,j) = &
 759                      albedo_visdir(i,j) * 0.35 + &
 760                      albedo_nirdir(i,j) * 0.35 + &
 761                      albedo_visdif(i,j) * 0.15 + &
 762                      albedo_nirdif(i,j) * 0.15
 763              end if
 764
 765              if (abs(albedo(i,j) - 1.) < 1.e-5) then
 766                 ! CTSM gives albedo values of 1 at night. To avoid problems (in case CTSM
 767                 ! thinks it's night but the albedo is still needed in WRF), use some more
 768                 ! reasonable value in this case.
 769                 albedo(i,j) = albedo_default
 770              end if
 771           end if
 772        end do
 773     end do
 774   end subroutine calculate_total_albedo
 775
 776 #endif
 777   ! endif WRF_USE_CTSM
 778
 779 end module module_sf_ctsm