Modifications for DAcycling=true

src/core_atmosphere/mpas_atm_core.F

@@ -58,8 +58,9 @@ function atm_core_init(domain, startTimeStamp) result(ierr)
       real (kind=RKIND), pointer :: dt
       type (block_type), pointer :: block
 
-      logical, pointer :: config_do_restart
 
+      logical, pointer :: config_do_restart
+      logical, pointer :: config_do_DAcycling
       type (mpas_pool_type), pointer :: state
       type (mpas_pool_type), pointer :: mesh
       type (mpas_pool_type), pointer :: diag
@@ -116,6 +117,7 @@ function atm_core_init(domain, startTimeStamp) result(ierr)
       end if
 
       call mpas_pool_get_config(domain % blocklist % configs, 'config_do_restart', config_do_restart)
+      call mpas_pool_get_config(domain % blocklist % configs, 'config_do_DAcycling', config_do_DAcycling)
       call mpas_pool_get_config(domain % blocklist % configs, 'config_dt', dt)
 
       !
@@ -137,6 +139,8 @@ function atm_core_init(domain, startTimeStamp) result(ierr)
       call mpas_timer_start('read_ICs')
       if (config_do_restart) then
          init_stream_name = 'restart'
+      elseif (config_do_DAcycling) then
+         init_stream_name = 'da_state'
       else
          init_stream_name = 'input'
       end if
@@ -158,6 +162,7 @@ function atm_core_init(domain, startTimeStamp) result(ierr)
 
       call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='input', direction=MPAS_STREAM_INPUT, ierr=ierr)
       call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_INPUT, ierr=ierr)
+      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='da_state', direction=MPAS_STREAM_INPUT, ierr=ierr)
       call mpas_log_write(' ----- done reading initial state -----')
 
 
@@ -642,6 +647,8 @@ function atm_core_run(domain) result(ierr)
 
       ! Avoid writing a restart file at the initial time
       call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='restart', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
