Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for disabling thermodynamics #389

Closed
wants to merge 2 commits into from
Closed

Add support for disabling thermodynamics #389

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5d5bc79
icepack_therm_vertical: let 'ktherm=-1' only disable thermodynamics
phil-blain Sep 28, 2021
2419c72
icedrv_RunMod: add support for disabling thermodynamics
phil-blain Apr 7, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
230 changes: 118 additions & 112 deletions columnphysics/icepack_therm_vertical.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2586,44 +2586,47 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
l_smliq = c0
if (present(smliqn)) l_smliq = smliqn

if (ktherm >= 0) then

!-----------------------------------------------------------------
! Initialize rate of snow loss to leads
!-----------------------------------------------------------------

l_fsloss = fsnow * (c1 - aice)
l_fsloss = fsnow * (c1 - aice)

!-----------------------------------------------------------------
! snow redistribution using snwlvlfac: precip factor
!-----------------------------------------------------------------

if (trim(snwredist) == 'bulk') then
worka = c0
if (aice > puny) then
do n = 1, ncat
worka = worka + alvl(n)*aicen(n)
enddo
worka = worka * (snwlvlfac/(c1+snwlvlfac)) / aice
endif
l_fsloss = l_fsloss + fsnow* worka
fsnow = fsnow*(c1-worka)
endif ! snwredist
if (trim(snwredist) == 'bulk') then
worka = c0
if (aice > puny) then
do n = 1, ncat
worka = worka + alvl(n)*aicen(n)
enddo
worka = worka * (snwlvlfac/(c1+snwlvlfac)) / aice
endif
l_fsloss = l_fsloss + fsnow* worka
fsnow = fsnow*(c1-worka)
endif ! snwredist

!-----------------------------------------------------------------
! solid and liquid components of snow mass
!-----------------------------------------------------------------

massicen(:,:) = c0
massliqn(:,:) = c0
if (snwgrain) then
rnslyr = c1 / real(nslyr, dbl_kind)
do n = 1, ncat
do k = 1, nslyr
massicen(k,n) = l_smice(k,n) * vsnon(n) * rnslyr ! kg/m^2
massliqn(k,n) = l_smliq(k,n) * vsnon(n) * rnslyr
massicen(:,:) = c0
massliqn(:,:) = c0
if (snwgrain) then
rnslyr = c1 / real(nslyr, dbl_kind)
do n = 1, ncat
do k = 1, nslyr
massicen(k,n) = l_smice(k,n) * vsnon(n) * rnslyr ! kg/m^2
massliqn(k,n) = l_smliq(k,n) * vsnon(n) * rnslyr
enddo
enddo
enddo
endif
endif

endif ! ktherm >= 0

!-----------------------------------------------------------------
! Update the neutral drag coefficients to account for form drag
Expand Down Expand Up @@ -2655,6 +2658,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
! Compute lateral and bottom heat fluxes.
!-----------------------------------------------------------------

