Skip to content

Commit

Permalink
Isotopes for Icepack (#307)
Browse files Browse the repository at this point in the history 
* initial isotope tracer mods

* calls to isotope routines

* correct isotope tracer indices

* update isotope interfaces to be backwards compatible
demonstrate prototype optional argument handling internal to the columnphysics
migrate to fortran use for n_iso and tr_iso, remove them from interfaces (to be done to other similar variables)
fix isotope hardcoded basal ice growth and isotope update logic
refactor optional closing argument to follow new approach (get ride of "extra" ridge_ice call)
fix isotope restart bug (read unit number)
add isotope restart test to base suite
add subnames to icepack_isotope.F90
rename public icepack_isotope variable frac to isotope_frac_method to improve clarity
add some coding standard documentation to developers guide under column physics section

* fix document rst

* fix comment, move declaration

* update acct on badger

* Add isotope documentation.

* Fix some typos

* updating documentation

* fixing format of Brady bib entry

* updating copyright date to 2020

* minor change to relaunch readthedocs check

Co-authored-by: apcraig <anthony.p.craig@gmail.com>
Co-authored-by: David Bailey <dbailey@ucar.edu>
  • Loading branch information
@eclare108213 @apcraig @dabail10
3 people authored Apr 6, 2020
1 parent 826dc88 commit 1ae0446
Show file tree
Hide file tree
Showing 35 changed files with 1,548 additions and 212 deletions.
Binary file modified LICENSE.pdf
View file
Binary file not shown.
106 changes: 90 additions & 16 deletions columnphysics/icepack_atmo.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@ module icepack_atmo
use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow
use icepack_parameters, only: atmbndy, calc_strair, formdrag
use icepack_parameters, only: highfreq, natmiter
use icepack_tracers, only: n_iso
use icepack_tracers, only: tr_iso
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted

Expand Down Expand Up @@ -61,6 +63,8 @@ subroutine atmo_boundary_layer (sfctype, &
lhcoef, shcoef, &
Cdn_atm, &
Cdn_atm_ratio_n, &
Qa_iso, Qref_iso, &
iso_flag, &
uvel, vvel, &
Uref )

Expand Down Expand Up @@ -104,6 +108,15 @@ subroutine atmo_boundary_layer (sfctype, &
shcoef , & ! transfer coefficient for sensible heat
lhcoef ! transfer coefficient for latent heat

logical (kind=log_kind), intent(in), optional :: &
iso_flag ! flag to trigger iso calculations

real (kind=dbl_kind), intent(in), optional, dimension(:) :: &
Qa_iso ! specific isotopic humidity (kg/kg)

real (kind=dbl_kind), intent(inout), optional, dimension(:) :: &
Qref_iso ! reference specific isotopic humidity (kg/kg)

real (kind=dbl_kind), intent(in) :: &
uvel , & ! x-direction ice speed (m/s)
vvel ! y-direction ice speed (m/s)
Expand All @@ -114,7 +127,7 @@ subroutine atmo_boundary_layer (sfctype, &
! local variables

integer (kind=int_kind) :: &
k ! iteration index
k,n ! iteration index

real (kind=dbl_kind) :: &
TsfK , & ! surface temperature in Kelvin (K)
Expand All @@ -136,6 +149,7 @@ subroutine atmo_boundary_layer (sfctype, &
ustar , & ! ustar (m/s)
tstar , & ! tstar
qstar , & ! qstar
ratio , & ! ratio
rdn , & ! sqrt of neutral exchange coefficient (momentum)
rhn , & ! sqrt of neutral exchange coefficient (heat)
ren , & ! sqrt of neutral exchange coefficient (water)
Expand All @@ -152,10 +166,18 @@ subroutine atmo_boundary_layer (sfctype, &
psixh ! stability function at zlvl (heat and water)

real (kind=dbl_kind), parameter :: &
zTrf = c2 ! reference height for air temp (m)
zTrf = c2 ! reference height for air temp (m)

logical (kind=log_kind) :: &
l_iso_flag ! local flag to trigger iso calculations

character(len=*),parameter :: subname='(atmo_boundary_layer)'

l_iso_flag = .false.
if (present(iso_flag)) then
l_iso_flag = iso_flag
endif

al2 = log(zref/zTrf)

!------------------------------------------------------------
Expand Down Expand Up @@ -356,6 +378,23 @@ subroutine atmo_boundary_layer (sfctype, &
Uref = vmag * rd / rdn
endif

if (l_iso_flag) then
if (present(Qref_iso) .and. present(Qa_iso)) then
Qref_iso(:) = c0
if (tr_iso) then
do n = 1, n_iso
ratio = c0
if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2)
Qref_iso(n) = Qa_iso(n) - ratio*delq*fac
enddo
endif
else
call icepack_warnings_add(subname//' l_iso_flag true but optional arrays missing')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
endif

end subroutine atmo_boundary_layer

!=======================================================================
Expand Down Expand Up @@ -468,7 +507,7 @@ end subroutine atmo_boundary_const
!=======================================================================

! Neutral drag coefficients for ocean and atmosphere also compute the
! intermediate necessary variables ridge height, distance, floe size
! intermediate necessary variables ridge height, distance, floe size
! based upon Tsamados et al. (2014), JPO, DOI: 10.1175/JPO-D-13-0215.1.
! Places where the code varies from the paper are commented.
!
Expand Down Expand Up @@ -574,7 +613,7 @@ subroutine neutral_drag_coeffs (apnd, hpnd, &
ctecaf, & ! constante
ctecwf, & ! constante
sca, & ! wind attenuation function
scw, & ! ocean attenuation function
scw, & ! ocean attenuation function
lp, & ! pond length (m)
ctecar, &
ctecwk, &
Expand Down Expand Up @@ -698,7 +737,7 @@ subroutine neutral_drag_coeffs (apnd, hpnd, &

!------------------------------------------------------------
! Skin drag (atmo)
!------------------------------------------------------------
!------------------------------------------------------------

Cdn_atm_skin = csa*(c1 - mrdg*tmp1/distrdg)
Cdn_atm_skin = max(min(Cdn_atm_skin,camax),c0)
Expand All @@ -725,7 +764,7 @@ subroutine neutral_drag_coeffs (apnd, hpnd, &

!------------------------------------------------------------
! Skin drag bottom ice (ocean)
!------------------------------------------------------------
!------------------------------------------------------------

Cdn_ocn_skin = csw * (c1 - mrdgo*tmp1/dkeel)
Cdn_ocn_skin = max(min(Cdn_ocn_skin,cwmax), c0)
Expand Down Expand Up @@ -800,7 +839,7 @@ end subroutine neutral_drag_coeffs
!autodocument_start icepack_atm_boundary
!

subroutine icepack_atm_boundary(sfctype, &
subroutine icepack_atm_boundary(sfctype, &
Tsf, potT, &
uatm, vatm, &
wind, zlvl, &
Expand All @@ -811,6 +850,7 @@ subroutine icepack_atm_boundary(sfctype, &
lhcoef, shcoef, &
Cdn_atm, &
Cdn_atm_ratio_n, &
Qa_iso, Qref_iso, &
uvel, vvel, &
Uref)

Expand Down Expand Up @@ -844,6 +884,12 @@ subroutine icepack_atm_boundary(sfctype, &
shcoef , & ! transfer coefficient for sensible heat
lhcoef ! transfer coefficient for latent heat

real (kind=dbl_kind), intent(in), optional, dimension(:) :: &
Qa_iso ! specific isotopic humidity (kg/kg)

real (kind=dbl_kind), intent(inout), optional, dimension(:) :: &
Qref_iso ! reference specific isotopic humidity (kg/kg)

real (kind=dbl_kind), optional, intent(in) :: &
uvel , & ! x-direction ice speed (m/s)
vvel ! y-direction ice speed (m/s)
Expand All @@ -856,18 +902,37 @@ subroutine icepack_atm_boundary(sfctype, &
! local variables

real (kind=dbl_kind) :: &
worku, workv, workr
l_uvel, l_vvel, l_Uref

real (kind=dbl_kind), dimension(:), allocatable :: &
l_Qa_iso, l_Qref_iso ! local copies of Qa_iso, Qref_iso

logical (kind=log_kind) :: &
iso_flag ! flag to turn on iso calcs in other subroutines

character(len=*),parameter :: subname='(icepack_atm_boundary)'

worku = c0
workv = c0
workr = c0
l_uvel = c0
l_vvel = c0
l_Uref = c0
if (present(uvel)) then
worku = uvel
l_uvel = uvel
endif
if (present(vvel)) then
workv = vvel
l_vvel = vvel
endif
if (present(Qa_iso) .and. present(Qref_iso)) then
iso_flag = .true.
allocate(l_Qa_iso(size(Qa_iso,dim=1)))
allocate(l_Qref_iso(size(Qref_iso,dim=1)))
l_Qa_iso = Qa_iso
l_Qref_iso = Qref_iso
else
iso_flag = .false.
allocate(l_Qa_iso(1))
allocate(l_Qref_iso(1))
l_Qa_iso = c0
l_Qref_iso = c0
endif

Cdn_atm_ratio_n = c1
Expand Down Expand Up @@ -896,15 +961,24 @@ subroutine icepack_atm_boundary(sfctype, &
lhcoef, shcoef, &
Cdn_atm, &
Cdn_atm_ratio_n, &
worku, workv, &
workr)
iso_flag = iso_flag, &
Qa_iso=l_Qa_iso, &
Qref_iso=l_Qref_iso, &
uvel=l_uvel, vvel=l_vvel, &
Uref=l_Uref)
if (icepack_warnings_aborted(subname)) return
endif ! atmbndy

if (present(Uref)) then
Uref = workr
Uref = l_Uref
endif

if (present(Qref_iso)) then
Qref_iso = l_Qref_iso
endif

deallocate(l_Qa_iso,l_Qref_iso)

end subroutine icepack_atm_boundary

!------------------------------------------------------------
Expand Down
29 changes: 28 additions & 1 deletion columnphysics/icepack_flux.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ module icepack_flux
use icepack_parameters, only: c1, emissivity
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
use icepack_tracers, only: tr_iso

implicit none
private
Expand Down Expand Up @@ -55,7 +56,10 @@ subroutine merge_fluxes (aicen, &
meltt, melts, &
meltb, &
congel, snoice, &
Uref, Urefn )
Uref, Urefn, &
Qref_iso, Qrefn_iso, &
fiso_ocn, fiso_ocnn, &
fiso_evap, fiso_evapn)

! single category fluxes
real (kind=dbl_kind), intent(in) :: &
Expand Down Expand Up @@ -119,6 +123,16 @@ subroutine merge_fluxes (aicen, &
real (kind=dbl_kind), optional, intent(inout):: &
Uref ! air speed reference level (m/s)

real (kind=dbl_kind), optional, dimension(:), intent(inout):: &
Qref_iso, & ! isotope air sp hum reference level (kg/kg)
fiso_ocn, & ! isotope fluxes to ocean (kg/m2/s)
fiso_evap ! isotope evaporation (kg/m2/s)

real (kind=dbl_kind), optional, dimension(:), intent(in):: &
Qrefn_iso, & ! isotope air sp hum reference level (kg/kg)
fiso_ocnn, & ! isotope fluxes to ocean (kg/m2/s)
fiso_evapn ! isotope evaporation (kg/m2/s)

character(len=*),parameter :: subname='(merge_fluxes)'

!-----------------------------------------------------------------
Expand Down Expand Up @@ -148,6 +162,19 @@ subroutine merge_fluxes (aicen, &
Tref = Tref + Trefn * aicen
Qref = Qref + Qrefn * aicen

! Isotopes
if (tr_iso) then
if (present(Qrefn_iso) .and. present(Qref_iso)) then
Qref_iso (:) = Qref_iso (:) + Qrefn_iso (:) * aicen
endif
if (present(fiso_ocnn) .and. present(fiso_ocn)) then
fiso_ocn (:) = fiso_ocn (:) + fiso_ocnn (:) * aicen
endif
if (present(fiso_evapn) .and. present(fiso_evap)) then
fiso_evap(:) = fiso_evap(:) + fiso_evapn(:) * aicen
endif
endif

! ocean fluxes
if (present(Urefn) .and. present(Uref)) then
Uref = Uref + Urefn * aicen
Expand Down
1 change: 1 addition & 0 deletions columnphysics/icepack_intfc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module icepack_intfc
use icepack_tracers, only: icepack_max_don => max_don
use icepack_tracers, only: icepack_max_fe => max_fe
use icepack_tracers, only: icepack_max_aero => max_aero
use icepack_tracers, only: icepack_max_iso => max_iso
use icepack_tracers, only: icepack_nmodal1 => nmodal1
use icepack_tracers, only: icepack_nmodal2 => nmodal2
use icepack_parameters, only: icepack_nspint => nspint
Expand Down
Loading

0 comments on commit 1ae0446

Please sign in to comment.