diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 6378596797..e32ce6df48 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -72,6 +72,11 @@ description="Model time step, seconds" possible_values="Positive real values"/> + + + + domain % clock block => domain % blocklist @@ -667,7 +671,7 @@ subroutine atm_timestep(domain, dt, nowTime, itimestep, exchange_halo_group) call mpas_pool_get_config(block % configs, 'config_apply_lbcs', config_apply_lbcs) if (trim(config_time_integration) == 'SRK3') then - call atm_srk3(domain, dt, itimestep, exchange_halo_group) + call atm_srk3(domain, dt, itimestep, xtime_s, exchange_halo_group) else call mpas_log_write('Unknown time integration option '//trim(config_time_integration), messageType=MPAS_LOG_ERR) call mpas_log_write('Currently, only ''SRK3'' is supported.', messageType=MPAS_LOG_CRIT) @@ -692,7 +696,7 @@ subroutine atm_timestep(domain, dt, nowTime, itimestep, exchange_halo_group) end subroutine atm_timestep - subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) + subroutine atm_srk3(domain, dt, itimestep, xtime_s, exchange_halo_group) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Advance model state forward in time by the specified time step using ! time-split RK3 scheme @@ -711,6 +715,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) type (domain_type), intent(inout) :: domain real (kind=RKIND), intent(in) :: dt integer, intent(in) :: itimestep + real (kind=RKIND), intent(in) :: xtime_s procedure (halo_exchange_routine) :: exchange_halo_group integer :: thread @@ -742,7 +747,6 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) logical, pointer :: config_scalar_advection logical, pointer :: config_positive_definite logical, pointer :: config_monotonic - real (kind=RKIND), pointer :: config_dt character (len=StrKIND), pointer :: config_microp_scheme character (len=StrKIND), pointer :: config_convection_scheme @@ -781,8 +785,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) call mpas_pool_get_config(block % configs, 'config_time_integration_order', config_time_integration_order) call mpas_pool_get_config(block % configs, 'config_scalar_advection', config_scalar_advection) call mpas_pool_get_config(block % configs, 'config_positive_definite', config_positive_definite) - call mpas_pool_get_config(block % configs, 'config_monotonic', config_monotonic) - call mpas_pool_get_config(block % configs, 'config_dt', config_dt) + call mpas_pool_get_config(block % configs, 'config_monotonic', config_monotonic) call mpas_pool_get_config(block % configs, 'config_IAU_option', config_IAU_option) ! config variables for dynamics-transport splitting, WCS 18 November 2014 call mpas_pool_get_config(block % configs, 'config_split_dynamics_transport', config_split_dynamics_transport) @@ -870,6 +873,8 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) ! Initialize RK weights ! + call mpas_log_write('--- atm_srk3 time-step DT = $r',realArgs=(/dt/)) + dynamics_split = config_dynamics_split if (config_split_dynamics_transport) then dt_dynamics = dt/real(dynamics_split) @@ -974,7 +979,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) ! IAU - Incremental Analysis Update ! if (trim(config_IAU_option) /= 'off') then - call atm_add_tend_anal_incr(block % configs, block % structs, itimestep, dt, & + call atm_add_tend_anal_incr(block % configs, block % structs, itimestep, dt, xtime_s, & tend_ru_physics, tend_rtheta_physics, tend_rho_physics) end if @@ -1482,7 +1487,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) !requires that the subroutine atm_advance_scalars_mono was called on the third Runge Kutta step, so that a halo !update for the scalars at time_levs(1) is applied. A halo update for the scalars at time_levs(2) is done above. if (config_monotonic) then - rqvdynten(:,:) = ( scalars_2(index_qv,:,:) - scalars_1(index_qv,:,:) ) / config_dt + rqvdynten(:,:) = ( scalars_2(index_qv,:,:) - scalars_1(index_qv,:,:) ) / dt else rqvdynten(:,:) = 0._RKIND end if diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 997d7ca8ba..c6d0f55f66 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -30,6 +30,8 @@ end subroutine halo_exchange_routine end interface type (MPAS_Clock_type), pointer :: clock + type (MPAS_TimeInterval_type) :: dt_incr, timeStep_init + type (MPAS_Time_Type) :: rampup_stopTime contains @@ -55,7 +57,7 @@ function atm_core_init(domain, startTimeStamp) result(ierr) character(len=*), intent(out) :: startTimeStamp integer :: ierr - real (kind=RKIND), pointer :: dt + real (kind=RKIND) :: dt type (block_type), pointer :: block logical, pointer :: config_do_restart @@ -68,6 +70,7 @@ function atm_core_init(domain, startTimeStamp) result(ierr) character (len=StrKIND), pointer :: xtime character (len=StrKIND), pointer :: initial_time1, initial_time2 type (MPAS_Time_Type) :: startTime + type (MPAS_TimeInterval_type) :: timeStepInterval !< the current time step as an interval real (kind=RKIND), pointer :: nominalMinDc real (kind=RKIND), pointer :: config_len_disp @@ -116,8 +119,9 @@ 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_dt', dt) + call mpas_get_timeInterval(mpas_get_clock_timestep(domain % clock, ierr), dt=dt) + call mpas_log_write('atm_core_init: Obtaining initial dt=$r', realArgs=(/dt/)) ! ! By default, the 'invariant' stream has an input_interval of "none", so ! the following stream read has no effect. However, in case the 'invariant' @@ -309,12 +313,18 @@ subroutine atm_simulation_clock_init(core_clock, configs, ierr) integer, intent(out) :: ierr type (MPAS_Time_Type) :: startTime, stopTime - type (MPAS_TimeInterval_type) :: runDuration, timeStep + type (MPAS_TimeInterval_type) :: runDuration, timeStep_config, timeStep, timeStep_sum + type (MPAS_TimeInterval_type) :: timestep_rampup_Duration, timestep_rampup_Duration_2, searchRemainder, tmp_interval integer :: local_err + integer (kind=I8KIND) :: nDivs + !integer :: nDivs + real (kind=RKIND) :: dt_incr_tmp, loc_rampup_duration, nDivs_r, loc_dt_incr, dt_rem, tot_remainder, tmp_b real (kind=RKIND), pointer :: config_dt + real (kind=RKIND), pointer :: config_dt_init character (len=StrKIND), pointer :: config_start_time character (len=StrKIND), pointer :: config_restart_timestamp_name character (len=StrKIND), pointer :: config_run_duration + character (len=StrKIND), pointer :: config_dt_rampup_duration character (len=StrKIND), pointer :: config_stop_time character (len=StrKIND) :: startTimeStamp integer :: iounit @@ -323,9 +333,11 @@ subroutine atm_simulation_clock_init(core_clock, configs, ierr) ierr = 0 call mpas_pool_get_config(configs, 'config_dt', config_dt) + call mpas_pool_get_config(configs, 'config_dt_init', config_dt_init) call mpas_pool_get_config(configs, 'config_start_time', config_start_time) call mpas_pool_get_config(configs, 'config_restart_timestamp_name', config_restart_timestamp_name) call mpas_pool_get_config(configs, 'config_run_duration', config_run_duration) + call mpas_pool_get_config(configs, 'config_dt_rampup_duration', config_dt_rampup_duration) call mpas_pool_get_config(configs, 'config_stop_time', config_stop_time) if(trim(config_start_time) == 'file') then @@ -338,7 +350,64 @@ subroutine atm_simulation_clock_init(core_clock, configs, ierr) startTimeStamp = config_start_time end if call mpas_set_time(curr_time=startTime, dateTimeString=startTimeStamp, ierr=local_err) - call mpas_set_timeInterval(timeStep, dt=config_dt, ierr=local_err) + + !call mpas_set_timeInterval(timeStep, dt=config_dt, ierr=local_err) + + call mpas_set_timeInterval(timeStep_config, dt=config_dt, ierr=local_err) + timeStep = timeStep_config + if ( config_dt_init > 0.0_RKIND ) then + call mpas_log_write('config_dt_init if clause met') + if ( config_dt_init > config_dt ) then + call mpas_log_write('config_dt_init should be smaller than config_dt.', messageType=MPAS_LOG_ERR) + ierr = 1 + end if + call mpas_set_timeInterval(timeStep_init, dt=config_dt_init, ierr=local_err) + if (trim(config_dt_rampup_duration) /= "none") then + + call mpas_set_timeInterval(timestep_rampup_Duration, timeString=config_dt_rampup_duration, ierr=local_err) + + call mpas_log_write('timestep_rampup_Duration is set to'//trim(config_dt_rampup_duration)) + timestep_rampup_Duration_2 = timestep_rampup_Duration * 2 + timeStep_sum = timeStep_init + timeStep_config + + ! Assuming an arithmetic sequence of time step sizes, the sum of elements in the sequence + ! is given by the formula Sn = nDivs/2(2*config_dt_init + (nDivs-1)dt_incr), where config_dt_init + ! is the first time step size, dt_incr is the timestep size increment, and nDivs is the number of + ! elements in the sequence. The sum of time step sizes should ideally equal the rampup duration, + ! except for a remainder searchRemainder. + + ! Given the above assumption, we can solve for nDivs as + ! nDivs_r = floor(2.0_RKIND*loc_rampup_duration/(config_dt + config_dt_init)) - 1 + call mpas_interval_division(startTime, timestep_rampup_Duration_2, timeStep_sum, nDivs, searchRemainder) + nDivs = nDivs - 1 + ! set the time step increment + dt_incr = (timeStep_config - timeStep_init)/int(nDivs) + ! calculate the remainder, which will be set as the first time step size, in order to + ! correctly reach the rampup duration + searchRemainder = timestep_rampup_Duration - (timeStep_init*2 + dt_incr*(nDivs))*(nDivs+1)/2 + + + ! nDivs_r = floor(2.0_RKIND*loc_rampup_duration/(config_dt + config_dt_init)) - 1 + ! loc_dt_incr = (config_dt - config_dt_init)/nDivs_r + ! tot_remainder = (loc_rampup_duration- 0.5*(nDivs_r+1)*(2.0*config_dt_init + nDivs_r*loc_dt_incr)) + ! call mpas_log_write('nDivs :$r, loc_dt_incr is set to $r', & + ! realArgs=(/nDivs_r, loc_dt_incr/)) + ! call mpas_set_timeInterval(searchRemainder, dt=tot_remainder, ierr=local_err) + + call mpas_get_timeInterval(dt_incr, dt=dt_incr_tmp, ierr=local_err) + call mpas_get_timeInterval(searchRemainder, dt=dt_rem, ierr=local_err) + call mpas_log_write('dt_rem is set to $r, dt_incr is set to $r, nDivs: $i', & + realArgs=(/dt_rem, dt_incr_tmp/), intArgs=(/int(nDivs)/)) + + else + call mpas_log_write('config_dt_init is set but config_dt_rampup_duration was not specified.', messageType=MPAS_LOG_ERR) + ierr = 1 + end if + !timeStep = timeStep_init + ! Setting the first timestep to the remainder + timeStep = searchRemainder + rampup_stopTime = startTime + timestep_rampup_Duration + end if if (trim(config_run_duration) /= "none") then call mpas_set_timeInterval(runDuration, timeString=config_run_duration, ierr=local_err) @@ -602,12 +671,15 @@ function atm_core_run(domain) result(ierr) type (domain_type), intent(inout) :: domain integer :: ierr - real (kind=RKIND), pointer :: dt + real (kind=RKIND) :: dt + real (kind=RKIND), pointer :: config_dt logical, pointer :: config_do_restart logical, pointer :: config_apply_lbcs type (block_type), pointer :: block_ptr type (MPAS_Time_Type) :: currTime + integer :: i_dt_incr = 0 + character(len=StrKIND) :: timeStamp character (len=StrKIND), pointer :: config_restart_timestamp_name integer :: itimestep @@ -623,16 +695,25 @@ function atm_core_run(domain) result(ierr) real (kind=R8KIND) :: diag_start_time, diag_stop_time real (kind=R8KIND) :: input_start_time, input_stop_time real (kind=R8KIND) :: output_start_time, output_stop_time - + type (MPAS_TimeInterval_type) :: timeStepInterval !< the current time step as an interval + ierr = 0 clock => domain % clock mpas_log_info => domain % logInfo + + call mpas_get_timeInterval(mpas_get_clock_timestep(clock, ierr), dt=dt) + call mpas_log_write('Obtaining initial dt=$r', realArgs=(/dt/)) + ! call mpas_log_write('Setting dt=60') + !call mpas_set_timeInterval(timeStepInterval, dt=60.0_RKIND, ierr=ierr) + !call mpas_set_clock_timestep(clock, timeStepInterval, ierr) + - ! Eventually, dt should be domain specific - call mpas_pool_get_config(domain % blocklist % configs, 'config_dt', dt) call mpas_pool_get_config(domain % blocklist % configs, 'config_do_restart', config_do_restart) call mpas_pool_get_config(domain % blocklist % configs, 'config_restart_timestamp_name', config_restart_timestamp_name) + !call mpas_pool_get_config(domain % blocklist % configs, 'config_dt', config_dt) + + !dt = config_dt ! 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) @@ -725,7 +806,10 @@ function atm_core_run(domain) result(ierr) do while (.not. mpas_is_clock_stop_time(clock)) currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr) - call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr) + call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr) + + !call mpas_get_timeInterval(mpas_get_clock_timestep(clock, ierr), dt=dt) + !call mpas_log_write('Obtaining dt=$r', realArgs=(/dt/)) call mpas_log_write('') call mpas_log_write('Begin timestep '//trim(timeStamp)) @@ -879,6 +963,24 @@ function atm_core_run(domain) result(ierr) call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_OUTPUT, ierr=ierr) + if (currTime < rampup_stopTime) then + + timeStepInterval = mpas_get_clock_timestep(clock, ierr) + call mpas_get_timeInterval(timeStepInterval, dt=dt) + call mpas_log_write('Obtaining current dt=$r', realArgs=(/dt/)) + !dt = real(floor(min(dt * (1.05_RKIND),config_dt))) + if (i_dt_incr == 0) then + timeStepInterval = timeStep_init + else + timeStepInterval = timeStepInterval + dt_incr + end if + i_dt_incr = i_dt_incr + 1 + call mpas_get_timeInterval(timeStepInterval, dt=dt) + call mpas_log_write('Incrementing dt to $r, incr count: $i', realArgs=(/dt/), intArgs=(/i_dt_incr/)) + !call mpas_set_timeInterval(timeStepInterval, dt=dt, ierr=ierr) + call mpas_set_clock_timestep(clock, timeStepInterval, ierr) + end if + end do end function atm_core_run @@ -1001,7 +1103,7 @@ subroutine atm_do_timestep(domain, dt, itimestep) !proceed with physics if moist_physics is set to true: if(moist_physics) then call physics_timetracker(domain,dt,clock,itimestep,xtime_s) - call physics_driver(domain,itimestep,xtime_s) + call physics_driver(domain,dt,itimestep,xtime_s) endif #endif diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver.F b/src/core_atmosphere/physics/mpas_atmphys_driver.F index 402517e98f..8ea1650297 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver.F @@ -107,10 +107,11 @@ module mpas_atmphys_driver !================================================================================================================= - subroutine physics_driver(domain,itimestep,xtime_s) + subroutine physics_driver(domain,dt,itimestep,xtime_s) !================================================================================================================= !input arguments: + real(kind=RKIND),intent(in):: dt integer,intent(in):: itimestep real(kind=RKIND),intent(in):: xtime_s @@ -282,7 +283,7 @@ subroutine physics_driver(domain,itimestep,xtime_s) call allocate_lsm !$OMP PARALLEL DO do thread=1,nThreads - call driver_lsm(itimestep,block%configs,mesh,diag_physics,sfc_input, & + call driver_lsm(itimestep,xtime_s,block%configs,mesh,diag_physics,sfc_input, & cellSolveThreadStart(thread),cellSolveThreadEnd(thread)) end do !$OMP END PARALLEL DO @@ -290,7 +291,7 @@ subroutine physics_driver(domain,itimestep,xtime_s) elseif(config_lsm_scheme == 'sf_noahmp') then do thread=1,nThreads - call driver_lsm_noahmp(block%configs,mesh,state,time_lev,diag,diag_physics, & + call driver_lsm_noahmp(block%configs,mesh,state,dt,time_lev,diag,diag_physics, & diag_physics_noahmp,output_noahmp,sfc_input,itimestep, & cellSolveThreadStart(thread),cellSolveThreadEnd(thread)) enddo diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F index 0116dcf561..a638e6fd7b 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F @@ -664,7 +664,7 @@ subroutine init_lsm(dminfo,mesh,configs,diag_physics,sfc_input) end subroutine init_lsm !================================================================================================================= - subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite) + subroutine driver_lsm(itimestep,xtime_s,configs,mesh,diag_physics,sfc_input,its,ite) !================================================================================================================= !input arguments: @@ -673,6 +673,7 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite) integer,intent(in):: its,ite integer,intent(in):: itimestep + real(kind=RKIND),intent(in):: xtime_s !inout arguments: type(mpas_pool_type),intent(inout):: diag_physics @@ -729,7 +730,7 @@ subroutine driver_lsm(itimestep,configs,mesh,diag_physics,sfc_input,its,ite) fvb_2d = fvbsnow_p , fbur_2d = fbursnow_p , fgsn_2d = fgsnsnow_p , & utype_urb2d = utype_urb_p , frc_urb2d = frc_urb_p , ust_urb2d = ust_urb_p , & swddir = swddir_p , swddif = swddif_p , fasdas = fasdas , & - julian = 0 , julyr = 0 , & + julian = 0 , julyr = 0 , xtime_s = xtime_s , & num_soil_layers = num_soils , & xice_threshold = xice_threshold , & usemonalb = config_sfc_albedo , & diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F index 8bbec89911..4d556c0662 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F @@ -30,7 +30,7 @@ module mpas_atmphys_driver_lsm_noahmp !================================================================================================================= subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input, & - state,time_lev,itimestep) + state,dt,time_lev,itimestep) !================================================================================================================= !--- input arguments: @@ -39,6 +39,7 @@ subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahm type(mpas_pool_type),intent(in):: diag type(mpas_pool_type),intent(in):: state + real(kind=RKIND),intent(in):: dt integer,intent(in):: time_lev integer,intent(in):: itimestep @@ -127,10 +128,11 @@ subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahm mpas_noahmp%day = day mpas_noahmp%julian = curr_julday - + mpas_noahmp%dtbl = dt !--- initialization of xice_threshold: mpas_noahmp%xice_threshold = xice_threshold + call mpas_log_write('--- lsm_noahmp_fromMPAS dtbl = $r',realArgs=(/dt/)) !--- initialization of INPUT surface variables: call mpas_pool_get_array(sfc_input,'shdmax',shdmax) @@ -1054,7 +1056,7 @@ subroutine lsm_noahmp_toMPAS(diag_physics,diag_physics_noahmp,output_noahmp,sfc_ end subroutine lsm_noahmp_toMPAS !================================================================================================================= - subroutine driver_lsm_noahmp(configs,mesh,state,time_lev,diag,diag_physics,diag_physics_noahmp,output_noahmp, & + subroutine driver_lsm_noahmp(configs,mesh,state,dt,time_lev,diag,diag_physics,diag_physics_noahmp,output_noahmp, & sfc_input,itimestep,its,ite) !================================================================================================================= @@ -1064,6 +1066,7 @@ subroutine driver_lsm_noahmp(configs,mesh,state,time_lev,diag,diag_physics,diag_ type(mpas_pool_type),intent(in):: diag type(mpas_pool_type),intent(in):: state + real(kind=RKIND),intent(in):: dt integer,intent(in):: itimestep,its,ite integer,intent(in):: time_lev @@ -1079,7 +1082,7 @@ subroutine driver_lsm_noahmp(configs,mesh,state,time_lev,diag,diag_physics,diag_ !call mpas_log_write('--- enter subroutine driver_lsm_noahmp:') call lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input, & - state,time_lev,itimestep) + state,dt,time_lev,itimestep) call NoahmpDriverMain(mpas_noahmp) diff --git a/src/core_atmosphere/physics/mpas_atmphys_init.F b/src/core_atmosphere/physics/mpas_atmphys_init.F index ce767fb6b3..15f9884c84 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_init.F +++ b/src/core_atmosphere/physics/mpas_atmphys_init.F @@ -377,7 +377,7 @@ subroutine physics_init(dminfo,clock,configs,mesh,diag,tend,state,time_lev,diag_ if(config_lsm_scheme .eq. 'sf_noah') then call init_lsm(dminfo,mesh,configs,diag_physics,sfc_input) elseif(config_lsm_scheme .eq. 'sf_noahmp') then - call init_lsm_noahmp(configs,mesh,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) + call init_lsm_noahmp(configs,mesh,clock,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) endif endif diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F index da1cead2c0..922bd98645 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F +++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F @@ -9,7 +9,7 @@ module mpas_atmphys_lsm_noahmpinit use mpas_log use mpas_pool_routines - + use mpas_timekeeping, only: mpas_get_timeInterval, mpas_get_clock_timestep use mpas_atmphys_utilities,only: physics_error_fatal use mpas_atmphys_vars,only : mpas_noahmp @@ -28,12 +28,13 @@ module mpas_atmphys_lsm_noahmpinit !================================================================================================================= - subroutine init_lsm_noahmp(configs,mesh,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) + subroutine init_lsm_noahmp(configs,mesh,clock,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) !================================================================================================================= !--- input arguments: type(mpas_pool_type),intent(in):: configs type(mpas_pool_type),intent(in):: mesh + type(MPAS_Clock_type),intent(in):: clock !--- inout arguments: type(mpas_pool_type),intent(inout):: diag_physics @@ -88,7 +89,7 @@ subroutine init_lsm_noahmp(configs,mesh,diag_physics,diag_physics_noahmp,output_ !--- initialize noahmp: - call noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) + call noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) !call mpas_log_write('--- end subroutine init_lsm_noahmp:') @@ -213,12 +214,13 @@ subroutine noahmp_read_namelist(configs) end subroutine noahmp_read_namelist !================================================================================================================= - subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) + subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,output_noahmp,sfc_input) !================================================================================================================= !--- input arguments: type(mpas_pool_type),intent(in):: configs type(mpas_pool_type),intent(in):: mesh + type(MPAS_Clock_type),intent(in):: clock !--- inout arguments: type(mpas_pool_type),intent(inout):: diag_physics @@ -234,7 +236,7 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah integer,dimension(:),pointer:: isnowxy integer,dimension(:),pointer:: irnumsi,irnummi,irnumfi - real(kind=RKIND),pointer:: dt + real(kind=RKIND) :: dt real(kind=RKIND),dimension(:),pointer:: soilcl1,soilcl2,soilcl3,soilcl4 real(kind=RKIND),dimension(:,:),pointer:: soilcomp @@ -271,7 +273,8 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah !--- initialization of Noah-MP run parameters: call mpas_pool_get_config(configs,'config_do_restart',do_restart) call mpas_pool_get_config(configs,'config_urban_physics',urban_physics) - call mpas_pool_get_config(configs,'config_dt',dt) + call mpas_get_timeInterval(mpas_get_clock_timestep(clock, ierr), dt=dt) + call mpas_log_write('--- LSM time-step DT = $r',realArgs=(/dt/)) mpas_noahmp%restart_flag = do_restart mpas_noahmp%sf_urban_physics = 0 diff --git a/src/core_atmosphere/physics/mpas_atmphys_manager.F b/src/core_atmosphere/physics/mpas_atmphys_manager.F index 1b88522786..25dbed881e 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_manager.F +++ b/src/core_atmosphere/physics/mpas_atmphys_manager.F @@ -156,10 +156,12 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) character(len=StrKIND),pointer:: config_convection_scheme, & config_radt_lw_scheme, & + config_microp_scheme, & config_radt_sw_scheme character(len=StrKIND),pointer:: config_conv_interval, & config_radtlw_interval, & + config_pbl_interval, & config_radtsw_interval type(block_type),pointer :: block @@ -187,9 +189,11 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) call mpas_pool_get_config(domain%blocklist%configs,'config_convection_scheme',config_convection_scheme) call mpas_pool_get_config(domain%blocklist%configs,'config_radt_lw_scheme' ,config_radt_lw_scheme ) call mpas_pool_get_config(domain%blocklist%configs,'config_radt_sw_scheme' ,config_radt_sw_scheme ) + call mpas_pool_get_config(domain%blocklist%configs,'config_microp_scheme' ,config_microp_scheme ) call mpas_pool_get_config(domain%blocklist%configs,'config_conv_interval' ,config_conv_interval ) call mpas_pool_get_config(domain%blocklist%configs,'config_radtlw_interval',config_radtlw_interval) + call mpas_pool_get_config(domain%blocklist%configs,'config_pbl_interval',config_pbl_interval ) call mpas_pool_get_config(domain%blocklist%configs,'config_radtsw_interval',config_radtsw_interval) call mpas_pool_get_config(domain%blocklist%configs,'config_frac_seaice' ,config_frac_seaice ) @@ -219,6 +223,22 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) !call mpas_log_write(' LEAP_YEAR = $l', logicArgs=(/LeapYear/)) !call mpas_log_write(' TIME STAMP = '//trim(timeStamp)) + dt_dyn = dt +!... cloud microphysics: + dt_microp = dt_dyn + n_microp = 1 + if(trim(config_microp_scheme)=='mp_thompson' .or. & + trim(config_microp_scheme)=='mp_thompson_aerosols') then + dt_microp = 90._RKIND + n_microp = max(nint(dt_dyn/dt_microp),1) + dt_microp = dt_dyn / n_microp + if(dt_dyn <= dt_microp) dt_microp = dt_dyn + endif + call mpas_log_write(' ') + call mpas_log_write('--- dt_microp = $r', realArgs=(/dt_microp/)) + call mpas_log_write('--- n_microp = $i', intArgs=(/n_microp/)) + + block => domain % blocklist do while(associated(block)) @@ -267,6 +287,8 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) endif elseif(config_radtlw_interval == "none") then l_radtlw = .true. + dt_radtlw = dt + call mpas_log_write('---- dt_radtlw = $r',realArgs=(/dt_radtlw/)) endif call mpas_log_write('--- time to run the LW radiation scheme L_RADLW = $l',logicArgs=(/l_radtlw/)) endif @@ -281,6 +303,8 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) endif elseif(config_radtsw_interval == "none") then l_radtsw = .true. + dt_radtsw = dt + call mpas_log_write('---- dt_radtsw = $r',realArgs=(/dt_radtsw/)) endif call mpas_log_write('--- time to run the SW radiation scheme L_RADSW = $l',logicArgs=(/l_radtsw/)) endif @@ -296,10 +320,19 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) endif elseif(config_conv_interval == "none") then l_conv = .true. + dt_cu = dt + n_cu = nint(dt_cu/dt_dyn) + n_cu = max(n_cu,1) + call mpas_log_write('---- dt_cu = $r',realArgs=(/dt_cu/)) endif call mpas_log_write('--- time to run the convection scheme L_CONV = $l',logicArgs=(/l_conv/)) endif + if(trim(config_pbl_interval) == "none") then + dt_pbl = dt + call mpas_log_write('---- dt_pbl = $r',realArgs=(/dt_pbl/)) + endif + !check to see if it is time to update ozone to the current julian day in the RRTMG radiation codes: if(config_o3climatology) then block => domain % blocklist @@ -410,7 +443,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) integer,pointer:: nAerosols,nAerLevels,nOznLevels integer,pointer:: nCellsSolve,nSoilLevels,nVertLevels - real(kind=RKIND),pointer:: config_dt + real(kind=RKIND):: dt !local variables: type(MPAS_Time_Type):: startTime,alarmStartTime @@ -439,8 +472,8 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call mpas_pool_get_config(configs,'config_frac_seaice' ,config_frac_seaice ) call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re ) - call mpas_pool_get_config(configs,'config_dt',config_dt) - + call mpas_get_timeInterval(mpas_get_clock_timestep(clock, ierr), dt=dt) + call mpas_log_write('--- physics time-step DT = $r',realArgs=(/dt/)) call mpas_pool_get_dimension(mesh,'cam_dim1' ,cam_dim1 ) call mpas_pool_get_dimension(mesh,'nMonths' ,nMonths ) @@ -482,7 +515,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call physics_error_fatal('subroutine physics_run_init: error defining dt_radtlw') elseif(trim(config_radtlw_interval) == "none") then - dt_radtlw = config_dt + dt_radtlw = dt else call physics_error_fatal('subroutine physics_run_init: dt_radtlw is not defined') @@ -501,7 +534,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call physics_error_fatal('subroutine physics_run_init: error defining radtswAlarmID') elseif(trim(config_radtsw_interval) == "none") then - dt_radtsw = config_dt + dt_radtsw = dt else call physics_error_fatal('subroutine physics_run_init: dt_radtsw is not defined') @@ -520,7 +553,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call physics_error_fatal('subroutine physics_run_init: error defining dt_cu') elseif(trim(config_conv_interval) == "none") then - dt_cu = config_dt + dt_cu = dt else call physics_error_fatal('subroutine physics_run_init: dt_cu is not defined') @@ -539,7 +572,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call physics_error_fatal('subroutine physics_run_init: error defining dt_pbl') elseif(trim(config_pbl_interval) == "none") then - dt_pbl = config_dt + dt_pbl = dt else call physics_error_fatal('subroutine physics_run_init: dt_pbl is not defined') @@ -578,7 +611,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) !set alarm to write the "CAM" local arrays absnst_p, absnxt_p, and emstot_p to the MPAS arrays !for writing to the restart file at the bottom of the time-step: if(trim(config_radt_lw_scheme) .eq. "cam_lw" ) then - call mpas_set_timeInterval(camlwTimeStep,dt=config_dt,ierr=ierr) + call mpas_set_timeInterval(camlwTimeStep,dt=dt,ierr=ierr) call MPAS_stream_mgr_get_property(stream_manager, 'restart', MPAS_STREAM_PROPERTY_RECORD_INTV, stream_interval, & direction=MPAS_STREAM_OUTPUT, ierr=ierr) if(trim(stream_interval) /= 'none') then @@ -593,7 +626,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) !set alarm to check if the accumulated rain due to cloud microphysics and convection is !greater than its maximum allowed value: if(config_bucket_update /= "none") then - call mpas_set_timeInterval(acrainTimeStep,dt=config_dt,ierr=ierr) + call mpas_set_timeInterval(acrainTimeStep,dt=dt,ierr=ierr) call mpas_set_timeInterval(alarmTimeStep,timeString=config_bucket_update,ierr=ierr) alarmStartTime = startTime + alarmTimeStep call mpas_add_clock_alarm(clock,acrainAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) @@ -604,7 +637,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) !set alarm to check if the accumulated radiation diagnostics due to long- and short-wave radiation !is greater than its maximum allowed value: if(config_bucket_update /= "none") then - call mpas_set_timeInterval(acradtTimeStep,dt=config_dt,ierr=ierr) + call mpas_set_timeInterval(acradtTimeStep,dt=dt,ierr=ierr) call mpas_set_timeInterval(alarmTimeStep,timeString=config_bucket_update,ierr=ierr) alarmStartTime = startTime + alarmTimeStep call mpas_add_clock_alarm(clock,acradtAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) @@ -681,7 +714,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) !initialization of local physics time-steps: !... dynamics: - dt_dyn = config_dt + dt_dyn = dt !... cloud microphysics: dt_microp = dt_dyn n_microp = 1 @@ -694,6 +727,7 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) endif call mpas_log_write(' ') call mpas_log_write('--- specifics on cloud microphysics option microp_scheme = '//trim(config_microp_scheme)) + call mpas_log_write('--- dt_dyn = $r', realArgs=(/dt_dyn/)) call mpas_log_write('--- dt_microp = $r', realArgs=(/dt_microp/)) call mpas_log_write('--- n_microp = $i', intArgs=(/n_microp/)) !... convection: diff --git a/src/core_atmosphere/physics/physics_wrf/module_sf_bep.F b/src/core_atmosphere/physics/physics_wrf/module_sf_bep.F index a0f2bdf645..27693a87a7 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_sf_bep.F +++ b/src/core_atmosphere/physics/physics_wrf/module_sf_bep.F @@ -64,7 +64,7 @@ MODULE module_sf_bep CONTAINS - subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, & + subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,xtime_s,dz8w,dt,u_phy,v_phy, & th_phy,rho,p_phy,swdown,glw, & gmt,julday,xlong,xlat, & declin_urb,cosz_urb2d,omg_urb2d, & @@ -141,6 +141,7 @@ subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, & ! integer nx,ny,nz ! Number of points in the mesocsale grid real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates REAL, INTENT(IN ):: DT ! Time step + REAL, INTENT(IN ):: xtime_s ! Time step ! real zr(ims:ime,jms:jme) ! Solar zenith angle ! real deltar(ims:ime,jms:jme) ! Declination of the sun ! real ah(ims:ime,jms:jme) ! Hour angle @@ -488,7 +489,8 @@ subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, & rs1D=swdown(ix,iy) rld1D=glw(ix,iy) - time_h=(itimestep*dt)/3600.+gmt + !time_h=(itimestep*dt)/3600.+gmt + time_h=xtime_s/3600.+gmt zr1D=acos(COSZ_URB2D(ix,iy)) deltar1D=DECLIN_URB diff --git a/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F b/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F index ab29ea29b8..069a9b5087 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F +++ b/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F @@ -45,7 +45,7 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,ISURBAN,ISICE,VEGFRA, & ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,EMISS,EMBCK, & SNOWC,QSFC,RAINBL,MMINLU, & - num_soil_layers,DT,DZS,ITIMESTEP, & + num_soil_layers,DT,DZS,ITIMESTEP, XTIME_S, & SMOIS,TSLB,SNOW,CANWAT, & CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,lai,qz0, & !H myj,frpcpn, & @@ -323,6 +323,8 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP + REAL, INTENT(IN):: XTIME_S + REAL, INTENT(IN ) :: DT,ROVCP REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS @@ -1604,7 +1606,7 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & b_q_bep(i,kts:kte,j)=0. end do end do - CALL BEP(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, & + CALL BEP(frc_urb2d,utype_urb2d,itimestep,xtime_s,dz8w,dt,u_phy,v_phy, & th_phy,rho,p_phy,swdown,glw, & gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, & num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &