Skip to content

Commit d73163f

Browse files
josephzhang8josephzhang8
committed
Added checks in Fairall air-sea exchange.
Turned this option on as default from now on
1 parent 8f7bcf1 commit d73163f

File tree

2 files changed

+67
-51
lines changed

2 files changed

+67
-51
lines changed

‎cmake/SCHISM.local.build‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ set (OLDIO OFF CACHE BOOLEAN "Old nc output (each rank dumps its own data)")
2727
set (USE_ATMOS OFF CACHE BOOLEAN "Coupling with atmospheric model via ESMF")
2828
set (USE_BMI OFF CACHE BOOLEAN "Build a library that can be used from BMI, with sources and some b.c. provided by BMI")
2929
set (PREC_EVAP OFF CACHE BOOLEAN "Include precipitation and evaporation calculation")
30-
set (USE_BULK_FAIRALL OFF CACHE BOOLEAN "Enable Fairall bulk scheme for air-sea exchange")
30+
set (USE_BULK_FAIRALL ON CACHE BOOLEAN "Enable Fairall bulk scheme for air-sea exchange")
3131
## GOTM5.2 is installed by default, and has been tested successfully, defined GOTM_BASE to use a non-standard GOTM source
3232
## Note that source code of GOTM5.2 is from GOTM web so make sure to use recursive clone (to get the submodule)
3333
##set( GOTM_BASE /work2/03473/seatonc/frontera/GOTM5.2/code CACHE STRING "Path to non-standard GOTM base directory")

‎src/Hydro/sflux_9c.F90‎

