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.

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

Isotopes for Icepack #307

Merged
merged 17 commits into from
Apr 6, 2020
Merged
Show file tree
Hide file tree
Changes from 10 commits
Commits
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
106 changes: 90 additions & 16 deletions columnphysics/icepack_atmo.F90
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
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
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