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

Modifications required for running with a global slab ocean #40

Merged
merged 8 commits into from
Feb 28, 2024
Merged

Modifications required for running with a global slab ocean #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
95d55da
Permit running with a global slab ocean
spencerkclark Feb 6, 2024
f35704e
Add option to disable quasi-sea-ice
spencerkclark Feb 6, 2024
234566f
A couple further refinements
spencerkclark Feb 8, 2024
bfd2d6b
Initialize ocean to con_tice instead of previous timestep temperature
spencerkclark Feb 14, 2024
383ee3b
Harmonize formatting
spencerkclark Feb 14, 2024
48a24fe
Improve comment on added ocean_nml.maxlat parameter
spencerkclark Feb 28, 2024
47e783f
Use more generic word than 'slab'
spencerkclark Feb 28, 2024
992580c
typo: occuring -> occurring
spencerkclark Feb 28, 2024
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
12 changes: 12 additions & 0 deletions FV3GFS/FV3GFS_io.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3093,6 +3093,18 @@ subroutine register_diag_manager_controlled_diagnostics(Time, Sfcprop, IntDiag,
Diag_diag_manager_controlled(index)%data(nb)%var2 => Sfcprop(nb)%mld(:)
enddo

index = index + 1
Diag_diag_manager_controlled(index)%axes = 2
Diag_diag_manager_controlled(index)%name = 'prescribed_mixed_layer_depth'
Diag_diag_manager_controlled(index)%desc = 'prescribed ocean mixed layer depth'
Diag_diag_manager_controlled(index)%unit = 'm'
Diag_diag_manager_controlled(index)%mod_name = 'gfs_phys'
Diag_diag_manager_controlled(index)%coarse_graining_method = AREA_WEIGHTED
allocate (Diag_diag_manager_controlled(index)%data(nblks))
do nb = 1,nblks
Diag_diag_manager_controlled(index)%data(nb)%var2 => Sfcprop(nb)%mldclim(:)
enddo

index = index + 1
Diag_diag_manager_controlled(index)%axes = 2
Diag_diag_manager_controlled(index)%name = 'prescribed_qflux'
Expand Down
3 changes: 2 additions & 1 deletion GFS_layer/GFS_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,8 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, &
si = (Init_parm%ak + Init_parm%bk * p_ref - Init_parm%ak(Model%levr+1)) &
/ (p_ref - Init_parm%ak(Model%levr+1))
call rad_initialize (si, Model%levr, Model%ictm, Model%isol, &
Model%ico2, Model%iaer, Model%ialb, Model%iems, &
Model%ico2, Model%iaer, Model%ialb, &
Model%disable_radiation_quasi_sea_ice, Model%iems, &
Model%ntcw, Model%num_p2d, Model%num_p3d, &
Model%npdf3d, Model%ntoz, &
Model%iovr_sw, Model%iovr_lw, Model%isubc_sw, &
Expand Down
17 changes: 17 additions & 0 deletions GFS_layer/GFS_typedefs.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -516,6 +516,13 @@ module GFS_typedefs
integer :: ico2 !< prescribed global mean value (old opernl)
integer :: ialb !< use climatology alb, based on sfc type
!< 1 => use modis based alb

logical :: disable_radiation_quasi_sea_ice
!< flag to disable
!< radiation code treating ocean grid
!< cells with temperature below
!< freezing as sea ice

integer :: iems !< use fixed value of 1.0
integer :: iaer !< default aerosol effect in sw only
integer :: iovr_sw !< sw: max-random overlap clouds
Expand Down Expand Up @@ -2137,6 +2144,13 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
integer :: ico2 = 0 !< prescribed global mean value (old opernl)
integer :: ialb = 0 !< use climatology alb, based on sfc type
!< 1 => use modis based alb

logical :: disable_radiation_quasi_sea_ice = .false.
!< flag to disable
!< radiation code treating ocean grid
!< cells with temperature below
!< freezing as sea ice

integer :: iems = 0 !< use fixed value of 1.0
integer :: iaer = 1 !< default aerosol effect in sw only
integer :: iovr_sw = 1 !< sw: max-random overlap clouds
Expand Down Expand Up @@ -2452,6 +2466,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
cplflx, cplwav, lsidea, &
!--- radiation parameters
fhswr, fhlwr, levr, nfxr, aero_in, iflip, isol, ico2, ialb, &
disable_radiation_quasi_sea_ice, &
isot, iems, iaer, iovr_sw, iovr_lw, ictm, isubc_sw, &
isubc_lw, crick_proof, ccnorm, lwhtr, swhtr, nkld, &
fixed_date, fixed_solhr, fixed_sollat, daily_mean, sollat, &
Expand Down Expand Up @@ -2615,6 +2630,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%isol = isol
Model%ico2 = ico2
Model%ialb = ialb
Model%disable_radiation_quasi_sea_ice = disable_radiation_quasi_sea_ice
Model%iems = iems
Model%iaer = iaer
Model%iovr_sw = iovr_sw
Expand Down Expand Up @@ -3310,6 +3326,7 @@ subroutine control_print(Model)
print *, ' isol : ', Model%isol
print *, ' ico2 : ', Model%ico2
print *, ' ialb : ', Model%ialb
print *, ' disable_radiation_quasi_sea_ice: ', Model%disable_radiation_quasi_sea_ice
print *, ' iems : ', Model%iems
print *, ' iaer : ', Model%iaer
print *, ' iovr_sw : ', Model%iovr_sw
Expand Down
5 changes: 4 additions & 1 deletion gsmphys/physparam.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,10 @@ module physparam !
!> surface albedo scheme control flag
!!\n =0:vegetation type based climatological albedo scheme
!!\n =1:seasonal albedo derived from MODIS measurements
integer, save :: ialbflg = 0
integer, save :: ialbflg = 0
!> flag to disable radiation code treating ocean grid cells with
!!\n temperature below freezing as sea ice
logical, save :: ldisable_radiation_quasi_sea_ice = .false.
!> surface emissivity scheme control flag
!!\n =0:black-body surface emissivity(=1.0)
!!\n =1:vegetation type based climatology emissivity(<1.0)
Expand Down
15 changes: 10 additions & 5 deletions gsmphys/rad_initialize.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@
subroutine rad_initialize &
!...................................
! --- inputs:
& ( si,levr,ictm,isol,ico2,iaer,ialb,iems,ntcw, num_p2d, &
& ( si,levr,ictm,isol,ico2,iaer,ialb, &
& disable_radiation_quasi_sea_ice,iems,ntcw, num_p2d, &
& num_p3d,npdf3d,ntoz,iovr_sw,iovr_lw,isubc_sw,isubc_lw, &
& crick_proof,ccnorm,norad_precip, &
& idate,iflip,me )
Expand Down Expand Up @@ -107,7 +108,7 @@ subroutine rad_initialize &
& iaermdl, laswflg, lalwflg, lavoflg, icldflg, icmphys,&
& iovrsw , iovrlw , lcrick , lcnorm , lnoprec, &
& ialbflg, iemsflg, isubcsw, isubclw, ivflip , ipsd0, &
& kind_phys
& kind_phys, ldisable_radiation_quasi_sea_ice

use module_radiation_driver, only : radinit
!
Expand All @@ -120,8 +121,8 @@ subroutine rad_initialize &

real (kind=kind_phys), intent(in) :: si(levr+1)

logical, intent(in) :: crick_proof, ccnorm, norad_precip

logical, intent(in) :: crick_proof, ccnorm, norad_precip, &
& disable_radiation_quasi_sea_ice
! --- output: ( none )

! --- local:
Expand Down Expand Up @@ -184,6 +185,7 @@ subroutine rad_initialize &
isubclw = isubc_lw ! sub-column cloud approx flag in lw radiation

ialbflg= ialb ! surface albedo control flag
ldisable_radiation_quasi_sea_ice = disable_radiation_quasi_sea_ice ! flag to disable treating below freezing ocean grid cells as sea ice
iemsflg= iems ! surface emissivity control flag

ivflip = iflip ! vertical index direction control flag
Expand All @@ -198,7 +200,10 @@ subroutine rad_initialize &
print *,' In rad_initialize, before calling radinit'
print *,' si =',si
print *,' levr=',levr,' ictm=',ictm,' isol=',isol,' ico2=',ico2,&
& ' iaer=',iaer,' ialb=',ialb,' iems=',iems,' ntcw=',ntcw
& ' iaer=',iaer,' ialb=',ialb, &
& ' disable_radiation_quasi_sea_ice=', &
& disable_radiation_quasi_sea_ice, &
& ' iems=',iems,' ntcw=',ntcw
print *,' np3d=',num_p3d,' ntoz=',ntoz,' iovr_sw=',iovr_sw, &
& ' iovr_lw=',iovr_lw,' isubc_sw=',isubc_sw, &
& ' isubc_lw=',isubc_lw,' iflip=',iflip,' me=',me
Expand Down
59 changes: 48 additions & 11 deletions gsmphys/radiation_surface.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,8 @@ module module_radiation_surface !
!........................................!
!
use physparam, only : ialbflg, iemsflg, semis_file, &
& kind_phys
& kind_phys, &
& ldisable_radiation_quasi_sea_ice
use physcons, only : con_t0c, con_ttp, con_pi, con_tice
use module_iounitdef, only : NIRADSF
!
Expand Down Expand Up @@ -207,6 +208,21 @@ subroutine sfc_init &
stop
endif ! end if_ialbflg_block

!! \n physparam::ldisable_radiation_quasi_sea_ice
!! - = .false.: use a sea-ice-like albedo and emissivity for below
!! freezing ocean grid cells.
!! - = .true.: treat all ocean grid cells as if they were
!! above freezing when determing the albedo and emissivity
if ( ldisable_radiation_quasi_sea_ice ) then
if ( me == 0 ) then
print *, '- Disabling radiation quasi-sea-ice'
endif
else
if ( me == 0 ) then
print *, '- Enabling radiation quasi-sea-ice'
endif
endif

!> - Initialization of surface emissivity section
!! \n physparam::iemsflg
!! - = 0: fixed SFC emissivity at 1.0
Expand Down Expand Up @@ -424,7 +440,13 @@ subroutine setalb &
argh = min(0.50, max(.025, 0.01*zorlf(i)))
hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) )
fsno0 = asnow / (argh + asnow) * hrgh
if (nint(slmsk(i))==0 .and. tsknf(i)>con_tice) fsno0 = f_zero

if (nint(slmsk(i))==0 .and.
& (tsknf(i)>con_tice .or.
& ldisable_radiation_quasi_sea_ice)) then
fsno0 = f_zero
endif

fsno1 = f_one - fsno0
flnd0 = min(f_one, facsf(i)+facwf(i))
fsea0 = max(f_zero, f_one-flnd0)
Expand All @@ -434,7 +456,8 @@ subroutine setalb &

!> - Calculate diffused sea surface albedo

if (tsknf(i) >= 271.5) then
if (tsknf(i) >= 271.5 .or.
& ldisable_radiation_quasi_sea_ice) then
asevd = 0.06
asend = 0.06
elseif (tsknf(i) < 271.1) then
Expand Down Expand Up @@ -486,7 +509,8 @@ subroutine setalb &
rfcs = 2.14 / (f_one + 1.48*coszf(i))
rfcw = rfcs

if (tsknf(i) >= con_t0c) then
if (tsknf(i) >= con_t0c .or.
& ldisable_radiation_quasi_sea_ice) then
asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) &
& + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
& * (coszf(i)-f_one))
Expand Down Expand Up @@ -525,7 +549,11 @@ subroutine setalb &

