diff --git a/components/mpas-framework/src/framework/mpas_forcing.F b/components/mpas-framework/src/framework/mpas_forcing.F index 950b103604c7..95de2ca61846 100644 --- a/components/mpas-framework/src/framework/mpas_forcing.F +++ b/components/mpas-framework/src/framework/mpas_forcing.F @@ -35,7 +35,8 @@ module mpas_forcing mpas_forcing_init_field_data, & mpas_forcing_get_forcing, & mpas_forcing_get_forcing_time, & - mpas_forcing_write_restart_times + mpas_forcing_write_restart_times, & + mpas_advance_forcing_clock contains @@ -1193,7 +1194,7 @@ subroutine mpas_forcing_get_forcing(&!{{{ if (trim(forcingGroup % forcingGroupName) == trim(forcingGroupName)) then ! advance the forcing time - call advance_forcing_clock(forcingGroup, dt) + call mpas_advance_forcing_clock(forcingGroup, dt) ! cycle the forcing clock if (forcingGroup % forcingCycleUse) then @@ -1230,7 +1231,7 @@ end subroutine mpas_forcing_get_forcing!}}} !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! advance_forcing_clock +! mpas_advance_forcing_clock ! !> \brief set the forcing clock !> \author Adrian K. Turner, LANL @@ -1241,7 +1242,7 @@ end subroutine mpas_forcing_get_forcing!}}} ! !----------------------------------------------------------------------- - subroutine advance_forcing_clock(&!{{{ + subroutine mpas_advance_forcing_clock(&!{{{ forcingGroup, & dt) @@ -1256,9 +1257,10 @@ subroutine advance_forcing_clock(&!{{{ ! increment clock with timestep call mpas_set_timeInterval(timeStep, dt=dt) + !call mpas_log_write('time step from mpas_set_timeInterval: $i $i $i', intArgs=(/int(timeStep%ti%basetime%S), int(timeStep%ti%basetime%Sn), int(timeStep%ti%basetime%Sd) /)) call mpas_advance_clock(forcingGroup % forcingClock, timeStep) - end subroutine advance_forcing_clock!}}} + end subroutine mpas_advance_forcing_clock!}}} !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -1597,8 +1599,10 @@ subroutine get_interpolants_linear(interpolants, forcingStream, currentTime)!{{{ call mpas_get_timeInterval(diff1, forcingStream % forcingTimes(1), dt=diffr1) call mpas_get_timeInterval(diff2, currentTime, dt=diffr2) + !call mpas_log_write('diffr2 $r, diffr, $r', realArgs=(/ diffr2, diffr /)) interpolants(1) = diffr2 / diffr interpolants(2) = 1.0_RKIND - interpolants(1) !diffr1 / diffr + end subroutine get_interpolants_linear!}}} @@ -1612,7 +1616,7 @@ end subroutine get_interpolants_linear!}}} !> \details !> Given the current time and forcing times calculate the correct !> interpolants with piecewise constant interpolation -! + !----------------------------------------------------------------------- subroutine get_interpolants_constant(interpolants, forcingStream, currentTime)!{{{ diff --git a/components/mpas-framework/src/framework/mpas_timekeeping.F b/components/mpas-framework/src/framework/mpas_timekeeping.F index 37a24d77358d..57d9651335a4 100644 --- a/components/mpas-framework/src/framework/mpas_timekeeping.F +++ b/components/mpas-framework/src/framework/mpas_timekeeping.F @@ -1427,6 +1427,8 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8, timeString_ = timeString end if + !call mpas_log_write('timeString_ '//trim(timeString_)) + numerator = 0 denominator = 1 @@ -1519,6 +1521,12 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8, return end if + if (sec < 0) then + numerator = -numerator + end if + + !call mpas_log_write('sec $i, num $i, denom $i', intArgs=(/int(sec),numerator,denominator/)) + call ESMF_TimeIntervalSet(interval % ti, YY=year, MM=month, D=day, H=hour, M=min, S_i8=sec, Sn=numerator, Sd=denominator, rc=ierr) else diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 2a21e22f40d2..382344585d87 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -4130,6 +4130,10 @@ description="Multiplication factors to smooth ssh at coastlines for SAL caculation" packages="tidalPotentialForcingPKG" /> + domain % blocklist do while (associated(block)) - call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err ) + call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool) + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) + forcingTimeIncrement = forcingTimeIncrementRK4(rk_step) + block => block % next + end do + + call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock) + block => domain % blocklist + do while (associated(block)) + call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err ) call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err ) block => block % next end do diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index 9529fbdc6ea1..103dadda7523 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -380,6 +380,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew + real (kind=RKIND), pointer :: forcingTimeIncrement + ! Remap variables real (kind=RKIND), dimension(:,:), pointer :: & layerThicknessLagNew @@ -507,6 +509,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2) + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', & + forcingTimeIncrement) + + forcingTimeIncrement = 0.0_RKIND + allocate(bottomDepthEdge(nEdgesAll+1)) if (config_transport_tests_flow_id > 0) then diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 653d88198131..d60098a8aeaf 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -280,6 +280,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew + real (kind=RKIND), pointer :: forcingTimeIncrement + ! Remap variables real (kind=RKIND), dimension(:,:), pointer :: & layerThicknessLagNew @@ -398,6 +400,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ allocate(baroclinicThickness(nEdgesAll+1)) + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', & + forcingTimeIncrement) + + forcingTimeIncrement = 0.0_RKIND + if (config_transport_tests_flow_id > 0) then ! This is a transport test. Write advection velocity from prescribed ! flow field. diff --git a/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F b/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F index e9ad6786103c..32efb9bf8a8a 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F @@ -76,7 +76,13 @@ module ocn_manufactured_solution ! !----------------------------------------------------------------------- - subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{ + subroutine ocn_manufactured_solution_tend_thick(forcingPool, tend, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: forcingPool !----------------------------------------------------------------- ! input/output variables @@ -97,6 +103,7 @@ subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{ integer :: iCell, kmin, kmax, k real (kind=RKIND) :: phase, time + real (kind=RKIND), pointer :: forcingTimeIncrement ! End preamble !------------- @@ -104,7 +111,9 @@ subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{ if (.not. config_use_manufactured_solution) return - time = daysSinceStartOfSim*86400.0_RKIND + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) + + time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement do iCell = 1,nCellsOwned @@ -137,7 +146,13 @@ end subroutine ocn_manufactured_solution_tend_thick!}}} !> This routine computes the velocity tendency for the manufactured solution ! !----------------------------------------------------------------------- - subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{ + subroutine ocn_manufactured_solution_tend_vel(forcingPool, tend, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: forcingPool !----------------------------------------------------------------- ! input/output variables @@ -158,6 +173,7 @@ subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{ integer :: iEdge, kmin, kmax, k real (kind=RKIND) :: phase, u, v, time + real (kind=RKIND), pointer :: forcingTimeIncrement ! End preamble !----------------------------------------------------------------- @@ -165,7 +181,9 @@ subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{ if (.not. config_use_manufactured_solution) return - time = daysSinceStartOfSim*86400.0_RKIND + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) + + time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement do iEdge = 1, nEdgesOwned diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index f85427913b09..e1c1509a1cd2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -264,7 +264,8 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ tendThick, err) ! Compute thickness tendency to manufactured solution - call ocn_manufactured_solution_tend_thick(tendThick, err) + call ocn_manufactured_solution_tend_thick(forcingPool, & + tendThick, err) #ifdef MPAS_OPENACC !$acc exit data copyout(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff, & @@ -493,7 +494,8 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & layerThickEdgeMean, tendVel, err) ! Compute tendency term for manufactured solution - call ocn_manufactured_solution_tend_vel(tendVel, err) + call ocn_manufactured_solution_tend_vel(forcingPool, & + tendVel, err) ! vertical mixing treated implicitly in a later routine ! adjust total velocity tendency based on wetting/drying diff --git a/components/mpas-ocean/src/shared/mpas_ocn_time_varying_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_time_varying_forcing.F index 8e64618a27ad..2afa3064050f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_time_varying_forcing.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_time_varying_forcing.F @@ -258,13 +258,21 @@ subroutine atmospheric_forcing(streamManager, domain, simulationClock)!{{{ character(len=StrKIND) :: timeStamp - real (kind=RKIND) :: dtSim + real (kind=RKIND) :: dtSim, dtSimReverse + + real (kind=RKIND), pointer :: forcingTimeIncrement + + type (mpas_pool_type), pointer :: forcingPool integer :: err + call mpas_pool_get_subpool(domain % blocklist % structs, 'forcing', forcingPool) + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) + ! convert config_dt to real call mpas_set_timeInterval(timeStepESMF, timeString=config_dt,ierr=err) - call mpas_get_timeInterval(timeStepESMF, dt=dtSim) + dtSim = forcingTimeIncrement + dtSimReverse = -dtSim ! use the forcing layer to get data call MPAS_forcing_get_forcing(& @@ -277,8 +285,18 @@ subroutine atmospheric_forcing(streamManager, domain, simulationClock)!{{{ forcingGroupHead, & ! forcingGroupHead "ocn_atmospheric_forcing", & ! forcingGroupName currentForcingTime) ! forcingTime - !call mpas_get_time(curr_time=currentForcingTime, dateTimeString=timeStamp, ierr=err) - !call mpas_log_write('Forcing time ' // trim(timeStamp)) + call mpas_get_time(curr_time=currentForcingTime, dateTimeString=timeStamp, ierr=err) + !call mpas_log_write('Forcing time for atmospheric forcing' // trim(timeStamp)) + + call mpas_advance_forcing_clock(forcingGroupHead, dtSimReverse) + + call MPAS_forcing_get_forcing_time(& + forcingGroupHead, & ! forcingGroupHead + "ocn_atmospheric_forcing", & ! forcingGroupName + currentForcingTime) ! forcingTime + + call mpas_get_time(curr_time=currentForcingTime, dateTimeString=timeStamp, ierr=err) + !call mpas_log_write('Forcing time reversed for atmospheric forcing' // trim(timeStamp)) ! perform post forcing block => domain % blocklist @@ -332,7 +350,10 @@ subroutine post_atmospheric_forcing(block)!{{{ windStressCoefficient, & rhoAir, & ramp, & - windStressCoefficientLimit + windStressCoefficientLimit, & + t + + real(kind=RKIND), pointer :: forcingTimeIncrement integer, pointer :: & nCells @@ -354,12 +375,14 @@ subroutine post_atmospheric_forcing(block)!{{{ call MPAS_pool_get_array(timeVaryingForcingPool, "atmosPressure", atmosPressure) call MPAS_pool_get_array(forcingPool, "atmosphericPressure", atmosphericPressure) + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) rhoAir = 1.225_RKIND windStressCoefficientLimit = 0.0035_RKIND if (daysSinceStartOfSim >= config_time_varying_atmospheric_forcing_ramp_delay) then - ramp = tanh((2.0_RKIND*(daysSinceStartOfSim-config_time_varying_atmospheric_forcing_ramp_delay)) & + t = (daysSinceStartOfSim*86400_RKIND + forcingTimeIncrement)/86400.0_RKIND + ramp = tanh((2.0_RKIND*(t-config_time_varying_atmospheric_forcing_ramp_delay)) & /config_time_varying_atmospheric_forcing_ramp) else ramp = 0.0_RKIND @@ -483,11 +506,13 @@ subroutine land_ice_forcing(streamManager, domain, simulationClock)!{{{ type (block_type), pointer :: block -! character(len=StrKIND), pointer :: config_dt + character(len=StrKIND) :: timeStamp type (MPAS_timeInterval_type) :: timeStepESMF - real (kind=RKIND) :: dtSim + real (kind=RKIND) :: dtSim, dtSimReverse + + real (kind=RKIND), pointer :: forcingTimeIncrement integer :: err @@ -512,14 +537,41 @@ subroutine land_ice_forcing(streamManager, domain, simulationClock)!{{{ integer :: & iCell + call mpas_pool_get_subpool(domain % blocklist % structs, 'forcing', forcingPool) + call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) + ! convert config_dt to real call mpas_set_timeInterval(timeStepESMF, timeString=config_dt,ierr=err) call mpas_get_timeInterval(timeStepESMF, dt=dtSim) + dtSim = forcingTimeIncrement + dtSimReverse = -dtSim + ! use the forcing layer to get data - call mpas_forcing_get_forcing(forcingGroupHead, "ocn_land_ice_forcing", streamManager, dtSim) + call MPAS_forcing_get_forcing(& + forcingGroupHead, & ! forcingGroupHead + "ocn_land_ice_forcing", & ! forcingGroupName + streamManager, & ! streamManager + dtSim) ! dt + + call MPAS_forcing_get_forcing_time(& + forcingGroupHead, & ! forcingGroupHead + "ocn_land_ice_forcing", & ! forcingGroupName + currentForcingTime) ! forcingTime + + call mpas_get_time(curr_time=currentForcingTime, dateTimeString=timeStamp, ierr=err) + call mpas_log_write('Forcing time for land ice forcing' // trim(timeStamp)) + + + call mpas_advance_forcing_clock(forcingGroupHead, dtSimReverse) + + call MPAS_forcing_get_forcing_time(& + forcingGroupHead, & ! forcingGroupHead + "ocn_land_ice_forcing", & ! forcingGroupName + currentForcingTime) ! forcingTime - call mpas_forcing_get_forcing_time(forcingGroupHead, "ocn_land_ice_forcing", currentForcingTime) + call mpas_get_time(curr_time=currentForcingTime, dateTimeString=timeStamp, ierr=err) + call mpas_log_write('Forcing time reversed for land ice forcing' // trim(timeStamp)) block => domain % blocklist do while (associated(block)) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F index 211538ac0a73..4eba6e2c0b45 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F @@ -80,6 +80,8 @@ module ocn_vel_tidal_potential tidalConstituentAstronomical, &! tidalConstituentNodalPhase ! + real (kind=RKIND), pointer :: forcingTimeIncrement + real (kind=RKIND), dimension(:,:), pointer :: & latitudeFunction @@ -273,8 +275,8 @@ subroutine ocn_compute_tidal_potential_forcing(err)!{{{ err = 0 if (tidalPotentialOff) return - ramp = tanh((2.0_RKIND*daysSinceStartOfSim)/tidalPotRamp) - t = daysSinceStartOfSim*86400.0_RKIND + t = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement + ramp = tanh((2.0_RKIND*t/86400.0_RKIND)/tidalPotRamp) do iCell = 1, nCellsAll tidalPotEta(iCell) = 0.0_RKIND @@ -416,6 +418,9 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{ call mpas_pool_get_array(forcingPool, & 'tidalPotentialLatitudeFunction', & latitudeFunction) + call mpas_pool_get_array(forcingPool, & + 'forcingTimeIncrement', & + forcingTimeIncrement) call mpas_set_time(refTime, & dateTimeString=config_tidal_potential_reference_time)