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

update icepack; switch to "no abort" when inconsistent ITD condition is found #4

Merged
merged 13 commits into from
Jun 4, 2021
Merged

update icepack; switch to "no abort" when inconsistent ITD condition is found #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
ee7cf8e
modify error tolerance for rare mushy thermo failure (#352)
eclare108213 Mar 17, 2021
7223366
Added catch for unusual case of (hicen_init(n+1) - hicen_init(n))>0 (…
proteanplanet Mar 17, 2021
e8df0ad
Addition of CompyMcNodeFace to Icepack machine files (#353)
proteanplanet Mar 17, 2021
fa15366
update Compy (#356)
apcraig Mar 18, 2021
f298152
Documentation on adding diagnostic output. (#331)
dabail10 Mar 18, 2021
c027a3c
move shortwave absorbed in negligible snow layers to surface (#355)
eclare108213 Mar 23, 2021
ed0a991
Update version number and revise reference (#358)
apcraig Mar 27, 2021
5cf2232
Update test results color implementation to point to files for colors…
apcraig Apr 7, 2021
43a909a
Add '.giattributes' and define diff attribute for Fortran files (#362)
phil-blain May 24, 2021
77c4e5f
machines: add ECCC machines 'daley' and 'banting' (#363)
phil-blain May 26, 2021
eb59c95
clean up ampersand continuation lines (#365)
apcraig May 27, 2021
37f2a17
Add support for staggered atmospheric levels (#364)
phil-blain May 28, 2021
2845c94
switch abort to false
DeniseWorthen Jun 2, 2021
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
1 change: 1 addition & 0 deletions .gitattributes
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
*.F90 diff=fortran
10 changes: 5 additions & 5 deletions .github/workflows/test-icepack.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ on:

defaults:
run:
shell: /bin/csh {0}
shell: /bin/csh -e {0}

jobs:
build:
Expand Down Expand Up @@ -83,9 +83,9 @@ jobs:
conda env create -f configuration/scripts/machines/environment.yml
- name: check conda env
run: |
conda activate icepack && which mpicc && which mpifort && which make
mpifort --version
mpicc --version
conda activate icepack && which clang && which gfortran && which make
gfortran --version
clang --version
make --version
- name: check setup case
run: |
Expand All @@ -104,7 +104,7 @@ jobs:
- name: download input data
run: |
cd $HOME/icepack-dirs/input
wget https://zenodo.org/record/3728287/files/Icepack_data-20200326.tar.gz && tar xvfz ICEPACK_data-20200326.tar.gz
wget --progress=dot:giga https://zenodo.org/record/3728287/files/Icepack_data-20200326.tar.gz && tar xvfz ICEPACK_data-20200326.tar.gz
pwd
ls -alR
# - name: run case
Expand Down
182 changes: 136 additions & 46 deletions columnphysics/icepack_atmo.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ subroutine atmo_boundary_layer (sfctype, &
Qa_iso, Qref_iso, &
iso_flag, &
uvel, vvel, &
Uref )
Uref, zlvs )

use icepack_parameters, only: highfreq, natmiter, atmiter_conv

Expand Down Expand Up @@ -118,14 +118,16 @@ subroutine atmo_boundary_layer (sfctype, &
real (kind=dbl_kind), intent(out) :: &
Uref ! reference height wind speed (m/s)

real (kind=dbl_kind), intent(in), optional :: &
zlvs ! atm level height (scalar quantities) (m)

! local variables

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

real (kind=dbl_kind) :: &
TsfK , & ! surface temperature in Kelvin (K)
xqq , & ! temporary variable
psimh , & ! stability function at zlvl (momentum)
tau , & ! stress at zlvl
fac , & ! interpolation factor
Expand Down Expand Up @@ -153,13 +155,15 @@ subroutine atmo_boundary_layer (sfctype, &
re , & ! sqrt of exchange coefficient (water)
rh , & ! sqrt of exchange coefficient (heat)
vmag , & ! surface wind magnitude (m/s)
alz , & ! ln(zlvl /z10)
alzm , & ! ln(zlvl /z10)
alzs , & ! ln(zlvs /z10) (if zlvs present)
thva , & ! virtual temperature (K)
cp , & ! specific heat of moist air
hol , & ! H (at zlvl ) over L
holm , & ! H (at zlvl ) over L
hols , & ! H (at zlvs ) over L (if zlvs present)
stable, & ! stability factor
cpvir , & ! defined as cp_wv/cp_air - 1.
psixh ! stability function at zlvl (heat and water)
psixh ! stability function at zlvl (at zlvs if present) (heat and water)

real (kind=dbl_kind), parameter :: &
zTrf = c2 ! reference height for air temp (m)
Expand Down Expand Up @@ -247,7 +251,12 @@ subroutine atmo_boundary_layer (sfctype, &

thva = potT * (c1 + zvir * Qa) ! virtual pot temp (K)
delq = Qa - ssq ! spec hum dif (kg/kg)
alz = log(zlvl/zref)
alzm = log(zlvl/zref)
if (present(zlvs)) then
alzs = log(zlvs/zref)
else
alzs = alzm
endif
cp = cp_air*(c1 + cpvir*ssq)

!------------------------------------------------------------
Expand Down Expand Up @@ -275,28 +284,24 @@ subroutine atmo_boundary_layer (sfctype, &
ustar_prev = ustar

! compute stability & evaluate all stability functions
hol = vonkar * gravit * zlvl &
* (tstar/thva &
+ qstar/(c1/zvir+Qa)) &
/ ustar**2
hol = sign( min(abs(hol),c10), hol )
stable = p5 + sign(p5 , hol)
xqq = max(sqrt(abs(c1 - c16*hol)) , c1)
xqq = sqrt(xqq)

! Jordan et al 1999
psimhs = -(0.7_dbl_kind*hol &
+ 0.75_dbl_kind*(hol-14.3_dbl_kind) &
* exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
psimh = psimhs*stable &
+ (c1 - stable)*psimhu(xqq)
psixh = psimhs*stable &
+ (c1 - stable)*psixhu(xqq)
holm = compute_stability_parameter(zlvl , thva , &
ustar, tstar, &
qstar, Qa)
if (present(zlvs)) then
hols = compute_stability_parameter(zlvs , thva , &
ustar, tstar, &
qstar, Qa)
else
hols = holm
endif

call compute_stability_function('momentum', holm, stable, psimh)
call compute_stability_function('scalar' , hols, stable, psixh)

! shift all coeffs to measurement height and stability
rd = rdn / (c1+rdn/vonkar*(alz-psimh))
rh = rhn / (c1+rhn/vonkar*(alz-psixh))
re = ren / (c1+ren/vonkar*(alz-psixh))
rd = rdn / (c1+rdn/vonkar*(alzm-psimh))
rh = rhn / (c1+rhn/vonkar*(alzs-psixh))
re = ren / (c1+ren/vonkar*(alzs-psixh))

! update ustar, tstar, qstar using updated, shifted coeffs
ustar = rd * vmag
Expand Down Expand Up @@ -361,16 +366,14 @@ subroutine atmo_boundary_layer (sfctype, &
! Compute diagnostics: 2m ref T, Q, U
!------------------------------------------------------------

hol = hol*zTrf/zlvl
xqq = max( c1, sqrt(abs(c1-c16*hol)) )
xqq = sqrt(xqq)
psix2 = -c5*hol*stable + (c1-stable)*psixhu(xqq)
hols = hols*zTrf/zlvl
psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols)
fac = (rh/vonkar) &
* (alz + al2 - psixh + psix2)
* (alzs + al2 - psixh + psix2)
Tref = potT - delt*fac
Tref = Tref - p01*zTrf ! pot temp to temp correction
fac = (re/vonkar) &
* (alz + al2 - psixh + psix2)
* (alzs + al2 - psixh + psix2)
Qref = Qa - delq*fac

if (highfreq .and. sfctype(1:3)=='ice') then
Expand Down Expand Up @@ -853,7 +856,7 @@ subroutine icepack_atm_boundary(sfctype, &
Cdn_atm_ratio_n, &
Qa_iso, Qref_iso, &
uvel, vvel, &
Uref)
Uref, zlvs)

character (len=3), intent(in) :: &
sfctype ! ice or ocean
Expand All @@ -864,7 +867,7 @@ subroutine icepack_atm_boundary(sfctype, &
uatm , & ! x-direction wind speed (m/s)
vatm , & ! y-direction wind speed (m/s)
wind , & ! wind speed (m/s)
zlvl , & ! atm level height (m)
zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m)
Qa , & ! specific humidity (kg/kg)
rhoa ! air density (kg/m^3)

Expand Down Expand Up @@ -893,7 +896,8 @@ subroutine icepack_atm_boundary(sfctype, &

real (kind=dbl_kind), optional, intent(in) :: &
uvel , & ! x-direction ice speed (m/s)
vvel ! y-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) :: &
Uref ! reference height wind speed (m/s)
Expand Down Expand Up @@ -965,7 +969,7 @@ subroutine icepack_atm_boundary(sfctype, &
Qa_iso=l_Qa_iso, &
Qref_iso=l_Qref_iso, &
uvel=l_uvel, vvel=l_vvel, &
Uref=l_Uref)
Uref=l_Uref, zlvs=zlvs)
if (icepack_warnings_aborted(subname)) return
endif ! atmbndy

Expand All @@ -981,31 +985,117 @@ subroutine icepack_atm_boundary(sfctype, &

end subroutine icepack_atm_boundary

!=======================================================================

function compute_stability_parameter(zlvl , thva , &
ustar, tstar, &
qstar, Qa) &
result(hol)

real (kind=dbl_kind), intent(in) :: &
zlvl , & ! atm level height (m)
thva , & ! virtual temperature (K)
ustar , & ! turbulent scale for momentum
tstar , & ! turbulent scale for temperature
qstar , & ! turbulent scale for humidity
Qa ! specific humidity (kg/kg)

real (kind=dbl_kind) :: &
hol ! H (at zlvl) over L

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

hol = vonkar * gravit * zlvl &
* (tstar/thva &
+ qstar/(c1/zvir+Qa)) &
/ ustar**2
hol = sign( min(abs(hol),c10), hol)

end function compute_stability_parameter

!=======================================================================

subroutine compute_stability_function(qty, hol, stable, psi)

character (len=*), intent(in) :: &
qty ! 'momentum' or 'scalar'

real (kind=dbl_kind), intent(in) :: &
hol ! H over L

real (kind=dbl_kind), intent(out) :: &
psi , & ! stability function at hol
stable ! unit step function at hol

! local variables

real (kind=dbl_kind) :: &
psi_stable , & ! stable stability funcion at hol
psi_unstable ! unstable stability funcion at hol

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

stable = p5 + sign(p5 , hol)

psi_stable = -(0.7_dbl_kind*hol &
+ 0.75_dbl_kind*(hol-14.3_dbl_kind) &
* exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)

if(trim(qty) == 'momentum') then
psi_unstable = psi_momentum_unstable(hol)
elseif(trim(qty) == 'scalar') then
psi_unstable = psi_scalar_unstable(hol)
else
call icepack_warnings_add(subname//' incorrect qty: ' // qty)
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
endif

psi = psi_stable*stable + (c1 - stable)*psi_unstable

end subroutine compute_stability_function

!------------------------------------------------------------
! Define functions
!------------------------------------------------------------

!=======================================================================

real(kind=dbl_kind) function psimhu(xd)
real(kind=dbl_kind) function psi_momentum_unstable(hol)

real(kind=dbl_kind), intent(in) :: hol

real(kind=dbl_kind) :: xd

xd = capital_X(hol)

psi_momentum_unstable = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
- c2*atan(xd) + pih

end function psi_momentum_unstable

!=======================================================================

real(kind=dbl_kind) function psi_scalar_unstable(hol)

real(kind=dbl_kind), intent(in) :: hol

real(kind=dbl_kind) :: xd

real(kind=dbl_kind), intent(in) :: xd
xd = capital_X(hol)

psimhu = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
- c2*atan(xd) + pih
!ech - c2*atan(xd) + 1.571_dbl_kind
psi_scalar_unstable = c2 * log((c1 + xd*xd)/c2)

end function psimhu
end function psi_scalar_unstable

!=======================================================================

real(kind=dbl_kind) function psixhu(xd)
real(kind=dbl_kind) function capital_X(hol)

real(kind=dbl_kind), intent(in) :: xd
real(kind=dbl_kind), intent(in) :: hol

psixhu = c2 * log((c1 + xd*xd)/c2)
capital_X = sqrt(max(sqrt(abs(c1 - c16*hol)) , c1))

end function psixhu
end function capital_X

!=======================================================================

Expand Down
Loading