if (ktherm >= 0) &
call frzmlt_bottom_lateral (dt, ncat, &
nilyr, nslyr, &
aice, frzmlt, &
Expand Down Expand Up @@ -2753,109 +2757,111 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
! Vertical thermodynamics: Heat conduction, growth and melting.
!-----------------------------------------------------------------

if (.not.(calc_Tsfc)) then
if (ktherm >= 0) then
if (.not.(calc_Tsfc)) then

! If not calculating surface temperature and fluxes, set
! surface fluxes (flatn, fsurfn, and fcondtopn) to be used
! in thickness_changes

! hadgem routine sets fluxes to default values in ice-only mode
call set_sfcflux(aicen (n), &
flatn_f (n), fsensn_f (n), &
fcondtopn_f(n), &
fsurfn_f (n), &
flatn (n), fsensn (n), &
fsurfn (n), &
fcondtopn (n))
if (icepack_warnings_aborted(subname)) return
endif
! If not calculating surface temperature and fluxes, set
! surface fluxes (flatn, fsurfn, and fcondtopn) to be used
! in thickness_changes
! hadgem routine sets fluxes to default values in ice-only mode
call set_sfcflux(aicen (n), &
flatn_f (n), fsensn_f (n), &
fcondtopn_f(n), &
fsurfn_f (n), &
flatn (n), fsensn (n), &
fsurfn (n), &
fcondtopn (n))
if (icepack_warnings_aborted(subname)) return
endif

call thermo_vertical(nilyr=nilyr, nslyr=nslyr, &
dt=dt, aicen=aicen (n), &
vicen=vicen (n), vsnon=vsnon (n), &
Tsf=Tsfc (n), zSin=zSin (:,n), &
zqin=zqin (:,n), zqsn=zqsn (:,n), &
apond=apnd (n), hpond=hpnd (n), &
flw=flw, potT=potT, &
Qa=Qa, rhoa=rhoa, &
fsnow=fsnow, fpond=fpond, &
fbot=fbot, Tbot=Tbot, &
Tsnice=Tsnice, sss=sss, &
rsnw=l_rsnw (:,n), &
lhcoef=lhcoef, shcoef=shcoef, &
fswsfc=fswsfcn (n), fswint=fswintn (n), &
Sswabs=Sswabsn(:,n), Iswabs=Iswabsn(:,n), &
fsurfn=fsurfn (n), fcondtopn=fcondtopn(n), &
fcondbotn=fcondbotn(n), &
fsensn=fsensn (n), flatn=flatn (n), &
flwoutn=flwoutn, evapn=evapn, &
evapsn=evapsn, evapin=evapin, &
freshn=freshn, fsaltn=fsaltn, &
fhocnn=fhocnn, frain=frain, &
meltt=melttn (n), melts=meltsn (n), &
meltb=meltbn (n), meltsliq=l_meltsliqn(n),&
smice=l_smice (:,n), massice=massicen(:,n), &
smliq=l_smliq (:,n), massliq=massliqn(:,n), &
congel=congeln (n), snoice=snoicen (n), &
mlt_onset=mlt_onset, frz_onset=frz_onset, &
yday=yday, dsnow=dsnown (n), &
prescribed_ice=prescribed_ice)

if (icepack_warnings_aborted(subname)) then
write(warnstr,*) subname, ' ice: Vertical thermo error, cat ', n
call icepack_warnings_add(warnstr)
return
endif
call thermo_vertical(nilyr=nilyr, nslyr=nslyr, &
dt=dt, aicen=aicen (n), &
vicen=vicen (n), vsnon=vsnon (n), &
Tsf=Tsfc (n), zSin=zSin (:,n), &
zqin=zqin (:,n), zqsn=zqsn (:,n), &
apond=apnd (n), hpond=hpnd (n), &
flw=flw, potT=potT, &
Qa=Qa, rhoa=rhoa, &
fsnow=fsnow, fpond=fpond, &
fbot=fbot, Tbot=Tbot, &
Tsnice=Tsnice, sss=sss, &
rsnw=l_rsnw (:,n), &
lhcoef=lhcoef, shcoef=shcoef, &
fswsfc=fswsfcn (n), fswint=fswintn (n), &
Sswabs=Sswabsn(:,n), Iswabs=Iswabsn(:,n), &
fsurfn=fsurfn (n), fcondtopn=fcondtopn(n), &
fcondbotn=fcondbotn(n), &
fsensn=fsensn (n), flatn=flatn (n), &
flwoutn=flwoutn, evapn=evapn, &
evapsn=evapsn, evapin=evapin, &
freshn=freshn, fsaltn=fsaltn, &
fhocnn=fhocnn, frain=frain, &
meltt=melttn (n), melts=meltsn (n), &
meltb=meltbn (n), meltsliq=l_meltsliqn(n),&
smice=l_smice (:,n), massice=massicen(:,n), &
smliq=l_smliq (:,n), massliq=massliqn(:,n), &
congel=congeln (n), snoice=snoicen (n), &
mlt_onset=mlt_onset, frz_onset=frz_onset, &
yday=yday, dsnow=dsnown (n), &
prescribed_ice=prescribed_ice)

if (icepack_warnings_aborted(subname)) then
write(warnstr,*) subname, ' ice: Vertical thermo error, cat ', n
call icepack_warnings_add(warnstr)
return
endif

!-----------------------------------------------------------------
! Total absorbed shortwave radiation
!-----------------------------------------------------------------

fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)
fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)

!-----------------------------------------------------------------
! Aerosol update
!-----------------------------------------------------------------

if (tr_aero) then
call update_aerosol (dt, &
nilyr, nslyr, n_aero, &
melttn (n), meltsn (n), &
meltbn (n), congeln (n), &
snoicen (n), fsnow, &
aerosno(:,:,n), aeroice(:,:,n), &
aicen_init (n), vicen_init (n), &
vsnon_init (n), &
vicen (n), vsnon (n), &
aicen (n), &
faero_atm , faero_ocn)
if (icepack_warnings_aborted(subname)) return
endif
if (tr_aero) then
call update_aerosol (dt, &
nilyr, nslyr, n_aero, &
melttn (n), meltsn (n), &
meltbn (n), congeln (n), &
snoicen (n), fsnow, &
aerosno(:,:,n), aeroice(:,:,n), &
aicen_init (n), vicen_init (n), &
vsnon_init (n), &
vicen (n), vsnon (n), &
aicen (n), &
faero_atm , faero_ocn)
if (icepack_warnings_aborted(subname)) return
endif