fsno0 = sncovr(i)

if (nint(slmsk(i))==0 .and. tsknf(i)>con_tice) fsno0 = f_zero
if (nint(slmsk(i))==0 .and.
& (tsknf(i)>con_tice .or.
& ldisable_radiation_quasi_sea_ice)) then
fsno0 = f_zero
endif

if (nint(slmsk(i)) == 2) then
asnow = 0.02*snowf(i)
Expand All @@ -543,7 +571,8 @@ subroutine setalb &

!> - Calculate diffused sea surface albedo.

if (tsknf(i) >= 271.5) then
if (tsknf(i) >= 271.5 .or.
& ldisable_radiation_quasi_sea_ice) then
asevd = 0.06
asend = 0.06
elseif (tsknf(i) < 271.1) then
Expand Down Expand Up @@ -602,7 +631,8 @@ subroutine setalb &
! & - 2.02*coszf(i)*coszf(i)*coszf(i)
rfcs = 1.775/(1.0+1.55*coszf(i))

if (tsknf(i) >= con_t0c) then
if (tsknf(i) >= con_t0c .or.
& ldisable_radiation_quasi_sea_ice) then
asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) &
& + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
& * (coszf(i)-f_one))
Expand Down Expand Up @@ -652,7 +682,8 @@ subroutine setalb &