Lines changed: 66 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -3241,14 +3241,15 @@ SUBROUTINE FAIRALL(num_nodes, &
32413241

32423242
use schism_glbl, only : rkind, uu2, vv2,tr_nd, & !tnd, snd, &
32433243
& idry, nvrt, ivcor,errmsg
3244-
use schism_glbl, only : grav
3244+
use schism_glbl, only : grav !GRAV
32453245
! use schism_glbl, only : rho0
32463246
use schism_msgp, only : myrank,parallel_abort
32473247

32483248
!================================================== Hernan G. Arango ===
32493249
! Copyright (c) 2002-2007 The ROMS/TOMS Group !
32503250
! Licensed under a MIT/X style license !
32513251
! See License_ROMS.txt !
3252+
! Edited by Jerome Lefevre and Joseph Zhang
32523253
!=======================================================================
32533254
! This routine computes the bulk parameterization of surface wind !
32543255
! stress and surface net heat fluxes. !
@@ -3354,7 +3355,7 @@ SUBROUTINE FAIRALL(num_nodes, &
33543355
! 4000, and varies only slightly (see Gill, 1982, Appendix 3).
33553356
REAL(rkind),parameter :: Cp=3985.0d0
33563357
! Functions:
3357-
REAL(rkind) :: bulk_psit,bulk_psiu
3358+
REAL(rkind) :: bulk_psit,bulk_psiu,tmp1,tmp2,tmp3
33583359

33593360

33603361
#ifdef DEBUG
@@ -3366,7 +3367,7 @@ SUBROUTINE FAIRALL(num_nodes, &
33663367

33673368
! set inverse of specific heat for seawater (kg-degC/Joule).
33683369
! cp is defined in scalars.h
3369-
cpi=1.0d0/cp
3370+
cpi=1.0d0/Cp
33703371
!
33713372
! Input bulk parameterization fields
33723373
!
@@ -3383,7 +3384,7 @@ SUBROUTINE FAIRALL(num_nodes, &
33833384
!$OMP rho0i,cpi,patmb,rhoSea,Qsea,TseaK,TseaC,TairC,TairK,rhoAir,Qair,Q, &
33843385
!$OMP VisAir,Hlv,delW,delT,delQ,u10,Zo10,Cd10,Ch10,Ct10,Cd,Ct,CC,Ri,Ribcu, &
33853386
!$OMP Zetu,L10,Wstar,Tstar,Qstar,ZoW,ZoT,ZoT10,ZoQ,ZoL,L,Rr,Bf,Wpsi,Tpsi, &
3386-
!$OMP Qpsi,Wgus,charn,upvel,evap,wspd0,dry,i_node, sfc_lev)
3387+
!$OMP Qpsi,Wgus,charn,upvel,evap,wspd0,dry,i_node, sfc_lev,tmp1,tmp2,tmp3)
33873388
do i_node = 1, num_nodes !=npa
33883389
!=================================================================
33893390
#ifdef DEBUG
@@ -3457,10 +3458,11 @@ SUBROUTINE FAIRALL(num_nodes, &
34573458
! Compute air saturation vapor pressure (mb), using Teten formula.
34583459
!
34593460
cff=(1.0007d0+3.46d-6*patmb)*6.1121d0*&
3460-
& exp(17.502d0*TairC/(240.97d0+TairC))
3461+
& exp(17.502d0*TairC/(240.97d0+TairC)) !check?
34613462
!
34623463
! Compute specific humidity at Saturation, Qair (kg/kg).
34633464
!
3465+
if(patmb-0.378d0*cff==0.d0) call parallel_abort('FAIRALL: Qair')
34643466
Qair=0.62197d0*(cff/(patmb-0.378d0*cff))
34653467
!
34663468
! Compute specific humidity, Q (kg/kg).
@@ -3478,14 +3480,15 @@ SUBROUTINE FAIRALL(num_nodes, &
34783480
! Compute water saturation vapor pressure (mb), using Teten formula.
34793481
!
34803482
cff=(1.0007d0+3.46d-6*patmb)*6.1121d0*&
3481-
& exp(17.502d0*TseaC/(240.97d0+TseaC))
3483+
& exp(17.502d0*TseaC/(240.97d0+TseaC)) !check?
34823484
!
34833485
! Vapor Pressure reduced for salinity (Kraus & Businger, 1994, pp 42).
34843486
!
34853487
cff=cff*0.98d0
34863488
!
34873489
! Compute Qsea (kg/kg) from vapor pressure.
34883490
!
3491+
if(patmb-0.378d0*cff==0.d0) call parallel_abort('FAIRALL: Qsea')
34893492
Qsea=0.62197d0*(cff/(patmb-0.378d0*cff))
34903493
!
34913494
!-----------------------------------------------------------------------
@@ -3495,8 +3498,7 @@ SUBROUTINE FAIRALL(num_nodes, &
34953498
!
34963499
! Moist air density (kg/m3).
34973500
!
3498-
rhoAir=patmb*100.0d0/(blk_Rgas*TairK*&
3499-
& (1.0d0+0.61d0*Q))
3501+
rhoAir=patmb*100.0d0/(blk_Rgas*TairK*(1.0d0+0.61d0*Q)) !Q>=0?
35003502
!
35013503
! Kinematic viscosity of dry air (m2/s), Andreas (1989).
35023504
!
@@ -3513,7 +3515,7 @@ SUBROUTINE FAIRALL(num_nodes, &
35133515
& v_air(i_node)*v_air(i_node) )
35143516
!
35153517
Wgus=0.5d0
3516-
delW=SQRT(wspd0*wspd0+Wgus*Wgus)
3518+
delW=SQRT(wspd0*wspd0+Wgus*Wgus) !>0
35173519
delQ=Qsea-Q
35183520
delT=TseaC-TairC
35193521

@@ -3529,27 +3531,32 @@ SUBROUTINE FAIRALL(num_nodes, &
35293531
! Neutral coefficients.
35303532
!
35313533
ZoW=1.0d-4
3532-
u10=delW*LOG(10.0d0/ZoW)/LOG(blk_ZW/ZoW)
3533-
Wstar=0.035d0*u10
3534+
u10=delW*LOG(10.0d0/ZoW)/LOG(blk_ZW/ZoW) !>0; actually blk_ZW=10
3535+
Wstar=0.035d0*u10 !>0
35343536
Zo10=0.011d0*Wstar*Wstar/GRAV+0.11d0*VisAir/Wstar
3535-
Cd10=(vonKar/LOG(10.0d0/Zo10))**2.d0
3537+
if(Zo10<=0.d0.or.Zo10==10.d0) call parallel_abort('FAIRALL: Zo10')
3538+
Cd10=(vonKar/LOG(10.0d0/Zo10))**2.d0 !>0; Zo10>0
35363539
Ch10=0.00115d0
3537-
Ct10=Ch10/sqrt(Cd10)
3538-
ZoT10=10.0d0/exp(vonKar/Ct10)
3539-
Cd=(vonKar/LOG(blk_ZW/Zo10))**2.d0
3540+
Ct10=Ch10/sqrt(Cd10) !>0
3541+
if(vonKar/Ct10>500.d0) call parallel_abort('FAIRALL: overflow, vonKar/Ct10>500')
3542+
ZoT10=10.0d0/exp(vonKar/Ct10) !>0
3543+
if(blk_ZW/Zo10<=0.d0.or.blk_ZW==Zo10.or.blk_ZT==ZoT10) call parallel_abort('FAIRALL: Zo10')
3544+
Cd=(vonKar/LOG(blk_ZW/Zo10))**2.d0 !>0
35403545
!
35413546
! Compute Richardson number.
3542-
!
3543-
Ct=vonKar/LOG(blk_ZT/ZoT10) ! T transfer coefficient
3544-
CC=vonKar*Ct/Cd
3545-
Ribcu=-blk_ZW/(blk_Zabl*0.004d0*blk_beta**3.d0)
3546-
Ri=-GRAV*blk_ZW*(delT+0.61d0*TairK*delQ)/&
3547-
(TairK*delW*delW)
3547+
! blk_ZT/ZoT10 checked above
3548+
Ct=vonKar/LOG(blk_ZT/ZoT10) ! T transfer coefficient; /=0
3549+
CC=vonKar*Ct/Cd !/=0
3550+
Ribcu=-blk_ZW/(blk_Zabl*0.004d0*blk_beta**3.d0) !/=0
3551+
Ri=-GRAV*blk_ZW*(delT+0.61d0*TairK*delQ)/(TairK*delW*delW) !delW>0
35483552
if (Ri.lt.0.0d0) then
3553+
if(1.0d0+Ri/Ribcu==0.d0) call parallel_abort('FAIRALL: Ri<0')
35493554
Zetu=CC*Ri/(1.0d0+Ri/Ribcu) ! Unstable
35503555
else
3556+
if(1.0d0+3.0d0*Ri/CC==0.d0) call parallel_abort('FAIRALL: Ri>=0')
35513557
Zetu=CC*Ri/(1.0d0+3.0d0*Ri/CC) ! Stable
35523558
endif
3559+
if(Zetu==0.d0) call parallel_abort('FAIRALL: Zetu=0')
35533560
L10=blk_ZW/Zetu
35543561
if (Zetu.gt.50.0d0) then
35553562
IterMax=1
@@ -3559,12 +3566,14 @@ SUBROUTINE FAIRALL(num_nodes, &
35593566
!
35603567
! First guesses for Monon-Obukhov similarity scales.
35613568
!
3562-
Wstar= delW*vonKar/(LOG(blk_ZW/Zo10)-&
3563-
& bulk_psiu(blk_ZW/L10))
3564-
Tstar=-delT*vonKar/(LOG(blk_ZT/ZoT10)-&
3565-
& bulk_psit(blk_ZT/L10))
3566-
Qstar=-delQ*vonKar/(LOG(blk_ZQ/ZoT10)-&
3567-
& bulk_psit(blk_ZQ/L10))
3569+
!Zo10>0
3570+
tmp1=LOG(blk_ZW/Zo10)-bulk_psiu(blk_ZW/L10)
3571+
tmp2=LOG(blk_ZT/ZoT10)-bulk_psit(blk_ZT/L10)
3572+
tmp3=LOG(blk_ZQ/ZoT10)-bulk_psit(blk_ZQ/L10)
3573+
if(tmp1==0.d0.or.tmp2==0.d0.or.tmp3==0.d0) call parallel_abort('FAIRALL: tmp=0(0)')
3574+
Wstar= delW*vonKar/tmp1 !(LOG(blk_ZW/Zo10)-bulk_psiu(blk_ZW/L10))
3575+
Tstar=-delT*vonKar/tmp2 !(LOG(blk_ZT/ZoT10)-bulk_psit(blk_ZT/L10))
3576+
Qstar=-delQ*vonKar/tmp3 !(LOG(blk_ZQ/ZoT10)-bulk_psit(blk_ZQ/L10))
35683577
!
35693578
! Modify Charnock for high wind speeds. The 0.125 factor below is for
35703579
! 1.0/(18.0-10.0).
@@ -3580,6 +3589,7 @@ SUBROUTINE FAIRALL(num_nodes, &
35803589
! Iterate until convergence. It usually converges within four
35813590
! iterations.
35823591
!
3592+
if(VisAir==0.d0) call parallel_abort('FAIRALL: VisAir=0')
35833593
do Iter=1,IterMax
35843594
ZoW=charn*Wstar*Wstar/GRAV+0.11d0*VisAir/(Wstar+eps)
35853595
Rr=ZoW*Wstar/VisAir
@@ -3591,7 +3601,8 @@ SUBROUTINE FAIRALL(num_nodes, &
35913601
ZoL=vonKar*GRAV*blk_ZW*&
35923602
& (Tstar*(1.0d0+0.61d0*Q)+0.61d0*TairK*Qstar)/&
35933603
& (TairK*Wstar*Wstar*(1.0d0+0.61d0*Q)+eps)
3594-
L=blk_ZW/(ZoL+eps)
3604+
if(ZoL+eps==0.d0) call parallel_abort('FAIRALL: ZoL+eps=0')
3605+
L=blk_ZW/(ZoL+eps) !/=0
35953606
!
35963607
! Evaluate stability functions at Z/L.
35973608
!
@@ -3601,20 +3612,25 @@ SUBROUTINE FAIRALL(num_nodes, &
36013612
!
36023613
! Compute wind scaling parameters, Wstar.
36033614
!
3604-
Wstar=MAX(eps,delW*vonKar/(LOG(blk_ZW/ZoW)-Wpsi))
3605-
Tstar=-delT*vonKar/(LOG(blk_ZT/ZoT)-Tpsi)
3606-
Qstar=-delQ*vonKar/(LOG(blk_ZQ/ZoQ)-Qpsi)
3615+
if(ZoW<=0.d0) call parallel_abort('FAIRALL: ZoW<=0')
3616+
tmp1=LOG(blk_ZW/ZoW)-Wpsi
3617+
tmp2=LOG(blk_ZT/ZoT)-Tpsi
3618+
tmp3=LOG(blk_ZQ/ZoQ)-Qpsi
3619+
if(tmp1==0.d0.or.tmp2==0.d0.or.tmp3==0.d0) call parallel_abort('FAIRALL: tmp=0')
3620+
Wstar=MAX(eps,delW*vonKar/tmp1) !(LOG(blk_ZW/ZoW)-Wpsi)) !>0
3621+
Tstar=-delT*vonKar/tmp2 !(LOG(blk_ZT/ZoT)-Tpsi)
3622+
Qstar=-delQ*vonKar/tmp3 !(LOG(blk_ZQ/ZoQ)-Qpsi)
36073623
!
36083624
! Compute gustiness in wind speed.
36093625
!
36103626
Bf=-GRAV/TairK*Wstar*(Tstar+0.61d0*TairK*Qstar)
36113627
if (Bf.gt.0.0d0) then
3612-
Wgus=blk_beta*(Bf*blk_Zabl)**r3
3628+
Wgus=blk_beta*(Bf*blk_Zabl)**r3 !>0
36133629
else
36143630
Wgus=0.2d0
36153631
endif
3616-
delW=SQRT(wspd0*wspd0+Wgus*Wgus)
3617-
enddo
3632+
delW=SQRT(wspd0*wspd0+Wgus*Wgus) !>0
3633+
enddo !Iter
36183634
!
36193635
!-----------------------------------------------------------------------
36203636
! Compute Atmosphere/Ocean fluxes.
@@ -3623,8 +3639,8 @@ SUBROUTINE FAIRALL(num_nodes, &
36233639
!
36243640
! Compute transfer coefficients for momentum (Cd).
36253641
!
3626-
wspd=SQRT(wspd0*wspd0+Wgus*Wgus)
3627-
Cd=Wstar*Wstar/(wspd*wspd+eps)
3642+
wspd=SQRT(wspd0*wspd0+Wgus*Wgus) !>0
3643+
Cd=Wstar*Wstar/(wspd*wspd+eps) !>0
36283644
!
36293645
! Compute turbulent sensible heat flux (W/m2), Hs.
36303646
!
@@ -3712,6 +3728,7 @@ END SUBROUTINE FAIRALL
37123728
REAL(rkind) FUNCTION bulk_psiu(ZoL)
37133729

37143730
use schism_glbl, only : rkind
3731+
use schism_msgp, only : myrank,parallel_abort
37153732

37163733
IMPLICIT NONE
37173734
!
@@ -3749,29 +3766,28 @@ REAL(rkind) FUNCTION bulk_psiu(ZoL)
37493766
! Unstable conditions.
37503767
!
37513768
if (ZoL.lt.0.0d0) then
3752-
x=(1.0d0-15.0d0*ZoL)**0.25d0
3769+
x=(1.0d0-15.0d0*ZoL)**0.25d0 !>0
37533770
psik=2.0d0*LOG(0.5d0*(1.0d0+x))+ &
37543771
& LOG(0.5d0*(1.0d0+x*x))-2.0d0*ATAN(x)+0.5d0*PI
37553772
!
37563773
! For very unstable conditions, use free-convection (Fairall).
37573774
!
37583775
cff=SQRT(3.0d0)
3759-
y=(1.0d0-10.15d0*ZoL)**r3
3776+
y=(1.0d0-10.15d0*ZoL)**r3 !>0
37603777
psic=1.5d0*LOG(r3*(1.0d0+y+y*y))- &
37613778
& cff*ATAN((1.0d0+2.0d0*y)/cff)+PI/cff
37623779
!
37633780
! Match Kansas and free-convection forms with weighting Fw.
37643781
!
3765-
cff=ZoL*ZoL
3766-
Fw=cff/(1.0d0+cff)
3782+
cff=ZoL*ZoL !>0
3783+
Fw=cff/(1.0d0+cff) !\in (0,1)
37673784
bulk_psiu=(1.0d0-Fw)*psik+Fw*psic
37683785
!
37693786
! Stable conditions.
37703787
!
3771-
else
3772-
cff=MIN(50.0d0,0.35d0*ZoL)
3773-
bulk_psiu=-((1.0d0+ZoL)+ &
3774-
& 0.6667d0*(ZoL-14.28d0)/EXP(cff)+8.525d0)
3788+
else !ZoL>=0
3789+
cff=MIN(50.0d0,0.35d0*ZoL) !>=0
3790+
bulk_psiu=-((1.0d0+ZoL)+0.6667d0*(ZoL-14.28d0)/EXP(cff)+8.525d0)
37753791
endif
37763792
return
37773793
END FUNCTION bulk_psiu
@@ -3813,26 +3829,26 @@ REAL(rkind) FUNCTION bulk_psit(ZoL)
38133829
! Unstable conditions.
38143830
!
38153831
if (ZoL.lt.0.0d0) then
3816-
x=(1.0d0-15.0d0*ZoL)**0.5d0
3832+
x=(1.0d0-15.0d0*ZoL)**0.5d0 !>0
38173833
psik=2.0d0*LOG(0.5d0*(1.0d0+x))
38183834
!
38193835
! For very unstable conditions, use free-convection (Fairall).
38203836
!
38213837
cff=SQRT(3.0d0)
3822-
y=(1.0-34.15*ZoL)**r3
3838+
y=(1.0-34.15*ZoL)**r3 !>0
38233839
psic=1.5d0*LOG(r3*(1.0d0+y+y*y))- &
38243840
& cff*ATAN((1.0d0+2.0d0*y)/cff)+PI/cff
38253841
!
38263842
! Match Kansas and free-convection forms with weighting Fw.
38273843
!
3828-
cff=ZoL*ZoL
3829-
Fw=cff/(1.0d0+cff)
3844+
cff=ZoL*ZoL !>0
3845+
Fw=cff/(1.0d0+cff) !\in (0,1)
38303846
bulk_psit=(1.0d0-Fw)*psik+Fw*psic
38313847
!
38323848
! Stable conditions.
38333849
!
3834-
else
3835-
cff=MIN(50.0d0,0.35d0*ZoL)
3850+
else !ZoL>=0
3851+
cff=MIN(50.0d0,0.35d0*ZoL) !>=0
38363852
bulk_psit=-((1.0d0+2.0d0*ZoL)**1.5d0+ &
38373853
& 0.6667d0*(ZoL-14.28d0)/EXP(cff)+8.525d0)
38383854
endif

0 commit comments

Comments
 (0)