if (tr_iso) then
call update_isotope (dt = dt, &
nilyr = nilyr, nslyr = nslyr, &
meltt = melttn(n),melts = meltsn(n), &
meltb = meltbn(n),congel=congeln(n), &
snoice=snoicen(n),evap=evapn, &
fsnow=fsnow, Tsfc=Tsfc(n), &
Qref_iso=Qrefn_iso(:), &
isosno=l_isosno(:,n),isoice=l_isoice(:,n), &
aice_old=aicen_init(n),vice_old=vicen_init(n), &
vsno_old=vsnon_init(n), &
vicen=vicen(n),vsnon=vsnon(n), &
aicen=aicen(n), &
fiso_atm=l_fiso_atm(:), &
fiso_evapn=fiso_evapn(:), &
fiso_ocnn=fiso_ocnn(:), &
HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, &
H2_18O_ocn=l_H2_18O_ocn)
if (icepack_warnings_aborted(subname)) return
endif
if (tr_iso) then
call update_isotope (dt = dt, &
nilyr = nilyr, nslyr = nslyr, &
meltt = melttn(n),melts = meltsn(n), &
meltb = meltbn(n),congel=congeln(n), &
snoice=snoicen(n),evap=evapn, &
fsnow=fsnow, Tsfc=Tsfc(n), &
Qref_iso=Qrefn_iso(:), &
isosno=l_isosno(:,n),isoice=l_isoice(:,n), &
aice_old=aicen_init(n),vice_old=vicen_init(n), &
vsno_old=vsnon_init(n), &
vicen=vicen(n),vsnon=vsnon(n), &
aicen=aicen(n), &
fiso_atm=l_fiso_atm(:), &
fiso_evapn=fiso_evapn(:), &
fiso_ocnn=fiso_ocnn(:), &
HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, &
H2_18O_ocn=l_H2_18O_ocn)
if (icepack_warnings_aborted(subname)) return
endif

endif ! ktherm >= 0
endif ! aicen_init

if (snwgrain .and. use_smliq_pnd) then
if (snwgrain .and. use_smliq_pnd .and. ktherm >= 0) then
call drain_snow (nslyr = nslyr, &
vsnon = vsnon(n), &
aicen = aicen(n), &
Expand All @@ -2873,7 +2879,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
!-----------------------------------------------------------------

!call ice_timer_start(timer_ponds)
if (tr_pond) then
if (tr_pond .and. ktherm >= 0) then

if (tr_pond_cesm) then
rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
Expand Down Expand Up @@ -3012,7 +3018,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
! reload snow mass tracers
!-----------------------------------------------------------------

if (snwgrain) then
if (snwgrain .and. ktherm >= 0) then
do n = 1, ncat
if (vsnon(n) > puny) then
do k = 1, nslyr
Expand Down Expand Up @@ -3074,7 +3080,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
! Calculate ponds from the topographic scheme
!-----------------------------------------------------------------
!call ice_timer_start(timer_ponds)
if (tr_pond_topo) then
if (tr_pond_topo .and. ktherm >= 0) then
call compute_ponds_topo(dt, ncat, nilyr, &
ktherm, heat_capacity, &
aice, aicen, &
Expand Down
12 changes: 8 additions & 4 deletions configuration/driver/icedrv_RunMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,7 @@ subroutine ice_step
use icedrv_flux, only: init_history_therm, init_history_bgc, &
daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd, init_history_dyn
use icedrv_history, only: history_cdf, history_write
use icepack_parameters, only: ktherm
use icedrv_restart, only: dumpfile, final_restart
use icedrv_restart_bgc, only: write_restart_bgc
use icedrv_step, only: prep_radiation, step_therm1, step_therm2, &
Expand Down Expand Up @@ -148,7 +149,7 @@ subroutine ice_step
! Scale radiation fields
!-----------------------------------------------------------------

if (calc_Tsfc) call prep_radiation ()
if (calc_Tsfc .and. ktherm >= 0) call prep_radiation ()

! call icedrv_diagnostics_debug ('post prep_radiation')

Expand All @@ -157,8 +158,10 @@ subroutine ice_step
!-----------------------------------------------------------------

call step_therm1 (dt) ! vertical thermodynamics
call biogeochemistry (dt) ! biogeochemistry
call step_therm2 (dt) ! ice thickness distribution thermo
if (ktherm >= 0) then
call biogeochemistry (dt) ! biogeochemistry
call step_therm2 (dt) ! ice thickness distribution thermo
endif

! clean up, update tendency diagnostics
offset = dt
Expand Down Expand Up @@ -204,7 +207,8 @@ subroutine ice_step
! albedo, shortwave radiation
!-----------------------------------------------------------------

call step_radiation (dt)
if (ktherm >= 0) &
call step_radiation (dt)

!-----------------------------------------------------------------
! get ready for coupling and the next time step
Expand Down