!> - Calculate diffused sea surface albedo.

if (tsknf(i) >= 271.5) then
if (tsknf(i) >= 271.5 .or.
& ldisable_radiation_quasi_sea_ice) then
asevd = 0.06
asend = 0.06
elseif (tsknf(i) < 271.1) then
Expand Down Expand Up @@ -708,7 +739,8 @@ subroutine setalb &
! & - 2.02*coszf(i)*coszf(i)*coszf(i)
rfcs = 1.775/(1.0+1.55*coszf(i))

if (tsknf(i) >= con_t0c) then
if (tsknf(i) >= con_t0c .or.
& ldisable_radiation_quasi_sea_ice) then
asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) &
& + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
& * (coszf(i)-f_one))
Expand Down Expand Up @@ -930,8 +962,13 @@ subroutine setemis &
argh = min(0.50, max(.025, 0.01*zorlf(i)))
hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) )
fsno0 = asnow / (argh + asnow) * hrgh
if (nint(slmsk(i)) == 0 .and. tsknf(i) > 271.2) &
& fsno0=f_zero

if (nint(slmsk(i)) == 0 .and.
& (tsknf(i) > 271.2 .or.
& ldisable_radiation_quasi_sea_ice)) then
fsno0=f_zero
endif

fsno1 = f_one - fsno0
sfcemis(i) = sfcemis(i)*fsno1 + emsref(8)*fsno0
endif
Expand Down
42 changes: 35 additions & 7 deletions gsmphys/som_mlm.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
module module_ocean !
!........................................!
!
use physcons, only : con_tice
use physparam, only : kind_phys
use GFS_typedefs, only : GFS_control_type, GFS_grid_type
! use constants_mod, only : omega, grav
Expand All @@ -35,9 +36,8 @@ module module_ocean !