+      ! Avoid writing a mpasout file at the initial time
+      call MPAS_stream_mgr_reset_alarms(domain % streamManager, streamID='da_state', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
 
       ! Also, for restart runs, avoid writing the initial history or diagnostics fields to avoid overwriting those from the preceding run
       if (config_do_restart) then

src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F

@@ -903,6 +903,7 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
 
 !local pointers:
  logical,pointer:: config_sfc_albedo
+ logical,pointer:: restart, cycling
  character(len=StrKIND),pointer:: lsm_scheme
  character(len=StrKIND),pointer:: mminlu
  integer,pointer:: isice
@@ -918,6 +919,8 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_config(configs,'config_sfc_albedo' ,config_sfc_albedo )
  call mpas_pool_get_config(configs,'config_lsm_scheme',lsm_scheme)
  call mpas_pool_get_config(configs,'config_lakemodel',lakemodel)
+ call mpas_pool_get_config(configs,'config_do_restart',restart)
+ call mpas_pool_get_config(configs,'config_do_DAcycling',cycling)
  call mpas_pool_get_array(sfc_input,'mminlu',mminlu)
  call mpas_pool_get_array(sfc_input,'isice' ,isice )
  call mpas_pool_get_array(sfc_input,'iswater',iswater)
@@ -989,6 +992,7 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
 
        call mpas_log_write('--- call subroutine ruc_land:')
        call ruc_land( spp_lsm      = spp_lsm_loc   , lakemodel    = lakemodel   , lakemask  = lakemask_p,  &
+                    restart        = restart       , cycling      = cycling     ,                          &
                     rhosnf         = rhosnf_p      , precipfr     = precipfr_p  , mosaic_lu = mosaic_lu,   &
                     qsg            = qsg_p         , qvg          = qvg_p       , dew       = dew_p,       &
                     soilt1         = soilt1_p      , tsnav        = tsnav_p     , acrunoff  = acrunoff_p,  &

src/core_atmosphere/physics/mpas_atmphys_driver_seaice.F

@@ -695,12 +695,16 @@ subroutine driver_seaice(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
  integer,pointer:: isice,iswater
  integer,pointer:: spp_lsm_loc 
  integer,pointer:: globalCells(:)
+ logical,pointer:: restart, cycling
 
 !-----------------------------------------------------------------------------------------------------------------
 !call mpas_log_write(' ')
 !call mpas_log_write('--- enter subroutine driver_seaice: xice_threshold = $r',realArgs=(/xice_threshold/))
 
  call mpas_pool_get_config(configs,'config_lsm_scheme',lsm_scheme)
+ call mpas_pool_get_config(configs,'config_do_restart',restart)
+ call mpas_pool_get_config(configs,'config_do_DAcycling',cycling)
+
 
 !copy MPAS arrays to local arrays:
  call seaice_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
@@ -778,6 +782,7 @@ subroutine driver_seaice(itimestep,configs,mesh,diag_physics,sfc_input,its,ite)
 
        call mpas_log_write('--- call subroutine ruc_ice:')
        call ruc_ice( spp_lsm       = spp_lsm_loc   , iswater      = iswater     , chklowq   = chklowq_p,   &
+                    restart        = restart       , cycling      = cycling     ,                          &
                     rhosnf         = rhosnf_p      , precipfr     = precipfr_p  , chs       = chs_p,       &
                     qsg            = qsg_p         , qvg          = qvg_p       , dew       = dew_p,       &
                     soilt1         = soilt1_p      , tsnav        = tsnav_p     , acrunoff  = acrunoff_p,  &

src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F

@@ -1162,6 +1162,7 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
 
 !local pointers:
  logical,pointer:: config_do_restart,config_frac_seaice
+ logical,pointer:: config_do_DAcycling     
  character(len=StrKIND),pointer:: sfclayer_scheme
  real(kind=RKIND),dimension(:),pointer:: areaCell
 
@@ -1181,10 +1182,11 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
  errmsg = ' '
  errflg = 0
 
- call mpas_pool_get_config(configs,'config_do_restart'     ,config_do_restart )
- call mpas_pool_get_config(configs,'config_frac_seaice'    ,config_frac_seaice)
- call mpas_pool_get_config(configs,'config_sfclayer_scheme',sfclayer_scheme   )
-
+ call mpas_pool_get_config(configs,'config_do_restart'     ,config_do_restart  )
+ call mpas_pool_get_config(configs,'config_frac_seaice'    ,config_frac_seaice )
+ call mpas_pool_get_config(configs,'config_sfclayer_scheme',sfclayer_scheme    )
+ call mpas_pool_get_config(configs,'config_do_DAcycling'   ,config_do_DAcycling) 
+
  call mpas_pool_get_array(mesh,'areaCell',areaCell)
 
 !copy all MPAS arrays to rectanguler grid:
@@ -1193,7 +1195,7 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
  dx = sqrt(maxval(areaCell))
 
  initflag = 1
- if(config_do_restart .or. itimestep > 1) initflag = 0
+ if(config_do_restart .or. config_do_DAcycling .or. itimestep > 1) initflag = 0
 
  sfclayer_select: select case (trim(sfclayer_scheme))
 
@@ -1413,7 +1415,10 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
                    ustm     = ustm_p     , ck        = ck_p     , cka      = cka_p       , &
                    cd       = cd_p       , cda       = cda_p    , ch       = ch_p        , &
                    qcg      = qcg_p      , spp_pbl   = spp_pbl  , isftcflx = isftcflx    , &
-                   iz0tlnd  = iz0tlnd    , itimestep = initflag                          , &
+                   iz0tlnd  = iz0tlnd                                                    , &
+                   itimestep= itimestep  , initflag  = initflag                          , &
+                   restart  = config_do_restart                                          , &
+                   cycling  = config_do_DAcycling                                        , &
                    errmsg   = errmsg     , errflg    = errflg                            , &
                    ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , &
                    ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , &
@@ -1445,7 +1450,10 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
                       ustm     = ustm_sea   , ck        = ck_sea   , cka      = cka_sea     , &
                       cd       = cd_sea     , cda       = cda_sea  , ch       = ch_sea      , &
                       qcg      = qcg_p      , spp_pbl   = spp_pbl  , isftcflx = isftcflx    , &
-                      iz0tlnd  = iz0tlnd    , itimestep = initflag                          , &
+                      iz0tlnd  = iz0tlnd                                                    , &
+                      itimestep= itimestep  , initflag  = initflag                          , &
+                      restart  = config_do_restart                                          , &
+                      cycling  = config_do_DAcycling                                        , &
                       errmsg   = errmsg     , errflg    = errflg                            , &
                       ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , &
                       ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , &

src/core_atmosphere/physics/mpas_atmphys_init.F

@@ -107,6 +107,7 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
 
 !local pointers:
  logical,pointer:: config_do_restart,    &
+                   config_do_DAcycling,  &
                    config_o3climatology, &
                    config_oml1d
 
@@ -162,6 +163,7 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
 !call mpas_log_write('--- enter subroutine physics_init:')
 
  call mpas_pool_get_config(configs,'config_do_restart'       ,config_do_restart       )
+ call mpas_pool_get_config(configs,'config_do_DAcycling'     ,config_do_DAcycling     )
  call mpas_pool_get_config(configs,'config_o3climatology'    ,config_o3climatology    )
  call mpas_pool_get_config(configs,'config_convection_scheme',config_convection_scheme)
  call mpas_pool_get_config(configs,'config_lsm_scheme'       ,config_lsm_scheme       )
@@ -308,7 +310,7 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
  endif
 
 !initialization of xicem:
- if(.not.config_do_restart) then
+ if(.not.config_do_restart .and. .not. config_do_DAcycling) then
 !   call mpas_log_write('--- initialization of xicem:')
     do iCell = 1, nCellsSolve
        xicem(iCell) = xice(iCell)
@@ -318,7 +320,7 @@ subroutine physics_init(dminfo,stream_manager,clock,configs,mesh,diag,tend,state
 !initialization of the local sea-surface temperature when a diurnal cycle of the
 !sea-surface temperature is applied. This avoids having the array sstsk equal to
 !zero over land:
- if(.not. config_do_restart) then
+ if(.not. config_do_restart .and. .not. config_do_DAcycling) then
 !   call mpas_log_write('--- initialization of sstsk:')
     do iCell = 1, nCellsSolve
        sstsk(iCell) = sst(iCell)

src/core_atmosphere/physics/mpas_atmphys_landuse.F

@@ -96,6 +96,7 @@ subroutine landuse_init_forMPAS(dminfo,julday,mesh,configs,diag_physics,sfc_inpu
 
 !local pointers:
  logical,pointer:: config_do_restart,  &
+                   config_do_DAcycling,&
                    config_frac_seaice, &
                    config_sfc_albedo
 
@@ -135,6 +136,7 @@ subroutine landuse_init_forMPAS(dminfo,julday,mesh,configs,diag_physics,sfc_inpu
 !call mpas_log_write('--- enter subroutine landuse_init_forMPAS:')
 
  call mpas_pool_get_config(configs,'config_do_restart' ,config_do_restart )
+ call mpas_pool_get_config(configs,'config_do_DAcycling' ,config_do_DAcycling )
  call mpas_pool_get_config(configs,'config_frac_seaice',config_frac_seaice)
  call mpas_pool_get_config(configs,'config_sfc_albedo' ,config_sfc_albedo )
  call mpas_pool_get_config(configs,'config_lsm_scheme' ,config_lsm_scheme )
@@ -164,7 +166,8 @@ subroutine landuse_init_forMPAS(dminfo,julday,mesh,configs,diag_physics,sfc_inpu
  call mpas_pool_get_array(diag_physics,'xicem'            ,xicem            )
  call mpas_pool_get_array(diag_physics,'z0'               ,z0               )
  call mpas_pool_get_array(diag_physics,'znt'              ,znt              )
- if (config_lsm_scheme .eq. 'ruc') then
+
+ if (config_lsm_scheme .eq. 'sf_ruc') then
     call mpas_pooL_get_array(diag_physics,'spp_lsm_loc'   , spp_lsm_loc)
     spp_lsm_loc = 0   !hardcode as 0 for now
     call mpas_pool_get_array(diag_physics,'mosaic_lu'     , mosaic_lu  )
@@ -179,6 +182,7 @@ subroutine landuse_init_forMPAS(dminfo,julday,mesh,configs,diag_physics,sfc_inpu
 !call mpas_log_write('--- config_frac_seaice = $1',logicArgs=(/config_frac_seaice/))
 !call mpas_log_write('--- xice_threshold     = $r',realArgs=(/xice_threshold/))
 
+ if((config_lsm_scheme .ne. 'sf_ruc')) then
 !reads in the landuse properties from landuse.tbl:
  if(dminfo % my_proc_id == IO_NODE) then
     !get a unit to open init file:
@@ -264,11 +268,19 @@ subroutine landuse_init_forMPAS(dminfo,julday,mesh,configs,diag_physics,sfc_inpu
 !call mpas_log_write('--- isice   =$i',intArgs=(/isice/))
 !call mpas_log_write('--- iswater =$i',intArgs=(/iswater/))
 !call mpas_log_write('--- isurban =$i',intArgs=(/isurban/))
- if(config_do_restart) then
+ endif ! not sf_ruc
+
+ if(config_do_restart ) then
     call mpas_log_write('--- config_do_restart =$l', logicArgs=(/config_do_restart/))
     call mpas_log_write('--- skip the end of landuse_init_forMPAS')
     return
  endif
+  if(config_do_DAcycling) then
+    call mpas_log_write('--- config_do_DAcycling =$l', logicArgs=(/config_do_restart/))
+    call mpas_log_write('--- skip the end of landuse_init_forMPAS')
+    return
+ endif
+
 
 !defines the surface properties over the entire domain:
  do iCell = 1, nCells

src/core_atmosphere/physics/mpas_atmphys_lsm_rucinit.F

@@ -88,7 +88,7 @@ subroutine ruclsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  type(mpas_pool_type),intent(inout):: sfc_input
 
 !local pointers::
- logical,pointer:: input_sfc_albedo,restart
+ logical,pointer:: input_sfc_albedo,restart,cycling
 
  character(len=StrKIND),pointer:: mminlu,mminsl
 
@@ -122,6 +122,8 @@ subroutine ruclsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  call mpas_pool_get_config(configs,'input_soil_data'      ,mminsl          )
  call mpas_pool_get_config(configs,'config_sfc_snowalbedo',input_sfc_albedo)
  call mpas_pool_get_config(configs,'config_do_restart'    ,restart         )
+ call mpas_pool_get_config(configs,'config_do_DAcycling'  ,cycling         )
+
 
  call mpas_pool_get_dimension(mesh,'nCells'     ,nCells     )
  call mpas_pool_get_dimension(mesh,'nSoilLevels',nSoilLevels)
@@ -171,8 +173,7 @@ subroutine ruclsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  call physics_message('--- end initialize RUC LSM tables' )
  call mpas_log_write('')
 
-
-if(.not.restart) then
+if(.not.restart .and. .not. cycling) then
 
    errflag = 0
    do iCell = 1, nCells

src/core_atmosphere/physics/physics_wrf/module_mynnsfclay.F90

@@ -139,7 +139,8 @@ subroutine mynnsfclay(                                 &
                      svp3,svpt0,ep1,ep2,karman,ch,qcg,   &
                      itimestep,wstar,qstar,ustm,ck,cka,  &
                      cd,cda,spp_pbl,rstoch1d,isftcflx,   &
-                     iz0tlnd,its,ite,errmsg,errflg       &
+                     iz0tlnd,its,ite,restart_or_cycle,   &
+                     errmsg,errflg                       &
                         )
  implicit none
 !=================================================================================================================
@@ -162,7 +163,7 @@ subroutine mynnsfclay(                                 &
 !-----------------------------
 ! namelist options
 !-----------------------------
- logical,intent(in):: spp_pbl
+ logical,intent(in):: spp_pbl,restart_or_cycle
 
  integer,intent(in):: isfflx
  integer,intent(in),optional:: isftcflx,iz0tlnd
@@ -385,7 +386,7 @@ subroutine mynnsfclay(                                 &
     ! according to Akb(1976), Eq(12).
     !--------------------------------------------------------
     br(i)=govrth(i)*za(i)*dthvdz/(wspd(i)*wspd(i))
-    if (itimestep == 1) then
+    if (.not. restart_or_cycle .and. itimestep == 1) then
        !set limits according to Li et al. (2010) boundary-layer meteorol (p.158)
        br(i)=max(br(i),-2.0)
        br(i)=min(br(i),2.0)
@@ -599,8 +600,8 @@ subroutine mynnsfclay(                                 &
           regime(i)=2.
        endif
 
-       !compute z/l first guess:
-       if (itimestep .le. 1) then
+       if (.not. restart_or_cycle .or. (restart_or_cycle .and. itimestep > 1) ) then
+          !compute z/l first guess:
           call li_etal_2010(zol(i),br(i),za(i)/zntstoch(i),zratio(i))
        else
           zol(i)=za(i)*karman*g*mol(i)/(th1d(i)*max(ust(i)*ust(i),0.0001))
@@ -684,8 +685,8 @@ subroutine mynnsfclay(                                 &
        !==========================================================
        regime(i)=4.
 
-       !compute z/l first guess:
-       if (itimestep .le. 1) then
+       if (.not. restart_or_cycle .or. (restart_or_cycle .and. itimestep > 1) ) then
+          !compute z/l first guess:
           call li_etal_2010(zol(i),br(i),za(i)/zntstoch(i),zratio(i))
        else
           zol(i)=za(i)*karman*g*mol(i)/(th1d(i)*max(ust(i)*ust(i),0.001))