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

Refactor optional argument implementation for isotopes, snwgrain, therm1 and therm2 #423

Merged
merged 4 commits into from
Jan 24, 2023
Merged

Refactor optional argument implementation for isotopes, snwgrain, therm1 and therm2 #423

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fb45b9e
Refactor isotope implementation with respect to optional arguments
apcraig Jan 11, 2023
7b2ad22
- Continue update of isotope implementation in Icepack including HDO_…
apcraig Jan 12, 2023
0c57334
- Update fswtrhu_[vi]d[fr] optional argument implementation
apcraig Jan 14, 2023
9e28d92
Create icepack_chkoptargflag function for optional argument checking
apcraig Jan 23, 2023
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
109 changes: 44 additions & 65 deletions columnphysics/icepack_atmo.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,13 +13,14 @@
module icepack_atmo

use icepack_kinds
use icepack_parameters, only: c0, c1, c2, c4, c5, c8, c10
use icepack_parameters, only: c16, c20, p001, p01, p2, p4, p5, p75, puny
use icepack_parameters, only: senscoef, latncoef
use icepack_parameters, only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
use icepack_parameters, only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow
use icepack_parameters, only: c0, c1, c2, c4, c5, c8, c10
use icepack_parameters, only: c16, c20, p001, p01, p2, p4, p5, p75, puny
use icepack_parameters, only: senscoef, latncoef
use icepack_parameters, only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
use icepack_parameters, only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow
use icepack_parameters, only: atmbndy, calc_strair, formdrag
use icepack_parameters, only: icepack_chkoptargflag
use icepack_tracers, only: n_iso
use icepack_tracers, only: tr_iso
use icepack_warnings, only: warnstr, icepack_warnings_add
Expand Down Expand Up @@ -61,7 +62,6 @@ subroutine atmo_boundary_layer (sfctype, &
Cdn_atm, &
Cdn_atm_ratio_n, &
Qa_iso, Qref_iso, &
iso_flag, &
uvel, vvel, &
Uref, zlvs )

Expand Down Expand Up @@ -103,13 +103,10 @@ 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(:) :: &
real (kind=dbl_kind), intent(in), dimension(:), optional :: &
Qa_iso ! specific isotopic humidity (kg/kg)

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

real (kind=dbl_kind), intent(in) :: &
Expand Down Expand Up @@ -167,16 +164,8 @@ subroutine atmo_boundary_layer (sfctype, &
real (kind=dbl_kind), parameter :: &
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 @@ -389,21 +378,13 @@ subroutine atmo_boundary_layer (sfctype, &
Uref = vmag * rd / rdn
endif

if (l_iso_flag) then
if (present(Qref_iso) .and. present(Qa_iso)) then
if (tr_iso .and. sfctype(1:3)=='ice') 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
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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Farther down in the comments for this module, you say

You can definitely get into trouble if you pass the optional argument down the calling tree and use it if it wasn't passed.

So are Qa_iso and Qref_iso safe to use here without checking their 'present' status because the call to this subroutine isn't the first call in which those arguments are declared optional? It appears that you are checking it in what would be the first call into icepack, icepack_atm_boundary. Is there some easy way to track when an optional argument needs to be checked 'if present' and when it doesn't? It's probably not as simple as just checking in icepack's interface routines.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's right, we have already checked the arguments in the highest level Icepack interface. So they are safe to use here without checking. That's our current design. As I mentioned the other day, there are pros and cons to doing it this way. A pro is that we check all the optional arguments early leveraging parameter flags and then we should be OK after that without having to check over and over. A negative is that you have to keep track of what's going on in different layers of the code.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, thank you.

enddo
endif

end subroutine atmo_boundary_layer
Expand Down Expand Up @@ -895,18 +876,18 @@ 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(:) :: &
real (kind=dbl_kind), intent(in), dimension(:), optional :: &
Qa_iso ! specific isotopic humidity (kg/kg)

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

real (kind=dbl_kind), optional, intent(in) :: &
real (kind=dbl_kind), intent(in), optional :: &
uvel , & ! x-direction ice speed (m/s)
vvel , & ! y-direction ice speed (m/s)
zlvs ! atm level height for scalars (if different than zlvl) (m)

real (kind=dbl_kind), optional, intent(out) :: &
real (kind=dbl_kind), intent(out), optional :: &
Uref ! reference height wind speed (m/s)

!autodocument_end
Expand All @@ -916,14 +897,28 @@ subroutine icepack_atm_boundary(sfctype, &
real (kind=dbl_kind) :: &
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
logical (kind=log_kind), save :: &
first_call_ice = .true. ! first call flag

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

!------------------------------------------------------------
! Check optional arguments
! Need separate first_call flags for 'ice' and 'ocn' sfctype
!------------------------------------------------------------

if (sfctype == 'ice') then
if (icepack_chkoptargflag(first_call_ice)) then
if (tr_iso) then
if (.not.(present(Qa_iso).and.present(Qref_iso))) then
call icepack_warnings_add(subname//' error in fiso_ocn argument, tr_iso=T')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
endif
endif
endif

l_uvel = c0
l_vvel = c0
l_Uref = c0
Expand All @@ -933,19 +928,6 @@ subroutine icepack_atm_boundary(sfctype, &
if (present(vvel)) then
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 All @@ -972,24 +954,21 @@ subroutine icepack_atm_boundary(sfctype, &
lhcoef, shcoef, &
Cdn_atm, &
Cdn_atm_ratio_n, &
iso_flag = iso_flag, &
Qa_iso=l_Qa_iso, &
Qref_iso=l_Qref_iso, &
uvel=l_uvel, vvel=l_vvel, &
Uref=l_Uref, zlvs=zlvs)
Qa_iso=Qa_iso, &
Qref_iso=Qref_iso, &
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens if present(Qref_iso)=F? Is this robust across all compilers/platforms? You probably tested that previously...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is robust. The optional attribute and "present" state is passed down the calling tree. As long as optional is set in the calling routine for this field if want to check it's "present" state (determined by the original caller) and/or the optional argument is not used unless it exists, all should be fine. You can definitely get into trouble if you pass the optional argument down the calling tree and use it if it wasn't passed. Hopefully, we'll avoid that with the logic we have setup.

uvel=l_uvel, vvel=l_vvel,&
Uref=l_Uref, zlvs=zlvs )
if (icepack_warnings_aborted(subname)) return
endif ! atmbndy

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

if (present(Qref_iso)) then
Qref_iso = l_Qref_iso
if (sfctype == 'ice') then
first_call_ice = .false.
endif

deallocate(l_Qa_iso,l_Qref_iso)

end subroutine icepack_atm_boundary

!=======================================================================
Expand Down
32 changes: 17 additions & 15 deletions columnphysics/icepack_flux.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,10 +89,6 @@ subroutine merge_fluxes (aicen, &
fsaltn , & ! salt flux to ocean (kg/m2/s)
fhocnn , & ! actual ocn/ice heat flx (W/m**2)
fswthrun, & ! sw radiation through ice bot (W/m**2)
fswthrun_vdr, & ! vis dir sw radiation through ice bot (W/m**2)
fswthrun_vdf, & ! vis dif sw radiation through ice bot (W/m**2)
fswthrun_idr, & ! nir dir sw radiation through ice bot (W/m**2)
fswthrun_idf, & ! nir dif sw radiation through ice bot (W/m**2)
melttn , & ! top ice melt (m)
meltbn , & ! bottom ice melt (m)
meltsn , & ! snow melt (m)
Expand All @@ -102,6 +98,10 @@ subroutine merge_fluxes (aicen, &
snoicen ! snow-ice growth (m)

real (kind=dbl_kind), optional, intent(in):: &
fswthrun_vdr, & ! vis dir sw radiation through ice bot (W/m**2)
fswthrun_vdf, & ! vis dif sw radiation through ice bot (W/m**2)
fswthrun_idr, & ! nir dir sw radiation through ice bot (W/m**2)
fswthrun_idf, & ! nir dif sw radiation through ice bot (W/m**2)
Urefn ! air speed reference level (m/s)

! cumulative fluxes
Expand Down Expand Up @@ -129,7 +129,6 @@ subroutine merge_fluxes (aicen, &
meltb , & ! bottom ice melt (m)
melts , & ! snow melt (m)
meltsliq, & ! mass of snow melt (kg/m^2)
dsnow , & ! change in snow depth (m)
congel , & ! congelation ice growth (m)
snoice ! snow-ice growth (m)

Expand All @@ -139,18 +138,19 @@ subroutine merge_fluxes (aicen, &
fswthru_idr , & ! nir dir sw radiation through ice bot (W/m**2)
fswthru_idf ! nir dif sw radiation through ice bot (W/m**2)

real (kind=dbl_kind), optional, intent(inout):: &
real (kind=dbl_kind), intent(inout), optional :: &
dsnow, & ! change in snow depth (m)
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), dimension(:), intent(inout), optional :: &
Qref_iso, & ! isotope air sp hum ref 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)
real (kind=dbl_kind), dimension(:), intent(in), optional :: &
Qrefn_iso, & ! isotope air sp hum ref 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 @@ -220,7 +220,9 @@ subroutine merge_fluxes (aicen, &
if (tr_snow) then
meltsliq = meltsliq + meltsliqn * aicen
endif
dsnow = dsnow + dsnown * aicen
if (present(dsnow)) then
dsnow = dsnow + dsnown * aicen
endif
congel = congel + congeln * aicen
snoice = snoice + snoicen * aicen

Expand Down
22 changes: 13 additions & 9 deletions columnphysics/icepack_isotope.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,21 +83,24 @@ subroutine update_isotope (dt, &
aicen, & ! ice area
aice_old, & ! beginning values
vice_old, &
vsno_old, &
HDO_ocn, & !
H2_16O_ocn, & !
H2_18O_ocn !

real (kind=dbl_kind), dimension(:), intent(in) :: &
fiso_atm, & ! isotopic snowfall (kg/m^2/s of water)
Qref_iso ! isotope reference humidity
vsno_old

real (kind=dbl_kind), dimension(:), intent(inout) :: &
fiso_ocnn, & ! isotopic freshwater (kg/m^2/s)
fiso_evapn ! evaporative water flux (kg/m^2/s)

real (kind=dbl_kind), dimension(:), intent(inout) :: &
isosno, isoice ! mass of isotopes (kg)
isosno, & ! mass of isotopes (kg)
isoice

real (kind=dbl_kind), dimension(:), intent(in) :: &
fiso_atm, & ! isotopic snowfall (kg/m^2/s of water)
Qref_iso ! isotope reference humidity

real (kind=dbl_kind), intent(in) :: &
HDO_ocn, & ! ocean concentration of HDO (kg/kg)
H2_16O_ocn, & ! ocean concentration of H2_16O (kg/kg)
H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)

! local variables

Expand Down Expand Up @@ -197,6 +200,7 @@ subroutine update_isotope (dt, &

if (congel > c0) then
do k = 1,n_iso
work = c0
if (k == 1) then
alpha = isoice_alpha(congel/dt,'HDO',isotope_frac_method)
work = alpha*HDO_ocn*rhoi*congel*aicen
Expand Down
2 changes: 1 addition & 1 deletion columnphysics/icepack_itd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -833,7 +833,7 @@ subroutine cleanup_itd (dt, ntrcr, &
real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
faero_ocn ! aerosol flux to ocean (kg/m^2/s)

real (kind=dbl_kind), dimension (:), intent(inout) :: &
real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
fiso_ocn ! isotope flux to ocean (kg/m^2/s)

logical (kind=log_kind), intent(in), optional :: &
Expand Down
Loading