public ocean_init, update_ocean
!
real (kind=kind_phys) :: maxlat, width_buffer, minmld, &
real (kind=kind_phys) :: width_buffer, minmld, &
cpwater, rhowater, omega, grav
parameter(maxlat = 60.) ! determine the maximum latitude band for SOM/MLM
parameter(minmld = 10.) ! minimum mixed layer depth
parameter(width_buffer = 15.) ! the width of a buffer band where SST is determined by both SOM/MLM
! and climatology (or climatology plus initial anomaly)
Expand Down Expand Up @@ -66,8 +66,8 @@ module module_ocean !
real(kind=kind_phys) :: eps_day = 10. ! damping time scale of ocean current (days)
real(kind=kind_phys) :: sst_restore_tscale = 3. ! restoring time scale for sst (day)
real(kind=kind_phys) :: mld_restore_tscale = 1. ! restoring time scale for mld (day)
real(kind=kind_phys) :: start_lat = -30. ! latitude starting from? Note that this value should not be smaller than -60.
real(kind=kind_phys) :: end_lat = 30. ! latitude ending with? Note that this value should not be bigger than 60.
real(kind=kind_phys) :: start_lat = -30. ! latitude starting from? Note that this value should not be smaller than -maxlat.
real(kind=kind_phys) :: end_lat = 30. ! latitude ending with? Note that this value should not be bigger than maxlat.
real(kind=kind_phys) :: tday1 = 3. !
real(kind=kind_phys) :: tday2 = 10. !
real(kind=kind_phys) :: sst_restore_tscale1= 3. ! restoring time scale for sst during the period from 1 to tday1
Expand All @@ -78,13 +78,20 @@ module module_ocean !
! climatological SST plus initial anomaly
logical :: use_tvar_restore_sst = .false.! using time varying restoring time scale for sst
logical :: use_tvar_restore_mld = .false.! using time varying restoring time scale for mld
real(kind=kind_phys) :: maxlat = 60. ! maximum latitudinal extent of the SOM/MLM. Generally set to <= 60, though
! it can be useful to set to > 60 if the desired start_lat and end_lat are
! poleward of 60 degrees. If set to > 60, it is recommended to also set
! gfs_physics_nml.disable_radiation_quasi_sea_ice to .true. to prevent
! an unphysical quasi-ice-albedo feedback from occurring. Set to 90 along
! with ocean_nml.start_lat = -90 and ocean_nml.end_lat = 90 to enable
! running with a global interactive ocean.

namelist /ocean_nml/ &
ocean_option, mld_option, mld_obs_ratio, stress_ratio, restore_method, &
use_old_mlm, use_rain_flux, use_qflux, do_mld_restore, const_mld, Gam, &
eps_day, sst_restore_tscale, mld_restore_tscale, start_lat, end_lat, &
tday1, tday2, sst_restore_tscale1, sst_restore_tscale2, mld_restore_tscale1, &
mld_restore_tscale2, use_tvar_restore_sst, use_tvar_restore_mld
mld_restore_tscale2, use_tvar_restore_sst, use_tvar_restore_mld, maxlat

! =================
contains
Expand Down Expand Up @@ -149,11 +156,11 @@ subroutine ocean_init &
#endif

if (start_lat < -maxlat) then
write(*,*) 'start_lat should not be smaller than -60.'
write(*,*) 'start_lat should not be smaller than', -maxlat
call abort
endif
if (end_lat > maxlat) then
write(*,*) 'end_lat should not be larger than 60.'
write(*,*) 'end_lat should not be larger than', maxlat
call abort
endif

Expand Down Expand Up @@ -295,6 +302,27 @@ subroutine update_ocean &
enddo
endif
endif

! Reset the slab ocean temperature or mixed layer ocean properties in
! any sea ice grid cells. This is required in the (rare) instances
! when one is running with the SOM/MLM with a latitudinal extent such
! that grid cells transition from ocean to sea ice and back.
where (islmsk .eq. 2)
ts_som = con_tice
endwhere

if (ocean_option == "MLM") then
where(islmsk .eq. 2)
tml = con_tice
tml0 = con_tice
mld = mldclim
mld0 = mld
huml = 0.
hvml = 0.
tmoml = con_tice - 5.
tmoml0 = tmoml
endwhere
endif
!
if (use_rain_flux) then
do i=1,im
Expand Down