diff --git a/ben02/mod_ben02.F b/ben02/mod_ben02.F
index d6fd01a5..2e6c63cb 100644
--- a/ben02/mod_ben02.F
+++ b/ben02/mod_ben02.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2002-2021 Mats Bentsen, Mehmet Ilicak
+! Copyright (C) 2002-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -26,7 +26,7 @@ module mod_ben02
c
use mod_types, only: i2, r4
use mod_config, only: expcnf
- use mod_constants, only: t0deg, spval
+ use mod_constants, only: t0deg, spval, L_mks2cgs
use mod_calendar, only: date_offset, calendar_noerr,
. calendar_errstr
use mod_time, only: date, calendar, nday_in_year, nday_of_year,
@@ -183,10 +183,17 @@ module mod_ben02
. atm_cswa_era ! short-wave radiation adjustment factor
! (NCEP)
c
+#ifdef MKS
+ data atm_ice_csmt_ncep,atm_rnf_csmt_ncep /2.e10,1.e9/,
+ . atm_crnf_ncep,atm_cswa_ncep /0.82073,0.88340/,
+ . atm_ice_csmt_era,atm_rnf_csmt_era /0.0,1.e9/,
+ . atm_crnf_era,atm_cswa_era /0.7234,0.9721/
+#else
data atm_ice_csmt_ncep,atm_rnf_csmt_ncep /2.e14,1.e13/,
. atm_crnf_ncep,atm_cswa_ncep /0.82073,0.88340/,
. atm_ice_csmt_era,atm_rnf_csmt_era /0.0,1.e13/,
. atm_crnf_era,atm_cswa_era /0.7234,0.9721/
+#endif
c
real ::
. zu, ! measurement height of wind [m]
@@ -2090,11 +2097,13 @@ subroutine inifrc_ben02clim
integer, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,12) :: smtmsk
real dx2,dy2,prc_sum,eva_sum,rnf_sum,swa_sum,lwa_sum,lht_sum,
. sht_sum,fwf_fac,dangle,garea,le,albedo,fac,swa_ave,lwa_ave,
- . lht_ave,sht_ave,crnf,cswa
+ . lht_ave,sht_ave,crnf,cswa,A_cgs2mks
real*4 rw4
integer i,j,k,l,il,jl
integer*2 rn2,ri2,rj2
c
+ A_cgs2mks=1./(L_mks2cgs**2)
+c
c --- Allocate memory for additional monthly forcing fields.
allocate(taud (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,12),
. tauxd (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,12),
@@ -2775,7 +2784,7 @@ subroutine inifrc_ben02clim
do k=1,12
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- garea=scp2(i,j)*1.e-4 ! [m^2]
+ garea=scp2(i,j)*A_cgs2mks ! [m^2]
c
c --- ----- freshwater fluxes [m/s]
util1(i,j)=util1(i,j)+precip(i,j,k)*fwf_fac*garea
@@ -2819,7 +2828,7 @@ subroutine inifrc_ben02clim
do j=1,jj
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- garea=scp2(i,j)*1.e-4 ! [m^2]
+ garea=scp2(i,j)*A_cgs2mks ! [m^2]
c
c --- ----- heat fluxes
albedo=albs_f*ricclm(i,j,k)+albw(i,j)*(1.-ricclm(i,j,k))
@@ -2838,7 +2847,7 @@ subroutine inifrc_ben02clim
call xcsum(lht_sum,util3,ip)
call xcsum(sht_sum,util4,ip)
c
- fac=1.e4/(12.*area)
+ fac=(L_mks2cgs*L_mks2cgs)/(12.*area)
swa_ave=swa_sum*fac
lwa_ave=lwa_sum*fac
lht_ave=lht_sum*fac
diff --git a/ben02/sfcstr_ben02.F b/ben02/sfcstr_ben02.F
index efd9e014..ee8161a6 100644
--- a/ben02/sfcstr_ben02.F
+++ b/ben02/sfcstr_ben02.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2004-2020 Mats Bentsen
+! Copyright (C) 2004-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,6 +24,7 @@ subroutine sfcstr_ben02(m,n,mm,nn,k1m,k1n)
c --- ------------------------------------------------------------------
c
use mod_xc
+ use mod_constants, only: P_mks2cgs
use mod_forcing, only: ztx, mty, taux, tauy
use mod_seaice, only: ficem, hicem, tauxice, tauyice
use mod_checksum, only: csdiag, chksummsk
@@ -44,14 +45,14 @@ subroutine sfcstr_ben02(m,n,mm,nn,k1m,k1n)
do i=max(1,ifu(j,l)),min(ii,ilu(j,l))
facice=(ficem(i,j)+ficem(i-1,j))
. *min(2.,hicem(i,j)+hicem(i-1,j))*.25
- taux(i,j)=10.*(ztx(i,j)*(1.-facice)+tauxice(i,j)*facice)
+ taux(i,j)=P_mks2cgs*(ztx(i,j)*(1.-facice)+tauxice(i,j)*facice)
enddo
enddo
do l=1,isv(j)
do i=max(1,ifv(j,l)),min(ii,ilv(j,l))
facice=(ficem(i,j)+ficem(i,j-1))
. *min(2.,hicem(i,j)+hicem(i,j-1))*.25
- tauy(i,j)=10.*(mty(i,j)*(1.-facice)+tauyice(i,j)*facice)
+ tauy(i,j)=P_mks2cgs*(mty(i,j)*(1.-facice)+tauyice(i,j)*facice)
enddo
enddo
enddo
diff --git a/ben02/thermf_ben02.F b/ben02/thermf_ben02.F
index 003e9c61..fcd5bd19 100644
--- a/ben02/thermf_ben02.F
+++ b/ben02/thermf_ben02.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2002-2021 Mats Bentsen
+! Copyright (C) 2002-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -21,7 +21,8 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c
c --- NERSC version of thermf.
c
- use mod_constants, only: spcifh, t0deg, epsil, onem
+ use mod_constants, only: spcifh, t0deg, alpha0, epsilt, onem,
+ . g2kg, kg2g, L_mks2cgs, M_mks2cgs
use mod_time, only: nday_in_year, nday_of_year, nstep,
. nstep_in_day, baclin,
. xmi, l1mi, l2mi, l3mi, l4mi, l5mi
@@ -71,7 +72,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
. fice,hice,hsnw,tsrf,fice0,hice0,hsnw0,qsww,qnsw,tice,albi,
. tsmlt,albi_h,qswi,dh,qsnwf,fcond,qdamp,qsmlt,qo2i,qbot,swfac,
. dtml,q,volice,df,dvi,dvs,fwflx,sstc,rice,trxflx,sssc,srxflx,
- . totsfl,totwfl,sflxc,totsrp,totsrn
+ . totsfl,totwfl,sflxc,totsrp,totsrn,A_cgs2mks
#ifdef TRC
integer nt
real, dimension(ntr,1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
@@ -82,10 +83,12 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c
real intp1d
external intp1d
+c
+ A_cgs2mks=1./(L_mks2cgs**2)
c
c --- Due to conservation, the ratio of ice and snow density must be
c --- equal to the ratio of ice and snow heat of fusion
- if (abs(fuss/fusi-rhosnw/rhoice).gt.epsil) then
+ if (abs(fuss/fusi-rhosnw/rhoice).gt.epsilt) then
if (mnproc.eq.1) then
write (lp,*)
. 'thermf: check consistency between snow/ice densities'
@@ -97,7 +100,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c
c --- Set various constants
dt=baclin ! Time step
- cpsw=spcifh*1.e3 ! Specific heat of seawater
+ cpsw=spcifh*M_mks2cgs ! Specific heat of seawater
rnf_fac=baclin/real(nrfets*86400) ! Runoff reservoar detrainment rate
sag_fac=exp(-sagets*dt) ! Snow aging rate
c
@@ -326,7 +329,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c --- ----- Ice volume that has to freeze to balance the heat budget
volice=-(qsww+qnsw-q)*(1.-fice)*dt/fusi
c
- if (volice.gt.epsil) then
+ if (volice.gt.epsilt) then
c
c --- ------- New ice in the lead is formed with a specified thickness.
c --- ------- Estimate the change in ice fraction
@@ -344,7 +347,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c --- ----- If the lead is warming, let the fraction (1 - fice) go to
c --- ----- warm the lead, and the fraction fice to melt ice laterally
fice=fice-(swfac*qsww+qnsw)*fice*dt
- . /max(hice*fusi+hsnw*fuss,epsil)
+ . /max(hice*fusi+hsnw*fuss,epsilt)
if (fice.lt.0.) then
fice=0.
hice=0.
@@ -398,14 +401,14 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
fwflx=eva(i,j)+lip(i,j)+sop(i,j)+rnf(i,j)+rfi(i,j)+fmltfz(i,j)
c
c --- --- Salt flux [kg m-2 s-1] (positive downwards)
- sfl(i,j)=-sice*dvi*rhoice/dt*1.e-3
+ sfl(i,j)=-sice*dvi*rhoice/dt*g2kg
c
c --- --- Salt flux due to brine rejection of freezing sea
c --- --- ice [kg m-2 m-1] (positive downwards)
- brnflx(i,j)=max(0.,-sotl*fmltfz(i,j)*1.e-3+sfl(i,j))
+ brnflx(i,j)=max(0.,-sotl*fmltfz(i,j)*g2kg+sfl(i,j))
c
c --- --- Virtual salt flux [kg m-2 s-1] (positive downwards)
- vrtsfl(i,j)=-sotl*fwflx*1.e-3
+ vrtsfl(i,j)=-sotl*fwflx*g2kg
c
c --- --- Store area weighted virtual salt flux and fresh water flux
util1(i,j)=vrtsfl(i,j)*scp2(i,j)
@@ -415,11 +418,11 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
hmltfz(i,j)=(dvi*fusi+dvs*fuss)/dt
c
c --- --- Total heat flux in BLOM units [W cm-2] (positive upwards)
- surflx(i,j)=-(swa(i,j)+nsf(i,j)+hmltfz(i,j))*1.e-4
+ surflx(i,j)=-(swa(i,j)+nsf(i,j)+hmltfz(i,j))*A_cgs2mks
c
c --- --- Short-wave heat flux in BLOM units [W cm-2] (positive
c --- --- upwards)
- sswflx(i,j)=-qsww*(1.-fice0)*1.e-4
+ sswflx(i,j)=-qsww*(1.-fice0)*A_cgs2mks
c
#ifdef TRC
c --- ------------------------------------------------------------------
@@ -452,7 +455,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
endif
# endif
# endif
- trflx(nt,i,j)=-trc(i,j,k1n,nt)*fwflx*1.e-3
+ trflx(nt,i,j)=-trc(i,j,k1n,nt)*fwflx*g2kg
ttrsf(nt,i,j)=trflx(nt,i,j)*scp2(i,j)
ttrav(nt,i,j)=trc(i,j,k1n,nt)*scp2(i,j)
enddo
@@ -465,7 +468,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
surrlx(i,j)=0.
c
c --- --- If trxday>0 , apply relaxation towards observed sst
- if (trxday.gt.epsil) then
+ if (trxday.gt.epsilt) then
sstc=intp1d(sstclm(i,j,l1mi),sstclm(i,j,l2mi),
. sstclm(i,j,l3mi),sstclm(i,j,l4mi),
. sstclm(i,j,l5mi),xmi)
@@ -473,8 +476,8 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
. ricclm(i,j,l3mi),ricclm(i,j,l4mi),
. ricclm(i,j,l5mi),xmi)
sstc=(1.-rice)*max(sstc,tice_f)+rice*tice_f
- trxflx=spcifh*100.*min(hmxl,trxdpt)/(trxday*86400.)
- . *min(trxlim,max(-trxlim,sstc-tmxl))
+ trxflx=spcifh*L_mks2cgs*min(hmxl,trxdpt)/(trxday*86400.)
+ . *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
surrlx(i,j)=-trxflx
else
trxflx=0.
@@ -496,12 +499,12 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
salrlx(i,j)=0.
c
c --- --- if srxday>0 , apply relaxation towards observed sss
- if (srxday.gt.epsil) then
+ if (srxday.gt.epsilt) then
sssc=intp1d(sssclm(i,j,l1mi),sssclm(i,j,l2mi),
. sssclm(i,j,l3mi),sssclm(i,j,l4mi),
. sssclm(i,j,l5mi),xmi)
- srxflx=100.*min(hmxl,srxdpt)/(srxday*86400.)
- . *min(srxlim,max(-srxlim,sssc-smxl))
+ srxflx=L_mks2cgs*min(hmxl,srxdpt)/(srxday*86400.)
+ . *min(srxlim,max(-srxlim,sssc-smxl))/alpha0
salrlx(i,j)=-srxflx
util3(i,j)=max(0.,salrlx(i,j))*scp2(i,j)
util4(i,j)=min(0.,salrlx(i,j))*scp2(i,j)
@@ -538,7 +541,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c --- -------------------------------------------------------------------
c
ustar(i,j)=(min(ustari(i,j),.8e-2)*fice0
- . +ustarw(i,j)*(1.-fice0))*1.e2
+ . +ustarw(i,j)*(1.-fice0))*L_mks2cgs
c
enddo
enddo
@@ -556,7 +559,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
call xcsum(totwfl,util2,ips)
c
c --- Correction for the virtual salt flux [kg m-2 s-1]
- sflxc=(-sref*totwfl*1.e-3-totsfl)/area
+ sflxc=(-sref*totwfl*g2kg-totsfl)/area
if (mnproc.eq.1) then
write (lp,*) 'thermf: totsfl/area,sflxc',totsfl/area,sflxc
endif
@@ -567,8 +570,10 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
do j=1,jj
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- salflx(i,j)=-(vrtsfl(i,j)+sflxc+sfl(i,j))*1.e2
- brnflx(i,j)=-brnflx(i,j)*1.e2
+ salflx(i,j)=-(vrtsfl(i,j)+sflxc+sfl(i,j))
+ . *(kg2g*(M_mks2cgs/L_mks2cgs**2))
+ brnflx(i,j)=-brnflx(i,j)
+ . *(kg2g*(M_mks2cgs/L_mks2cgs**2))
enddo
enddo
enddo
@@ -577,7 +582,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c --- if srxday>0 and srxbal=.true. , balance the sss relaxation flux
c --- so the net input of salt in grid cells connected to the world
c --- ocean is zero
- if (srxday.gt.epsil.and.srxbal) then
+ if (srxday.gt.epsilt.and.srxbal) then
call xcsum(totsrp,util3,ipwocn)
call xcsum(totsrn,util4,ipwocn)
if (abs(totsrp).gt.abs(totsrn)) then
@@ -632,14 +637,14 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
c tottrav=tottrav/area
cc
c trflxc=(-tottrsf)/area
-c trflxc=(-tottrav*totwfl*1.e-3-tottrsf)/area
+c trflxc=(-tottrav*totwfl*g2kg-tottrsf)/area
trflxc=0.
c
c$OMP PARALLEL DO PRIVATE(l,i)
do j=1,jj
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- trflx(nt,i,j)=-(trflx(nt,i,j)+trflxc)*1.e2
+ trflx(nt,i,j)=-(trflx(nt,i,j)+trflxc)*L_mks2cgs
enddo
enddo
enddo
diff --git a/cesm/sfcstr_cesm.F b/cesm/sfcstr_cesm.F
index b352cbfd..d0d047b7 100644
--- a/cesm/sfcstr_cesm.F
+++ b/cesm/sfcstr_cesm.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2020 Mats Bentsen
+! Copyright (C) 2015-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,6 +24,7 @@ subroutine sfcstr_cesm(m,n,mm,nn,k1m,k1n)
c --- ------------------------------------------------------------------
c
use mod_xc
+ use mod_constants, only: P_mks2cgs
use mod_forcing, only: ztx, mty, taux, tauy
use mod_checksum, only: csdiag, chksummsk
c
@@ -37,12 +38,12 @@ subroutine sfcstr_cesm(m,n,mm,nn,k1m,k1n)
do j=1,jj
do l=1,isu(j)
do i=max(1,ifu(j,l)),min(ii,ilu(j,l))
- taux(i,j)=10.*ztx(i,j)
+ taux(i,j)=P_mks2cgs*ztx(i,j)
enddo
enddo
do l=1,isv(j)
do i=max(1,ifv(j,l)),min(ii,ilv(j,l))
- tauy(i,j)=10.*mty(i,j)
+ tauy(i,j)=P_mks2cgs*mty(i,j)
enddo
enddo
enddo
diff --git a/cesm/thermf_cesm.F b/cesm/thermf_cesm.F
index 63b6a4e8..9b9740a0 100644
--- a/cesm/thermf_cesm.F
+++ b/cesm/thermf_cesm.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2008-2021 Mats Bentsen
+! Copyright (C) 2008-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -21,7 +21,8 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
c
c --- NERSC version of thermf. To be used when coupled to CESM
c
- use mod_constants, only: g, spcifh, t0deg, epsil, onem
+ use mod_constants, only: g, spcifh, t0deg, alpha0, epsilt, onem,
+ . g2kg, kg2g, L_mks2cgs, M_mks2cgs
use mod_time, only: nstep, nstep_in_day, nday_in_year,
. nday_of_year, baclin,
. xmi, l1mi, l2mi, l3mi, l4mi, l5mi
@@ -61,7 +62,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
integer i,j,k,l,m1,m2,m3,m4,m5
real y,dpotl,hotl,totl,sotl,dpmxl,hmxl,tmxl,smxl,tice_f,fwflx,
. sstc,rice,trxflx,sssc,srxflx,totsfl,totwfl,sflxc,totsrp,
- . totsrn,qp,qn
+ . totsrn,qp,qn,A_cgs2mks
#ifdef TRC
integer nt
real, dimension(ntr,1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
@@ -72,6 +73,8 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
real intp1d
external intp1d
c
+ A_cgs2mks=1./(L_mks2cgs**2)
+c
c --- Set parameters for time interpolation when applying diagnosed heat
c --- and salt relaxation fluxes
y=(nday_of_year-1+mod(nstep,nstep_in_day)/real(nstep_in_day))*48.
@@ -132,10 +135,10 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
c
c --- --- Salt flux due to brine rejection of freezing sea
c --- --- ice [kg m-2 m-1] (positive downwards)
- brnflx(i,j)=max(0.,-sotl*fmltfz(i,j)*1.e-3+sfl(i,j))
+ brnflx(i,j)=max(0.,-sotl*fmltfz(i,j)*g2kg+sfl(i,j))
c
c --- --- Virtual salt flux [kg m-2 s-1] (positive downwards)
- vrtsfl(i,j)=-sotl*fwflx*1.e-3
+ vrtsfl(i,j)=-sotl*fwflx*g2kg
c
c --- --- Store area weighted virtual salt flux and fresh water flux
util1(i,j)=vrtsfl(i,j)*scp2(i,j)
@@ -150,20 +153,21 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
c --- --- be heated. Note the freezing potential is multiplied by 1/2
c --- --- due to the leap-frog time stepping. The melting potential uses
c --- --- time averaged quantities since it is not accumulated.
- frzpot(i,j)=max(0.,tice_f-totl)*spcifh*dpotl/(2.*g)*1.e4
+ frzpot(i,j)=max(0.,tice_f-totl)*spcifh*dpotl
+ . /(2.*g)*(L_mks2cgs**2)
mltpot(i,j)=
. min(0.,tfrzm(i,j)-.5*(temp(i,j,k1m)+temp(i,j,k1n)))
- . *spcifh*.5*(dp(i,j,k1m)+dp(i,j,k1n))/g*1.e4
+ . *spcifh*.5*(dp(i,j,k1m)+dp(i,j,k1n))/g*(L_mks2cgs**2)
c
c --- --- Heat flux due to melting/freezing [W m-2] (positive downwards)
hmltfz(i,j)=hmlt(i,j)+frzpot(i,j)/baclin
c
c --- --- Total heat flux in BLOM units [W cm-2] (positive upwards)
- surflx(i,j)=-(swa(i,j)+nsf(i,j)+hmltfz(i,j))*1.e-4
+ surflx(i,j)=-(swa(i,j)+nsf(i,j)+hmltfz(i,j))*A_cgs2mks
c
c --- --- Short-wave heat flux in BLOM units [W cm-2] (positive
c --- --- upwards)
- sswflx(i,j)=-swa(i,j)*1.e-4
+ sswflx(i,j)=-swa(i,j)*A_cgs2mks
c
#ifdef TRC
c --- ------------------------------------------------------------------
@@ -182,7 +186,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
if (nt.eq.itrgls) then
trflx(nt,i,j)=-gls_n*difdia(i,j,1)*(gls_cmu0**gls_p)
. *(trc(i,j,k1n,itrtke)**gls_m)
- . *(vonKar**gls_n)*Zos**(gls_n-1.)
+ . *(vonKar**gls_n)*zos**(gls_n-1.)
ttrsf(nt,i,j)=0.
ttrav(nt,i,j)=0.
cycle
@@ -196,11 +200,12 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
endif
# endif
# endif
- trflx(nt,i,j)=-trc(i,j,k1n,nt)*fwflx*1.e-3
+ trflx(nt,i,j)=-trc(i,j,k1n,nt)*fwflx*g2kg
ttrsf(nt,i,j)=trflx(nt,i,j)*scp2(i,j)
ttrav(nt,i,j)=trc(i,j,k1n,nt)*scp2(i,j)
enddo
#endif
+c
c --- ------------------------------------------------------------------
c --- --- Relaxation fluxes
c --- ------------------------------------------------------------------
@@ -208,7 +213,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
surrlx(i,j)=0.
c
c --- --- If trxday>0 , apply relaxation towards observed sst
- if (trxday.gt.epsil) then
+ if (trxday.gt.epsilt) then
sstc=intp1d(sstclm(i,j,l1mi),sstclm(i,j,l2mi),
. sstclm(i,j,l3mi),sstclm(i,j,l4mi),
. sstclm(i,j,l5mi),xmi)
@@ -216,8 +221,8 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
. ricclm(i,j,l3mi),ricclm(i,j,l4mi),
. ricclm(i,j,l5mi),xmi)
sstc=(1.-rice)*max(sstc,tice_f)+rice*tice_f
- trxflx=spcifh*100.*min(hmxl,trxdpt)/(trxday*86400.)
- . *min(trxlim,max(-trxlim,sstc-tmxl))
+ trxflx=spcifh*L_mks2cgs*min(hmxl,trxdpt)/(trxday*86400.)
+ . *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
surrlx(i,j)=-trxflx
else
trxflx=0.
@@ -239,12 +244,12 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
salrlx(i,j)=0.
c
c --- --- if srxday>0 , apply relaxation towards observed sss
- if (srxday.gt.epsil) then
+ if (srxday.gt.epsilt) then
sssc=intp1d(sssclm(i,j,l1mi),sssclm(i,j,l2mi),
. sssclm(i,j,l3mi),sssclm(i,j,l4mi),
. sssclm(i,j,l5mi),xmi)
- srxflx=100.*min(hmxl,srxdpt)/(srxday*86400.)
- . *min(srxlim,max(-srxlim,sssc-smxl))
+ srxflx=L_mks2cgs*min(hmxl,srxdpt)/(srxday*86400.)
+ . *min(srxlim,max(-srxlim,sssc-smxl))/alpha0
salrlx(i,j)=-srxflx
util3(i,j)=max(0.,salrlx(i,j))*scp2(i,j)
util4(i,j)=min(0.,salrlx(i,j))*scp2(i,j)
@@ -269,7 +274,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
c --- --- Friction velocity (cm/s)
c --- -------------------------------------------------------------------
c
- ustar(i,j)=ustarw(i,j)*1.e2
+ ustar(i,j)=ustarw(i,j)*L_mks2cgs
c
enddo
enddo
@@ -287,7 +292,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
call xcsum(totwfl,util2,ips)
c
c --- Correction for the virtual salt flux [kg m-2 s-1]
- sflxc=(-sref*totwfl*1.e-3-totsfl)/area
+ sflxc=(-sref*totwfl*g2kg-totsfl)/area
if (mnproc.eq.1) then
write (lp,*) 'thermf: totsfl/area,sflxc',totsfl/area,sflxc
endif
@@ -298,8 +303,10 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
do j=1,jj
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- salflx(i,j)=-(vrtsfl(i,j)+sflxc+sfl(i,j))*1.e2
- brnflx(i,j)=-brnflx(i,j)*1.e2
+ salflx(i,j)=-(vrtsfl(i,j)+sflxc+sfl(i,j))
+ . *(kg2g*(M_mks2cgs/L_mks2cgs**2))
+ brnflx(i,j)=-brnflx(i,j)
+ . *(kg2g*(M_mks2cgs/L_mks2cgs**2))
enddo
enddo
enddo
@@ -308,7 +315,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
c --- if srxday>0 and srxbal=.true. , balance the sss relaxation flux
c --- so the net input of salt in grid cells connected to the world
c --- ocean is zero
- if (srxday.gt.epsil.and.srxbal) then
+ if (srxday.gt.epsilt.and.srxbal) then
call xcsum(totsrp,util3,ipwocn)
call xcsum(totsrn,util4,ipwocn)
if (abs(totsrp-totsrn).gt.0.) then
@@ -352,14 +359,14 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
tottrav=tottrav/area
c
trflxc=(-tottrsf)/area
-c trflxc=(-tottrav*totwfl*1.e-3-tottrsf)/area
+c trflxc=(-tottrav*totwfl*g2kg-tottrsf)/area
c trflxc=0.
c
c$OMP PARALLEL DO PRIVATE(l,i)
do j=1,jj
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- trflx(nt,i,j)=-(trflx(nt,i,j)+trflxc)*1.e2
+ trflx(nt,i,j)=-(trflx(nt,i,j)+trflxc)*L_mks2cgs
enddo
enddo
enddo
diff --git a/channel/thermf_channel.F b/channel/thermf_channel.F
index bc0ee447..41a565f9 100644
--- a/channel/thermf_channel.F
+++ b/channel/thermf_channel.F
@@ -24,7 +24,7 @@ subroutine thermf_channel(m,n,mm,nn,k1m,k1n)
use mod_xc
use mod_types, only: r8
use mod_ben02, only: ntda
- use mod_constants, only: spcifh, t0deg, epsil, onem
+ use mod_constants, only: spcifh, t0deg, epsilt, onem
use mod_time, only: nday_in_year, nday_of_year, nstep,
. nstep_in_day, baclin,
. xmi, l1mi, l2mi, l3mi, l4mi, l5mi
@@ -217,7 +217,7 @@ subroutine thermf_channel(m,n,mm,nn,k1m,k1n)
surrlx(i,j)=0._r8
!
! --- --- If trxday>0 , apply relaxation towards observed sst
- if (trxday.gt.epsil) then
+ if (trxday.gt.epsilt) then
sstc=intp1d(sstclm(i,j,l1mi),sstclm(i,j,l2mi),
. sstclm(i,j,l3mi),sstclm(i,j,l4mi),
. sstclm(i,j,l5mi),xmi)
@@ -248,7 +248,7 @@ subroutine thermf_channel(m,n,mm,nn,k1m,k1n)
salrlx(i,j)=0._r8
!
! --- --- if srxday>0 , apply relaxation towards observed sss
- if (srxday.gt.epsil) then
+ if (srxday.gt.epsilt) then
sssc=intp1d(sssclm(i,j,l1mi),sssclm(i,j,l2mi),
. sssclm(i,j,l3mi),sssclm(i,j,l4mi),
. sssclm(i,j,l5mi),xmi)
@@ -319,7 +319,7 @@ subroutine thermf_channel(m,n,mm,nn,k1m,k1n)
! --- if srxday>0 and srxbal=.true. , balance the sss relaxation flux
! --- so the net input of salt in grid cells connected to the world
! --- ocean is zero
- if (srxday.gt.epsil.and.srxbal) then
+ if (srxday.gt.epsilt.and.srxbal) then
call xcsum(totsrp,util3,ipwocn)
call xcsum(totsrn,util4,ipwocn)
if (abs(totsrp).gt.abs(totsrn)) then
diff --git a/cime_config/buildcpp b/cime_config/buildcpp
index a5765864..ca879e86 100644
--- a/cime_config/buildcpp
+++ b/cime_config/buildcpp
@@ -88,6 +88,7 @@ def buildcpp(case):
hamocc_sedbypass = case.get_value("HAMOCC_SEDBYPASS")
hamocc_ciso = case.get_value("HAMOCC_CISO")
hamocc_vsls = case.get_value("HAMOCC_VSLS")
+ blom_unit = case.get_value("BLOM_UNIT")
expect(blom_vcoord != "cntiso_hybrid" or not turbclo, "BLOM_VCOORD == {} and BLOM_TURBULENT_CLOSURE == {} is not a valid combination".format(blom_vcoord, turbclo))
@@ -149,6 +150,11 @@ def buildcpp(case):
else:
expect(False, "tracer module {} is not recognized".format(module))
+ if blom_unit == "mks":
+ blom_cppdefs = blom_cppdefs + " -DMKS"
+ else:
+ expect(blom_unit == "cgs", "Unit system {} is not recognized".format(blom_unit))
+
blom_cppdefs = "-DMPI" + blom_cppdefs
# update the xml variable BLOM_CPPDEFS with the above definition
diff --git a/cime_config/buildnml b/cime_config/buildnml
index 57cc45eb..56c2ad5a 100755
--- a/cime_config/buildnml
+++ b/cime_config/buildnml
@@ -7,6 +7,7 @@
set CASEROOT = `./xmlquery CASEROOT --value`
set OCN_GRID = `./xmlquery OCN_GRID --value`
set BLOM_VCOORD = `./xmlquery BLOM_VCOORD --value`
+set BLOM_UNIT = `./xmlquery BLOM_UNIT --value`
set DIN_LOC_ROOT = `./xmlquery DIN_LOC_ROOT --value`
set RUN_TYPE = `./xmlquery RUN_TYPE --value`
set CONTINUE_RUN = `./xmlquery CONTINUE_RUN --value`
@@ -74,20 +75,37 @@ set EXPCNF = "'cesm'"
set RUNTYP = "'$RUN_TYPE'"
set GRFILE = "'unset'"
set ICFILE = "'unset'"
-set PREF = 2000.e5
+if ($BLOM_UNIT == cgs) then
+ set PREF = 2000.e5
+else
+ set PREF = 2000.e4
+endif
set BACLIN = 1800.
set BATROP = 36.
-set MDV2HI = 2.
-set MDV2LO = .4
-set MDV4HI = 0.
-set MDV4LO = 0.
-set MDC2HI = 5000.e4
-set MDC2LO = 300.e4
+if ($BLOM_UNIT == cgs) then
+ set MDV2HI = 2.
+ set MDV2LO = .4
+ set MDV4HI = 0.
+ set MDV4LO = 0.
+ set MDC2HI = 5000.e4
+ set MDC2LO = 300.e4
+else
+ set MDV2HI = .02
+ set MDV2LO = .004
+ set MDV4HI = 0.
+ set MDV4LO = 0.
+ set MDC2HI = 5000.
+ set MDC2LO = 300.
+endif
set VSC2HI = .5
set VSC2LO = .5
set VSC4HI = 0.
set VSC4LO = 0.
-set CBAR = 5.
+if ($BLOM_UNIT == cgs) then
+ set CBAR = 5.
+else
+ set CBAR = .05
+endif
set CB = .002
set CWBDTS = 5.e-5
set CWBDLS = 25.
@@ -161,14 +179,25 @@ set EDWMTH = "'smooth'"
set EDSPRS = .true.
set EGC = 0.85
set EGGAM = 200.
-set EGLSMN = 4000.e2
-set EGMNDF = 100.e4
-set EGMXDF = 1500.e4
+if ($BLOM_UNIT == cgs) then
+ set EGLSMN = 4000.e2
+ set EGMNDF = 100.e4
+ set EGMXDF = 1500.e4
+else
+ set EGLSMN = 4000.
+ set EGMNDF = 100.
+ set EGMXDF = 1500.
+endif
set EGIDFQ = 1.
set RI0 = 1.2
set BDMTYP = 2
-set BDMC1 = 5.e-4
-set BDMC2 = .1
+if ($BLOM_UNIT == cgs) then
+ set BDMC1 = 5.e-4
+ set BDMC2 = .1
+else
+ set BDMC1 = 5.e-8
+ set BDMC2 = 1.e-5
+endif
set TKEPF = .006
if ($BLOM_VCOORD == isopyc_bulkml) then
set LTEDTP = "'layer'"
@@ -624,7 +653,11 @@ else if ($OCN_GRID == tnx2v1 ) then
set BACLIN = 4800.
set BATROP = 96.
set EGC = 0.5
- set EGMXDF = 1000.e4
+ if ($BLOM_UNIT == cgs) then
+ set EGMXDF = 1000.e4
+ else
+ set EGMXDF = 1000.
+ endif
set CWMTAG = "'Gibraltar','Gibraltar'"
set CWMEDG = " 'u', 'u'"
set CWMI = " 53, 54"
@@ -634,7 +667,11 @@ else if ($OCN_GRID == tnx1.5v1 ) then
set BACLIN = 4800.
set BATROP = 96.
set EGC = 0.5
- set EGMXDF = 1000.e4
+ if ($BLOM_UNIT == cgs) then
+ set EGMXDF = 1000.e4
+ else
+ set EGMXDF = 1000.
+ endif
else if ($OCN_GRID == tnx1v1 || $OCN_GRID == tnx1v3 || $OCN_GRID == tnx1v4) then
if ($OCN_NCPL == 24) then
set BACLIN = 3600.
@@ -655,33 +692,58 @@ else if ($OCN_GRID == tnx1v1 || $OCN_GRID == tnx1v3 || $OCN_GRID == tnx1v4) then
else if ($OCN_GRID == tnx0.25v1 || $OCN_GRID == tnx0.25v3 || $OCN_GRID == tnx0.25v4) then
set BACLIN = 900.
set BATROP = 15.
- set MDV2HI = .15
- set MDV2LO = .15
+ if ($BLOM_UNIT == cgs) then
+ set MDV2HI = .15
+ set MDV2LO = .15
+ set MDV4HI = 0.
+ set MDV4LO = 0.
+ set MDC2HI = 300.e4
+ set MDC2LO = 300.e4
+ else
+ set MDV2HI = .0015
+ set MDV2LO = .0015
+ set MDV4HI = 0.
+ set MDV4LO = 0.
+ set MDC2HI = 300.
+ set MDC2LO = 300.
+ endif
set VSC2HI = .15
set VSC2LO = .15
set VSC4HI = 0.0625
set VSC4LO = 0.0625
- set MDC2HI = 300.e4
set CWBDTS = 0.75e-4
set CWBDLS = 25.
set EDWMTH = "'step'"
set EGC = 0.85
- set EGMXDF = 1500.e4
+ if ($BLOM_UNIT == cgs) then
+ set EGMXDF = 1500.e4
+ else
+ set EGMXDF = 1500.
+ endif
set CE = 1.0
else if ($OCN_GRID == tnx0.125v4) then
set BACLIN = 300.
set BATROP = 6.
- set EGMNDF = 0.0
- set EGMXDF = 0.0
+ set EGMNDF = 0.
+ set EGMXDF = 0.
set EDWMTH = "'step'"
set CWBDTS = .75e-4
set CWBDLS = 25
- set MDV2HI = .5
- set MDV2LO = .1
- set MDV4HI = 0.
- set MDV4LO = 0.
- set MDC2HI = 300.e4
- set MDC2LO = 100.e4
+ if ($BLOM_UNIT == cgs) then
+ set MDV2HI = .5
+ set MDV2LO = .1
+ set MDV4HI = 0.
+ set MDV4LO = 0.
+ set MDC2HI = 300.e4
+ set MDC2LO = 100.e4
+ else
+ set MDV2HI = .005
+ set MDV2LO = .001
+ set MDV4HI = 0.
+ set MDV4LO = 0.
+ set MDC2HI = 300.
+ set MDC2LO = 100.
+ endif
set VSC2HI = .5
set VSC2LO = .5
set VSC4HI = 0.0
diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml
index 84973cb0..8ed52757 100644
--- a/cime_config/config_component.xml
+++ b/cime_config/config_component.xml
@@ -227,6 +227,15 @@
Optional turbulent closure. Valid values one of: twoeq oneeq. Additional values: advection isodif
+
+ char
+
+ cgs
+ build_component_blom
+ env_build.xml
+ Unit system. Valid values one of: cgs mks.
+
+
BLOM default:
BLOM/Ecosystem:
diff --git a/drivers/mct/domain_mct.F b/drivers/mct/domain_mct.F
index 34e07813..4c088b69 100644
--- a/drivers/mct/domain_mct.F
+++ b/drivers/mct/domain_mct.F
@@ -27,6 +27,7 @@ subroutine domain_mct(gsMap_ocn, dom_ocn, lsize, perm, jjcpl)
use mod_types, only: r8
use mod_xc
use mod_grid, only: scp2, plon, plat
+ use mod_constants, only: L_mks2cgs
implicit none
@@ -105,7 +106,7 @@ subroutine domain_mct(gsMap_ocn, dom_ocn, lsize, perm, jjcpl)
enddo
call mct_gGrid_importRattr(dom_ocn, "lat", rdata, lsize)
- radius = SHR_CONST_REARTH*1.e2_r8 ! Earth's radius in cm
+ radius = SHR_CONST_REARTH*L_mks2cgs ! Earth's radius in cm
n = 0
do j = 1, jjcpl
diff --git a/drivers/mct/export_mct.F b/drivers/mct/export_mct.F
index d5ca668e..cff6c4c2 100644
--- a/drivers/mct/export_mct.F
+++ b/drivers/mct/export_mct.F
@@ -23,6 +23,7 @@ subroutine export_mct(o2x_o, lsize, perm, jjcpl, nsend, sbuff,
! Uses modules
use mct_mod
+ use mod_constants, only: L_mks2cgs
use shr_const_mod, only: SHR_CONST_TKFRZ
use mod_types, only: r8
use blom_cpl_indices
@@ -47,8 +48,10 @@ subroutine export_mct(o2x_o, lsize, perm, jjcpl, nsend, sbuff,
integer i, j, n
real(r8) :: tfac, utmp, vtmp
+ real(r8) :: iL_mks2cgs
tfac = 1._r8/tlast_coupled
+ iL_mks2cgs = 1._r8/L_mks2cgs
! ----------------------------------------------------------------
! Interpolate onto scalar points, rotate, and pack surface
@@ -73,9 +76,9 @@ subroutine export_mct(o2x_o, lsize, perm, jjcpl, nsend, sbuff,
vtmp = .5_r8*( sbuff(i,j ,index_o2x_So_v)
. + sbuff(i,j+1,index_o2x_So_v))
o2x_o%rattr(index_o2x_So_u,n) =
- . (utmp*cosang(i,j) - vtmp*sinang(i,j))*tfac*1.e-2_r8
+ . (utmp*cosang(i,j) - vtmp*sinang(i,j))*tfac*iL_mks2cgs
o2x_o%rattr(index_o2x_So_v,n) =
- . (utmp*sinang(i,j) + vtmp*cosang(i,j))*tfac*1.e-2_r8
+ . (utmp*sinang(i,j) + vtmp*cosang(i,j))*tfac*iL_mks2cgs
utmp = ( sbuff(i ,j,index_o2x_So_dhdx)*iu(i ,j)
. + sbuff(i+1,j,index_o2x_So_dhdx)*iu(i+1,j))
. /max(1,iu(i,j) + iu(i+1,j))
diff --git a/fuk95/mod_fuk95.F90 b/fuk95/mod_fuk95.F90
index 6dd0a7bc..9a2b7054 100644
--- a/fuk95/mod_fuk95.F90
+++ b/fuk95/mod_fuk95.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2021 Mats Bentsen
+! Copyright (C) 2021-2022 Mats Bentsen
!
! This file is part of BLOM.
!
@@ -25,7 +25,8 @@ module mod_fuk95
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: g, rearth, pi, radian, epsil
+ use mod_constants, only: g, rearth, rho0, pi, radian, epsilz, &
+ L_mks2cgs, R_mks2cgs
use mod_xc
use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml, cntiso_hybrid, sigmar
use mod_grid, only: qclon, qclat, pclon, pclat, uclon, uclat, vclon, vclat, &
@@ -45,19 +46,32 @@ module mod_fuk95
private
+! real(r8), parameter :: &
+! u0 = 30._r8, & ! Maximum jet velocity [cm s-1].
+! h1 = 1.e4_r8, & ! Depth of active layer [cm].
+! h0 = 2.e4_r8, & ! Depth of water column [cm].
+! l0 = 2.e6_r8, & ! Half-width of the jet [cm].
+! drho = 0.19e-3_r8, & ! Active layer density difference [g cm-3].
+! rhoc = 1.0259_r8, & ! Density at the center of active layer [g cm-3].
+! rhob = 1.0270_r8, & ! Density of water beneath active layer [g cm-3].
+! f = 1.e-4_r8, & ! Coriolis parameter [1 s-1].
+! lat0 = 45._r8, & ! Center latitude of grid domain [deg].
+! lambda = 20.8e5, & ! Channel length [cm].
+! mindz = 1.e2_r8, & ! Minimum interior layer thickness [cm].
+! saln0 = 35._r8 ! Constant salinity value [g kg-1].
real(r8), parameter :: &
- u0 = 30._r8, & ! Maximum jet velocity [cm s-1].
- h1 = 1.e4_r8, & ! Depth of active layer [cm].
- h0 = 2.e4_r8, & ! Depth of water column [cm].
- l0 = 2.e6_r8, & ! Half-width of the jet [cm].
- drho = 0.19e-3_r8, & ! Active layer density difference [g cm-3].
- rho1 = 1.0259_r8, & ! Density at the center of active layer [g cm-3].
- rho0 = 1.0270_r8, & ! Density of water beneath active layer [g cm-3].
- f = 1.e-4_r8, & ! Coriolis parameter [1 s-1].
- lat0 = 45._r8, & ! Center latitude of grid domain [deg].
- lambda = 20.8e5, & ! Channel length [cm].
- mindz = 1.e2_r8, & ! Minimum interior layer thickness [cm].
- saln0 = 35._r8 ! Constant salinity value [g kg-1].
+ u0 = .3_r8*L_mks2cgs, & ! Maximum jet velocity [m s-1].
+ h1 = 1.e2_r8*L_mks2cgs, & ! Depth of active layer [m].
+ h0 = 2.e2_r8*L_mks2cgs, & ! Depth of water column [m].
+ l0 = 2.e4_r8*L_mks2cgs, & ! Half-width of the jet [m].
+ drho = 0.19_r8*R_mks2cgs, & ! Active layer density difference [kg m-3].
+ rhoc = 1025.9_r8*R_mks2cgs, & ! Density at the center of active layer [kg m-3].
+ rhob = 1027.0_r8*R_mks2cgs, & ! Density of water beneath active layer [kg m-3].
+ f = 1.e-4_r8, & ! Coriolis parameter [1 s-1].
+ lat0 = 45._r8, & ! Center latitude of grid domain [deg].
+ lambda = 20.8e3*L_mks2cgs, & ! Channel length [m].
+ mindz = 1._r8*L_mks2cgs, & ! Minimum interior layer thickness [m].
+ saln0 = 35._r8 ! Constant salinity value [g kg-1].
public :: geoenv_fuk95, inifrc_fuk95, ictsz_fuk95
@@ -132,7 +146,7 @@ subroutine geoenv_fuk95
tmpg(1 , j) = 0._r8
tmpg(itdm, j) = 0._r8
do i = 2, itdm - 1
- tmpg(i, j) = h0*1.e-2
+ tmpg(i, j) = h0*L_mks2cgs**(-1)
enddo
enddo
!$omp end parallel do
@@ -281,10 +295,10 @@ subroutine ictsz_fuk95
! and corresponding isopycnic layer structure. The bulk mixed layer
! is set to the minimum mixed layer thickness.
- drhojet = rho1*f*u0*l0/(g*h1)
+ drhojet = rhoc*f*u0*l0/(g*h1)
dsig = (drho + drhojet)/(kk - 4)
- sigref(kk) = rho0 - 1._r8
- sigref(kk - 1) = rho1 + .5_r8*(drho + drhojet) - 1._r8
+ sigref(kk) = rhob - rho0
+ sigref(kk - 1) = rhoc + .5_r8*(drho + drhojet) - rho0
do k = kk - 2, 1, -1
sigref(k) = sigref(k + 1) - dsig
enddo
@@ -310,11 +324,11 @@ subroutine ictsz_fuk95
do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
x = x_nudge(real(i, r8), real(j, r8))
z(i, j, 1) = 0._r8
- z(i, j, 2) = .5_r8*mltmin*1.e2_r8
- z(i, j, 3) = mltmin*1.e2_r8
+ z(i, j, 2) = .5_r8*mltmin*L_mks2cgs
+ z(i, j, 3) = mltmin*L_mks2cgs
z(i, j, kk ) = h1
z(i, j, kk + 1) = h0
- sigm = rho1*(1._r8 + f*u0*x_psi(x)/(g*h1)) - 1._r8
+ sigm = rhoc*(1._r8 + f*u0*x_psi(x)/(g*h1)) - rho0
sigma(i, j, 1) = sigm &
+ .5_r8*drho*(z(i, j, 2) + z(i, j, 1) - h1)/h1
sigma(i, j, 2) = sigm &
@@ -327,7 +341,7 @@ subroutine ictsz_fuk95
do l = 1, isp(j)
do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
x = x_nudge(real(i, r8), real(j, r8))
- sigm = rho1*(1._r8 + f*u0*x_psi(x)/(g*h1)) - 1._r8
+ sigm = rhoc*(1._r8 + f*u0*x_psi(x)/(g*h1)) - rho0
sigi = .5_r8*(sigref(k - 1) + sigref(k))
z(i, j, k) = ((sigi - sigm)/drho + .5_r8)*h1
z(i, j, k) = min(z(i, j, kk) - mindz*(kk - k), &
@@ -347,20 +361,20 @@ subroutine ictsz_fuk95
! active layer is distributed equally among the remaining model
! layers using constant z-level interfaces.
-! drhojet = rho1*f*u0*l0/(g*h1)
+! drhojet = rhoc*f*u0*l0/(g*h1)
! dsig = (drho + drhojet)/(kk - 4)
-! sigref(kk) = .5_r8*(rho0 + rho1) + .25_r8*(drho + drhojet) - 1._r8
-! sigref(kk - 1) = rho1 + .5_r8*(drho + drhojet - dsig) - 1._r8
+! sigref(kk) = .5_r8*(rhob + rhoc) + .25_r8*(drho + drhojet) - rho0
+! sigref(kk - 1) = rhoc + .5_r8*(drho + drhojet - dsig) - rho0
! do k = kk - 2, 1, -1
! sigref(k) = sigref(k + 1) - dsig
! enddo
- drhojet = rho1*f*u0*l0/(g*h1)
+ drhojet = rhoc*f*u0*l0/(g*h1)
dsig = (drho + drhojet)/(kk - 5)
- sigref(kk - 2) = rho1 + .5_r8*(drho + drhojet - dsig) - 1._r8
+ sigref(kk - 2) = rhoc + .5_r8*(drho + drhojet - dsig) - rho0
do k = kk - 3, 1, -1
sigref(k) = sigref(k + 1) - dsig
enddo
- sigref(kk ) = rho0 - 1._r8
+ sigref(kk ) = rhob - rho0
sigref(kk - 1) = (2._r8*sigref(kk - 2) + sigref(kk))/3._r8
sigref(kk ) = (sigref(kk - 2) + 2._r8*sigref(kk))/3._r8
@@ -383,14 +397,14 @@ subroutine ictsz_fuk95
enddo
!$omp end parallel do
- s0 = rho0 - 1._r8
+ s0 = rhob - rho0
!$omp parallel do private(k, l, i, x, s1)
do j = 1, jj
do k = 1, kk
do l = 1, isp(j)
do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
x = x_nudge(real(i, r8), real(j, r8))
- s1 = rho1*(1._r8 + f*u0*x_psi(x)/(g*h1)) - 1._r8 &
+ s1 = rhoc*(1._r8 + f*u0*x_psi(x)/(g*h1)) - rho0 &
+ .5_r8*drho*(z(i, j, k + 1) + z(i, j, k) - h1)/h1
sigma(i, j, k) = &
( s1*max(0._r8, min(z(i, j, k + 1), h1) - z(i, j, k)) &
@@ -426,7 +440,7 @@ subroutine ictsz_fuk95
zl = .5_r8*(z(i, j - 1, k + 1) + z(i, j, k + 1))
v1 = u0*psi(x)*(h1 - .5*(zu + zl))/h1
v1 = 0._r8
- if (abs(zl - zu) < epsil) then
+ if (abs(zl - zu) < epsilz) then
v(i, j, k) = v1
else
v(i, j, k) = ( v1*max(0._r8, min(zl, h1) - zu) &
diff --git a/meson.build b/meson.build
index baf58d7d..9f237b7e 100644
--- a/meson.build
+++ b/meson.build
@@ -71,6 +71,10 @@ subdir('pkgs/')
# Handle options and add necessary flags and subfolders with source files
+if get_option('mks')
+ add_project_arguments('-DMKS', language: 'fortran')
+endif
+
turbclo = get_option('turbclo')
if turbclo.length() > 0 and get_option('vcoord') == 'cntiso_hybrid'
message('Setting turbclo = [] for vcoord == \'cntiso_hybrid\'')
diff --git a/meson_options.txt b/meson_options.txt
index af50852a..8e48383f 100644
--- a/meson_options.txt
+++ b/meson_options.txt
@@ -13,6 +13,8 @@ option('vcoord', type: 'combo',
option('driver', type: 'combo',
choices: ['nocoupler', 'noforc'], value: 'nocoupler')
# List of BLOM options
+option('mks', type: 'boolean',
+ description: 'Enable MKS units', value: false)
option('turbclo', type: 'array',
choices: ['oneeq', 'twoeq', 'advection', 'isodif'],
description: 'Turbulent closure options', value: ['oneeq', 'advection'])
diff --git a/phy/convec.F b/phy/convec.F
index 83d12baa..9b68bfcb 100644
--- a/phy/convec.F
+++ b/phy/convec.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2009-2021 Mats Bentsen
+! Copyright (C) 2009-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,7 +24,7 @@ subroutine convec(m,n,mm,nn,k1m,k1n)
c --- layers
c --- ------------------------------------------------------------------
c
- use mod_constants, only: epsil
+ use mod_constants, only: epsilp
use mod_xc
use mod_vcoord, only: sigmar
use mod_eos, only: rho, sig, sofsig
@@ -84,7 +84,7 @@ subroutine convec(m,n,mm,nn,k1m,k1n)
c
k=3
dps=0.
- do while (delp(k).lt.epsil)
+ do while (delp(k).lt.epsilp)
dps=dps+delp(k)
delp(k)=0.
k=k+1
@@ -212,7 +212,7 @@ subroutine convec(m,n,mm,nn,k1m,k1n)
k=kfpl
do while (rho(dps,ttmp,stmp).gt.
. rho(dps,ttem(k),ssal(k)).or.
- . delp(k).lt.epsil)
+ . delp(k).lt.epsilp)
tdps=tdps+ttem(k)*delp(k)
sdps=sdps+ssal(k)*delp(k)
dps=dps+delp(k)
diff --git a/phy/diapfl.F b/phy/diapfl.F
index 0ba42b2b..6d6badc5 100644
--- a/phy/diapfl.F
+++ b/phy/diapfl.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2009-2021 Mats Bentsen, Mehmet Ilicak
+! Copyright (C) 2009-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -23,7 +23,7 @@ subroutine diapfl(m,n,mm,nn,k1m,k1n)
c --- Diapycnal mixing
c --- ------------------------------------------------------------------
c
- use mod_constants, only: g, alpha0, spval, epsil, onem
+ use mod_constants, only: g, alpha0, spval, epsilp, onem, L_mks2cgs
use mod_time, only: delt1
use mod_xc
use mod_vcoord, only: sigmar
@@ -64,7 +64,7 @@ subroutine diapfl(m,n,mm,nn,k1m,k1n)
c --- scale bottom boundary layer mixing [cm/s]
real dsgmnr,fcmxr,dsgcr0,dfeps,gbbl,kappa,ustmin
parameter (dsgmnr=.1,fcmxr=.25,dsgcr0=.25,dfeps=1.e-12,gbbl=.2,
- . kappa=.4,ustmin=.01)
+ . kappa=.4,ustmin=.0001*L_mks2cgs)
c
real, save, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm) ::
. fpug=spval,fplg=spval
@@ -133,7 +133,7 @@ subroutine diapfl(m,n,mm,nn,k1m,k1n)
kmin=kfpl-2
kmax=1
do k=2,kk
- if (delp(k).gt.epsil) kmax=k
+ if (delp(k).gt.epsilp) kmax=k
enddo
c
if (kmin.lt.kmax) then
@@ -314,7 +314,7 @@ subroutine diapfl(m,n,mm,nn,k1m,k1n)
. i+i0,j+j0
open (10,file='diapfl.uf',form='unformatted')
write (10) kk,kfpl
- write (10) g,alpha0,epsil,onem,delt1,dsgmnr,q,q
+ write (10) g,alpha0,epsilp,onem,delt1,dsgmnr,q,q
write (10) ttem0,ssal0,delp0,dens0,sigr0,nu0
close (10)
call xchalt('(diapfl)')
@@ -334,7 +334,7 @@ subroutine diapfl(m,n,mm,nn,k1m,k1n)
f(k)=min(fmax(k),
. .5*sqrt(c*nu(k)*dsg(k)
. *(dsgui(k)+dsgli(k)))*dsghm(k),
- . c*nu(k)*dsg(k)/max(epsil,delp(k)))
+ . c*nu(k)*dsg(k)/max(epsilp,delp(k)))
fold(k)=f(k)
h(k)=fcu(k )*dsgui(k )-fcl(k )*dsgli(k )
. +fcl(k-1)*dsgli(k-1)-fcu(k+1)*dsgui(k+1)
@@ -519,7 +519,7 @@ subroutine diapfl(m,n,mm,nn,k1m,k1n)
. i+i0,j+j0,maxdf,dflim
open (10,file='diapfl.uf',form='unformatted')
write (10) kk,kfpl
- write (10) g,alpha0,epsil,onem,delt1,dsgmnr,q,q
+ write (10) g,alpha0,epsilp,onem,delt1,dsgmnr,q,q
write (10) ttem0,ssal0,delp0,dens0,sigr0,nu0
close (10)
call xchalt('(diapfl)')
diff --git a/phy/diffus.F b/phy/diffus.F
index e99a4983..fc590fd7 100644
--- a/phy/diffus.F
+++ b/phy/diffus.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2006-2022 Mats Bentsen
+! Copyright (C) 2006-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,6 +24,7 @@ subroutine diffus(m,n,mm,nn,k1m,k1n)
c --- ------------------------------------------------------------------
c
use mod_time, only: delt1
+ use mod_constants, only: P_mks2cgs
use mod_xc
use mod_grid, only: scuy, scvx, scp2, scuxi, scvyi
use mod_eos, only: sig
@@ -47,7 +48,7 @@ subroutine diffus(m,n,mm,nn,k1m,k1n)
#endif
c
real dpeps
- parameter (dpeps=1.e-4)
+ parameter (dpeps=1.e-5*P_mks2cgs)
c
call xctilr(dp(1-nbdy,1-nbdy,k1n), 1,kk, 3,3, halo_ps)
if (ltedtp_opt.eq.ltedtp_neutral) then
diff --git a/phy/geoenv_file.F b/phy/geoenv_file.F
index 11ebf88c..c130fd50 100644
--- a/phy/geoenv_file.F
+++ b/phy/geoenv_file.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2020 Mats Bentsen, Ping-Gin Chiu
+! Copyright (C) 2015-2022 Mats Bentsen, Ping-Gin Chiu, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -25,7 +25,7 @@ subroutine geoenv_file
c --- ------------------------------------------------------------------
c
use mod_config, only: inst_suffix
- use mod_constants, only: rearth, pi, radian
+ use mod_constants, only: rearth, pi, radian, L_mks2cgs
use mod_xc
use mod_grid, only: grfile, qclon, qclat, pclon, pclat, uclon,
. uclat, vclon, vclat, scqx, scqy, scpx, scpy,
@@ -797,18 +797,18 @@ subroutine geoenv_file
do j=1,jj
do i=1,ii
c
- scqx(i,j)=scqx(i,j)*1.e2
- scqy(i,j)=scqy(i,j)*1.e2
- scpx(i,j)=scpx(i,j)*1.e2
- scpy(i,j)=scpy(i,j)*1.e2
- scux(i,j)=scux(i,j)*1.e2
- scuy(i,j)=scuy(i,j)*1.e2
- scvx(i,j)=scvx(i,j)*1.e2
- scvy(i,j)=scvy(i,j)*1.e2
- scq2(i,j)=scq2(i,j)*1.e4
- scp2(i,j)=scp2(i,j)*1.e4
- scu2(i,j)=scu2(i,j)*1.e4
- scv2(i,j)=scv2(i,j)*1.e4
+ scqx(i,j)=scqx(i,j)*L_mks2cgs
+ scqy(i,j)=scqy(i,j)*L_mks2cgs
+ scpx(i,j)=scpx(i,j)*L_mks2cgs
+ scpy(i,j)=scpy(i,j)*L_mks2cgs
+ scux(i,j)=scux(i,j)*L_mks2cgs
+ scuy(i,j)=scuy(i,j)*L_mks2cgs
+ scvx(i,j)=scvx(i,j)*L_mks2cgs
+ scvy(i,j)=scvy(i,j)*L_mks2cgs
+ scq2(i,j)=scq2(i,j)*L_mks2cgs**2
+ scp2(i,j)=scp2(i,j)*L_mks2cgs**2
+ scu2(i,j)=scu2(i,j)*L_mks2cgs**2
+ scv2(i,j)=scv2(i,j)*L_mks2cgs**2
c
cosang(i,j)=cos(angle(i,j))
sinang(i,j)=sin(angle(i,j))
diff --git a/phy/mod_cmnfld.F90 b/phy/mod_cmnfld.F90
index 6fe55510..4b9b8890 100644
--- a/phy/mod_cmnfld.F90
+++ b/phy/mod_cmnfld.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2022 Mats Bentsen
+! Copyright (C) 2015-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,7 +24,7 @@ module mod_cmnfld
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: spval
+ use mod_constants, only: spval, onem, L_mks2cgs
use mod_xc
implicit none
@@ -33,7 +33,7 @@ module mod_cmnfld
! Parameters:
real(r8) :: &
- sls0 = 10._r8*98060._r8, & ! Minimum smoothing length scale in the
+ sls0 = 10._r8*onem, & ! Minimum smoothing length scale in the
! computation of filtered BFSQ [g cm-1 s-2].
slsmfq = 2._r8, & ! Factor to be multiplied with the mixed
! layer depth to find the smoothing length
@@ -45,7 +45,7 @@ module mod_cmnfld
! computation of filtered BFSQ [].
bfsqmn = 1.e-7_r8, & ! Minimum value of BFSQ used in the
! computation of neutral slope [s-2].
- dbcrit = .03_r8 ! Critical buoyancy difference used in the
+ dbcrit = .0003_r8*L_mks2cgs! Critical buoyancy difference used in the
! mixed layer thickness estimation (Levitus,
! 1982) [cm s-2].
diff --git a/phy/mod_cmnfld_routines.F90 b/phy/mod_cmnfld_routines.F90
index 761ac198..ef6f7cb1 100644
--- a/phy/mod_cmnfld_routines.F90
+++ b/phy/mod_cmnfld_routines.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2022 Mats Bentsen
+! Copyright (C) 2015-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,7 +24,7 @@ module mod_cmnfld_routines
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: g, alpha0, epsil, onem, onecm, onemm
+ use mod_constants, only: g, alpha0, rho0, epsilp, onem, onecm, onemm
use mod_xc
use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml, cntiso_hybrid
use mod_grid, only: scuxi, scvyi
@@ -125,7 +125,7 @@ subroutine cmnfld_bfsqf_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
sup = saln(i, j, 2 + nn)
do k = kfpl, kk
kn = k + nn
- if (p(i, j, kk + 1) - p(i, j, k) < epsil) then
+ if (p(i, j, kk + 1) - p(i, j, k) < epsilp) then
delp(k) = onemm
bfsqi(i, j, k) = bfsqi(i, j, k - 1)
bfsq(k) = bfsqmn
@@ -133,7 +133,7 @@ subroutine cmnfld_bfsqf_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
q = max(sls0, delp(kfpl - 1)*slsmfq*q + sls0*(1._r8 - q))
sls2(k) = q*q
else
- if (p(i, j, kk + 1) - p(i, j, k + 1) < epsil) then
+ if (p(i, j, kk + 1) - p(i, j, k + 1) < epsilp) then
plo = p(i, j, kk + 1)
else
plo = .5_r8*(p(i, j, k) + p(i, j, k + 1))
@@ -262,13 +262,13 @@ subroutine cmnfld_bfsqf_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
sup = saln(i, j, 1 + nn)
do k = 2, kk
kn = k + nn
- if (p(i, j, kk + 1) - p(i, j, k) < epsil) then
+ if (p(i, j, kk + 1) - p(i, j, k) < epsilp) then
delp(k) = onemm
bfsqi(i, j, k) = bfsqi(i, j, k - 1)
bfsq(k) = bfsqmn
sls2(k) = sls0*sls0
else
- if (p(i, j, kk + 1) - p(i, j, k + 1) < epsil) then
+ if (p(i, j, kk + 1) - p(i, j, k + 1) < epsilp) then
plo = p(i, j, kk + 1)
else
plo = .5_r8*(p(i, j, k) + p(i, j, k + 1))
@@ -376,10 +376,10 @@ subroutine cmnfld_bfsqi_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
sup = saln(i, j, 1 + nn)
do k = 2, kk
kn = k + nn
- if (p(i, j, kk + 1) - p(i, j, k) < epsil) then
+ if (p(i, j, kk + 1) - p(i, j, k) < epsilp) then
bfsqi(i, j, k) = bfsqi(i, j, k - 1)
else
- if (p(i, j, kk + 1) - p(i, j, k + 1) < epsil) then
+ if (p(i, j, kk + 1) - p(i, j, k + 1) < epsilp) then
plo = p(i, j, kk + 1)
else
plo = .5_r8*(p(i, j, k) + p(i, j, k + 1))
@@ -420,7 +420,7 @@ subroutine cmnfld_nslope_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
integer, intent(in) :: m, n, mm, nn, k1m, k1n
- real(r8) :: rho0, pm, rho_x, phi_x, bfsqm, rho_y, phi_y
+ real(r8) :: pm, rho_x, phi_x, bfsqm, rho_y, phi_y
integer :: i, j, k, l, kn, kintr, kmax, knnsl
! ------------------------------------------------------------------------
@@ -433,7 +433,7 @@ subroutine cmnfld_nslope_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
kn = k + nn
do l = 1, isp(j)
do i = max(- 1, ifp(j, l)), min(ii + 2, ilp(j, l))
- if (dp(i, j, kn) < epsil) then
+ if (dp(i, j, kn) < epsilp) then
phi(i, j, k) = phi(i, j, k + 1)
else
phi(i, j, k) = phi(i, j, k + 1) &
@@ -454,8 +454,6 @@ subroutine cmnfld_nslope_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
! bathymetry and in this case values are extrapolated from above.
! ------------------------------------------------------------------------
- rho0 = 1._r8/alpha0
-
!$omp parallel do private(l, i, k, kmax, kn, kintr, knnsl, pm, rho_x, &
!$omp phi_x, bfsqm)
do j = - 1, jj + 2
@@ -475,7 +473,7 @@ subroutine cmnfld_nslope_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
kmax = 1
do k = 3, kk
kn = k + nn
- if (dp(i - 1, j, kn) > epsil .or. dp(i, j, kn) > epsil) &
+ if (dp(i - 1, j, kn) > epsilp .or. dp(i, j, kn) > epsilp) &
kmax = k
enddo
@@ -564,7 +562,7 @@ subroutine cmnfld_nslope_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
kmax = 1
do k = 3, kk
kn = k + nn
- if (dp(i, j - 1, kn) > epsil .or. dp(i, j, kn) > epsil) &
+ if (dp(i, j - 1, kn) > epsilp .or. dp(i, j, kn) > epsilp) &
kmax = k
enddo
@@ -653,7 +651,7 @@ subroutine cmnfld_nslope_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
integer, intent(in) :: m, n, mm, nn, k1m, k1n
- real(r8) :: rho0, pm, rho_x, phi_x, bfsqm, rho_y, phi_y
+ real(r8) :: pm, rho_x, phi_x, bfsqm, rho_y, phi_y
integer :: i, j, k, l, kn, kmax, knnsl
! ------------------------------------------------------------------------
@@ -666,7 +664,7 @@ subroutine cmnfld_nslope_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
kn = k + nn
do l = 1, isp(j)
do i = max(- 1, ifp(j, l)), min(ii + 2, ilp(j, l))
- if (dp(i, j, kn) < epsil) then
+ if (dp(i, j, kn) < epsilp) then
phi(i, j, k) = phi(i, j, k + 1)
else
phi(i, j, k) = phi(i, j, k + 1) &
@@ -687,8 +685,6 @@ subroutine cmnfld_nslope_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
! bathymetry and in this case values are extrapolated from above.
! ------------------------------------------------------------------------
- rho0 = 1._r8/alpha0
-
!$omp parallel do private(l, i, k, kmax, kn, knnsl, pm, rho_x, phi_x, bfsqm)
do j = - 1, jj + 2
do l = 1, isu(j)
@@ -705,7 +701,7 @@ subroutine cmnfld_nslope_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
kmax = 1
do k = 2, kk
kn = k + nn
- if (dp(i - 1, j, kn) > epsil .or. dp(i, j, kn) > epsil) kmax = k
+ if (dp(i - 1, j, kn) > epsilp .or. dp(i, j, kn) > epsilp) kmax=k
enddo
! Index of last interface where slope vector times Brunt-Vaisala
@@ -758,7 +754,7 @@ subroutine cmnfld_nslope_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
kmax = 1
do k = 2, kk
kn = k + nn
- if (dp(i, j - 1, kn) > epsil .or. dp(i, j, kn) > epsil) kmax = k
+ if (dp(i, j - 1, kn) > epsilp .or. dp(i, j, kn) > epsilp) kmax=k
enddo
! Index of last interface where slope vector times Brunt-Vaisala
@@ -881,7 +877,7 @@ subroutine cmnfld_z(m, n, mm, nn, k1m, k1n)
km = k + mm
do l = 1, isp(j)
do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
- if (dp(i, j, km) < epsil) then
+ if (dp(i, j, km) < epsilp) then
z(i, j, k) = z(i, j, k + 1)
else
z(i, j, k) = z(i, j, k + 1) &
@@ -934,7 +930,7 @@ subroutine cmnfld_mlts(m, n, mm, nn, k1m, k1n)
zup = zlo
dbup = dblo
else
- dbup = min(dbup, dbcrit - epsil)
+ dbup = min(dbup, dbcrit - epsilp)
mlts(i, j) = ( zup*(dblo - dbcrit) &
+ zlo*(dbcrit - dbup))/(dblo - dbup) &
- z(i, j, 1)
diff --git a/phy/mod_constants.F90 b/phy/mod_constants.F90
index 48c17b47..adc6ce9c 100644
--- a/phy/mod_constants.F90
+++ b/phy/mod_constants.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2020-2021 Mats Bentsen
+! Copyright (C) 2020-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -28,6 +28,39 @@ module mod_constants
private
+#ifdef MKS
+ ! MKS unit
+ real(r8), parameter :: &
+ g = 9.806_r8, & ! Gravitational acceleration [m s-2].
+ rearth = 6.37122e6_r8, & ! Radius of the Earth [m].
+ spcifh = 3990._r8, & ! Specific heat capacity of sea water
+ ! [J kg-1 K-1].
+ t0deg = 273.15_r8, & ! Zero degrees Celsius in Kelvin [K].
+ alpha0 = 1.e-3_r8, & ! Reference value of specific volume
+ ! [m3 kg-1].
+ rho0 = 1.e3_r8, & ! Reference value of density [kg m-3].
+ pi = 3.1415926536_r8, & ! pi [].
+ radian = 57.295779513_r8, & ! 180/pi [].
+ epsilpl = 1.e-14_r8, & ! Small value for pressure*dx [].
+ epsilp = 1.e-12_r8, & ! Small value for pressure [].
+ epsilz = 1.e-9_r8, & ! Small value for depth [].
+ epsilt = 1.e-11_r8, & ! Small value for time [].
+ epsilk = 1.e-15_r8, & ! Small value for kappa [].
+ spval = 1.e33_r8, & ! Large value [].
+ tenm = 98060._r8, & ! 10 m in units of pressure [kg m-1 s-2].
+ onem = 9806._r8, & ! 1 m in units of pressure [kg m-1 s-2].
+ tencm = 980.6_r8, & ! 10 cm in units of pressure [kg m-1 s-2].
+ onecm = 98.06_r8, & ! 1 cm in units of pressure [kg m-1 s-2].
+ onemm = 9.806_r8, & ! 1 mm in units of pressure [kg m-1 s-2].
+ onemu = .009806_r8, & ! 1 micrometer in units of pressure
+ ! [kg m-1 s-2].
+ g2kg = 1.e-3_r8, & ! convert g to kg coeff
+ kg2g = 1.e3_r8, & ! convert kg to g coeff
+ L_mks2cgs = 1._r8, & ! length coefficient converting CGS to MKS
+ M_mks2cgs = 1._r8, & ! mass coefficient converting CGS to MKS
+ P_mks2cgs = 1._r8, & ! pressure coefficient converting CGS to MKS
+ R_mks2cgs = 1._r8 ! rho coefficient converting CGS to MKS
+#else
real(r8), parameter :: &
g = 980.6_r8, & ! Gravitational acceleration [cm s-2].
rearth = 6.37122e8_r8, & ! Radius of the Earth [cm].
@@ -36,20 +69,33 @@ module mod_constants
t0deg = 273.15_r8, & ! Zero degrees Celsius in Kelvin [K].
alpha0 = 1._r8, & ! Reference value of specific volume
! [cm3 g-1].
+ rho0 = 1._r8, & ! Reference value of density [kg m-3].
pi = 3.1415926536_r8, & ! pi [].
radian = 57.295779513_r8, & ! 180/pi [].
- epsil = 1.e-11_r8, & ! Small value [].
+ epsilpl = 1.e-11_r8, & ! Small value for pressure*dx [].
+ epsilp = 1.e-11_r8, & ! Small value for pressure [].
+ epsilz = 1.e-11_r8, & ! Small value for depth [].
+ epsilt = 1.e-11_r8, & ! Small value for time [].
+ epsilk = 1.e-11_r8, & ! Small value for kappa [].
spval = 1.e33_r8, & ! Large value [].
tenm = 980600._r8, & ! 10 m in units of pressure [g cm-1 s-2].
onem = 98060._r8, & ! 1 m in units of pressure [g cm-1 s-2].
tencm = 9806._r8, & ! 10 cm in units of pressure [g cm-1 s-2].
onecm = 980.6_r8, & ! 1 cm in units of pressure [g cm-1 s-2].
onemm = 98.06_r8, & ! 1 mm in units of pressure [g cm-1 s-2].
- onemu = .09806_r8 ! 1 micrometer in units of pressure
+ onemu = .09806_r8, & ! 1 micrometer in units of pressure
! [g cm-1 s-2].
+ g2kg = 1.e-3_r8, & ! convert g to kg coeff
+ kg2g = 1.e3_r8, & ! convert kg to g coeff
+ L_mks2cgs = 1.e2_r8, & ! length coefficient converting CGS to MKS
+ M_mks2cgs = 1.e3_r8, & ! mass coefficient converting CGS to MKS
+ P_mks2cgs = 1.e1_r8, & ! pressure coefficient converting CGS to MKS
+ R_mks2cgs = 1.e-3_r8 ! rho coefficient converting CGS to MKS
+#endif
-
- public :: g, rearth, spcifh, t0deg, alpha0, pi, radian, epsil, spval, &
- tenm, onem, tencm, onecm, onemm, onemu
+ public :: g, rearth, spcifh, t0deg, alpha0, rho0, pi, radian, &
+ epsilpl, epsilp, epsilz, epsilt, epsilk, spval, &
+ tenm, onem, tencm, onecm, onemm, onemu, g2kg, kg2g, &
+ L_mks2cgs, M_mks2cgs, P_mks2cgs, R_mks2cgs
end module mod_constants
diff --git a/phy/mod_dia.F b/phy/mod_dia.F
index e824a5fd..1edf48fd 100644
--- a/phy/mod_dia.F
+++ b/phy/mod_dia.F
@@ -25,8 +25,10 @@ module mod_dia
use mod_calendar, only: date_type, date_offset, calendar_noerr
use mod_time, only: date0, date, calendar, nstep, nstep_in_day,
. nday_of_year, time, time0, baclin, dlt
- use mod_constants, only: g, spcifh, t0deg, alpha0, epsil, spval,
- . onem, onecm, onemm
+ use mod_constants, only: g, spcifh, t0deg, alpha0, epsilp, spval,
+ . onem, onecm, onemm,
+ . L_mks2cgs, M_mks2cgs, P_mks2cgs,
+ . R_mks2cgs, g2kg
use mod_xc
use mod_nctools
use netcdf, only : nf90_fill_double
@@ -169,7 +171,15 @@ module mod_dia
c
c --- Pressure thickness [g cm-1 s-2] of region for bottom salinity and
c --- temperature diagnostics
- real, parameter :: dpbot=98060.
+ real, parameter :: dpbot=onem
+c
+ real, parameter ::
+ . L_cgs2mks=1./L_mks2cgs,
+ . A_cgs2mks=1./(L_mks2cgs**2),
+ . V_cgs2mks=1./(L_mks2cgs**3),
+ . M_cgs2mks=1./M_mks2cgs,
+ . P_cgs2mks=1./P_mks2cgs,
+ . R_cgs2mks=1./R_mks2cgs
c
c --- Namelist
integer, dimension(nphymax), save ::
@@ -323,7 +333,7 @@ subroutine diafnm(ctag,diagfq,diagmon,diagann,fname)
c
date_tmp=date
c
- if (diagfq+epsil.gt.1.) then
+ if (diagfq+epsilp.gt.1.) then
errstat=date_offset(calendar,date_tmp,-1)
if (errstat.ne.calendar_noerr) then
if (mnproc.eq.1) then
@@ -968,7 +978,7 @@ subroutine diasg1
call xcbcst(j1)
do k=1,kk
call xceget(sigmar1(k),sigmar(1-nbdy,1-nbdy,k),i1,j1)
- sigmar1(k)=sigmar1(k)*1.e3 ! Convert units from g cm-3 to kg m-3
+ sigmar1(k)=sigmar1(k)*M_mks2cgs ! Convert units from g cm-3 to kg m-3
enddo
if (mnproc.eq.1) then
write(lp,*) 'sigma layers=',sigmar1
@@ -1201,7 +1211,7 @@ subroutine diaacc(m,n,mm,nn,k1m,k1n)
zup=z(i,j,kup)+.5*dz(i,j,kup)
zlo=z(i,j,k )+.5*dz(i,j,k )
tup=temp(i,j,kup+mm)
- tlo=min(temp(i,j,km),tup-epsil)
+ tlo=min(temp(i,j,km),tup-epsilp)
t20d(i,j)=(zup*(tlo-20.)+zlo*(20.-tup))/(tlo-tup)
endif
enddo
@@ -1880,7 +1890,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
enddo
c$OMP END PARALLEL DO
call xcsum(volgs(1),util1,ips)
- volgs(1)=rnacc*1e-6*volgs(1)/g
+ volgs(1)=rnacc*V_cgs2mks*volgs(1)/g
endif
if (MSC_SALNGA(iogrp).ne.0) then
c$OMP PARALLEL DO PRIVATE(l,i)
@@ -1949,7 +1959,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
tempga(1)=tempga(1)/massgs(1)
endif
if (MSC_MASSGS(iogrp).ne.0) then
- massgs(1)=rnacc*1e-3*massgs(1)/g
+ massgs(1)=rnacc*M_cgs2mks*massgs(1)/g
endif
if (MSC_SSSGA(iogrp).ne.0) then
c$OMP PARALLEL DO PRIVATE(l,i)
@@ -2004,30 +2014,30 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
c
c --- compute log10 of diffusivities
if (LYR_DIFDIA(iogrp).eq.2)
- . call loglyr(ACC_DIFDIA(iogrp),'p',1e-4,0.)
+ . call loglyr(ACC_DIFDIA(iogrp),'p',A_cgs2mks,0.)
if (LYR_DIFVMO(iogrp).eq.2)
- . call loglyr(ACC_DIFVMO(iogrp),'p',1e-4,0.)
+ . call loglyr(ACC_DIFVMO(iogrp),'p',A_cgs2mks,0.)
if (LYR_DIFVHO(iogrp).eq.2)
- . call loglyr(ACC_DIFVHO(iogrp),'p',1e-4,0.)
+ . call loglyr(ACC_DIFVHO(iogrp),'p',A_cgs2mks,0.)
if (LYR_DIFVSO(iogrp).eq.2)
- . call loglyr(ACC_DIFVSO(iogrp),'p',1e-4,0.)
+ . call loglyr(ACC_DIFVSO(iogrp),'p',A_cgs2mks,0.)
if (LYR_DIFINT(iogrp).eq.2)
- . call loglyr(ACC_DIFINT(iogrp),'p',1e-4,0.)
+ . call loglyr(ACC_DIFINT(iogrp),'p',A_cgs2mks,0.)
if (LYR_DIFISO(iogrp).eq.2)
- . call loglyr(ACC_DIFISO(iogrp),'p',1e-4,0.)
+ . call loglyr(ACC_DIFISO(iogrp),'p',A_cgs2mks,0.)
c
if (LVL_DIFDIA(iogrp).eq.2)
- . call loglvl(ACC_DIFDIALVL(iogrp),'p',1e-4*rnacc,0.)
+ . call loglvl(ACC_DIFDIALVL(iogrp),'p',A_cgs2mks*rnacc,0.)
if (LVL_DIFVMO(iogrp).eq.2)
- . call loglvl(ACC_DIFVMOLVL(iogrp),'p',1e-4*rnacc,0.)
+ . call loglvl(ACC_DIFVMOLVL(iogrp),'p',A_cgs2mks*rnacc,0.)
if (LVL_DIFVHO(iogrp).eq.2)
- . call loglvl(ACC_DIFVHOLVL(iogrp),'p',1e-4*rnacc,0.)
+ . call loglvl(ACC_DIFVHOLVL(iogrp),'p',A_cgs2mks*rnacc,0.)
if (LVL_DIFVSO(iogrp).eq.2)
- . call loglvl(ACC_DIFVSOLVL(iogrp),'p',1e-4*rnacc,0.)
+ . call loglvl(ACC_DIFVSOLVL(iogrp),'p',A_cgs2mks*rnacc,0.)
if (LVL_DIFINT(iogrp).eq.2)
- . call loglvl(ACC_DIFINTLVL(iogrp),'p',1e-4*rnacc,0.)
+ . call loglvl(ACC_DIFINTLVL(iogrp),'p',A_cgs2mks*rnacc,0.)
if (LVL_DIFISO(iogrp).eq.2)
- . call loglvl(ACC_DIFISOLVL(iogrp),'p',1e-4*rnacc,0.)
+ . call loglvl(ACC_DIFISOLVL(iogrp),'p',A_cgs2mks*rnacc,0.)
c
c --- mask sea floor of level fields
call msklvl(ACC_BFSQLVL(iogrp),'p')
@@ -2224,32 +2234,34 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
endif
c
c --- write 2d fields
- call wrth2d(ACC_SIGMX(iogrp),H2D_SIGMX(iogrp),rnacc*1e3,
+ call wrth2d(ACC_SIGMX(iogrp),H2D_SIGMX(iogrp),rnacc*R_cgs2mks,
, 0.,cmpflg,ip,'p','sigmx','Mixed layer density',' ','kg m-3')
c
- call wrth2d(ACC_UB(iogrp),H2D_UB(iogrp),rnacc*1e-2,
+ call wrth2d(ACC_UB(iogrp),H2D_UB(iogrp),rnacc*L_cgs2mks,
. 0.,cmpflg,iuu,'u','ubaro','Barotropic velocity x-component',
. ' ','m s-1')
c
- call wrth2d(ACC_VB(iogrp),H2D_VB(iogrp),rnacc*1e-2,
+ call wrth2d(ACC_VB(iogrp),H2D_VB(iogrp),rnacc*L_cgs2mks,
. 0.,cmpflg,ivv,'v','vbaro','Barotropic velocity y-component',
. ' ','m s-1')
c
call wrth2d(ACC_PSRF(iogrp),H2D_PSRF(iogrp),
- . rnacc*.1,0.,cmpflg,ip,'p','psrf','Surface pressure',' ','Pa')
+ . rnacc*P_cgs2mks,0.,cmpflg,ip,'p','psrf','Surface pressure',
+ . ' ','Pa')
c
call wrth2d(ACC_PBOT(iogrp),H2D_PBOT(iogrp),
- . rnacc*.1,0.,cmpflg,ip,'p','pbot','Bottom pressure',' ','Pa')
+ . rnacc*P_cgs2mks,0.,cmpflg,ip,'p','pbot','Bottom pressure',
+ . ' ','Pa')
c
call wrth2d(ACC_SEALV(iogrp),H2D_SEALV(iogrp),
- . -rnacc*1e-2,0.,cmpflg,ip,'p','sealv','Sea level',' ','m')
+ . -rnacc*L_cgs2mks,0.,cmpflg,ip,'p','sealv','Sea level',' ','m')
c
call wrth2d(ACC_SLVSQ(iogrp),H2D_SLVSQ(iogrp),
- . rnacc*1e-4,0.,cmpflg,ip,'p','slvsq','Sea level squared',' ',
- . 'm2')
+ . rnacc*A_cgs2mks,0.,cmpflg,ip,'p','slvsq','Sea level squared',
+ . ' ','m2')
c
call wrth2d(ACC_UTILH2D(1),H2D_BTMSTR(iogrp),
- . rnacc*0.5e-3*dlt/(g*baclin),0.,cmpflg,ip,'p','btmstr',
+ . rnacc*0.5*M_cgs2mks*dlt/(g*baclin),0.,cmpflg,ip,'p','btmstr',
. 'Barotropic mass streamfunction',' ','kg s-1')
c
call wrth2d(ACC_HICE(iogrp),H2D_HICE(iogrp),1.,0.,
@@ -2270,10 +2282,10 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
call wrth2d(ACC_IAGE(iogrp),H2D_IAGE(iogrp),1.,0.,
. cmpflg,ip,'p','iage','Ice age',' ','day')
c
- call wrth2d(ACC_UICE(iogrp),H2D_UICE(iogrp),1e-2,0.,
+ call wrth2d(ACC_UICE(iogrp),H2D_UICE(iogrp),L_cgs2mks,0.,
. cmpflg,iuu,'u','uice','Ice velocity x-component',' ','m s-1')
c
- call wrth2d(ACC_VICE(iogrp),H2D_VICE(iogrp),1e-2,0.,
+ call wrth2d(ACC_VICE(iogrp),H2D_VICE(iogrp),L_cgs2mks,0.,
. cmpflg,ivv,'v','vice','Ice velocity y-component',' ','m s-1')
c
call wrth2d(ACC_SWA(iogrp),H2D_SWA(iogrp),rnacc,0.,
@@ -2291,11 +2303,11 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'W m-2 K-1')
c
call wrth2d(ACC_SURFLX(iogrp),H2D_SURFLX(iogrp),
- . -rnacc*1e4,0.,cmpflg,ip,'p','hflx',
+ . -rnacc*L_mks2cgs*L_mks2cgs,0.,cmpflg,ip,'p','hflx',
. 'Heat flux received by ocean',' ','W m-2')
c
call wrth2d(ACC_SURRLX(iogrp),H2D_SURRLX(iogrp),
- . -rnacc*1e4,0.,cmpflg,ip,'p','hrflx',
+ . -rnacc*L_mks2cgs*L_mks2cgs,0.,cmpflg,ip,'p','hrflx',
. 'Restoring heat flux received by ocean',' ','W m-2')
c
call wrth2d(ACC_LIP(iogrp),H2D_LIP(iogrp),rnacc,0.,
@@ -2318,16 +2330,16 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. rnacc,0.,cmpflg,ip,'p','rfi','Frozen runoff',' ','kg m-2 s-1')
c
call wrth2d(ACC_SALFLX(iogrp),H2D_SALFLX(iogrp),
- . -rnacc*1e-2,0.,cmpflg,ip,'p','sflx',
+ . -rnacc*(g2kg*M_cgs2mks/A_cgs2mks),0.,cmpflg,ip,'p','sflx',
. 'Salt flux received by ocean',' ','kg m-2 s-1')
c
call wrth2d(ACC_SALRLX(iogrp),H2D_SALRLX(iogrp),
- . -rnacc*1e-2,0.,cmpflg,ip,'p','srflx',
+ . -rnacc*(g2kg*M_cgs2mks/A_cgs2mks),0.,cmpflg,ip,'p','srflx',
. 'Restoring salt flux received by ocean',' ','kg m-2 s-1')
c
call wrth2d(ACC_BRNFLX(iogrp),H2D_BRNFLX(iogrp),
- . rnacc*(-1e-2),0.,cmpflg,ip,'p','bflx','Brine flux',' ',
- . 'kg m-2 s-1')
+ . rnacc*(-g2kg*M_cgs2mks/A_cgs2mks),0.,cmpflg,ip,'p','bflx',
+ . 'Brine flux',' ','kg m-2 s-1')
c
call wrth2d(ACC_ZTX(iogrp),H2D_ZTX(iogrp),rnacc,0.,
. cmpflg,iuu,'u','ztx','Wind stress x-component',' ','N m-2')
@@ -2344,16 +2356,16 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'Momentum flux received by ocean y-component',' ','N m-2')
c
call wrth2d(ACC_IDKEDT(iogrp),H2D_IDKEDT(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','idkedt',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','idkedt',
. 'Mixed layer inertial kinetic energy tendency per unit area',
. ' ','kg s-3')
c
call wrth2d(ACC_USTAR(iogrp),H2D_USTAR(iogrp),
- . rnacc*1e-2,0.,cmpflg,ip,'p','ustar','Friction velocity',' ',
- . 'm s-1')
+ . rnacc*L_cgs2mks,0.,cmpflg,ip,'p','ustar','Friction velocity',
+ . ' ','m s-1')
c
call wrth2d(ACC_USTAR3(iogrp),H2D_USTAR3(iogrp),
- . rnacc*1.e-6,0.,cmpflg,ip,'p','ustar3',
+ . rnacc*V_cgs2mks,0.,cmpflg,ip,'p','ustar3',
. 'Friction velocity cubed',' ','m3 s-3')
c
call wrth2d(ACC_ABSWND(iogrp),H2D_ABSWND(iogrp),
@@ -2361,37 +2373,37 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'm s-1')
c
call wrth2d(ACC_MTKEUS(iogrp),H2D_MTKEUS(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','mtkeus',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','mtkeus',
. 'Mixed layer turbulent kinetic energy tendency '//
. 'per unit area related to friction velocity',
. ' ','kg s-3')
c
call wrth2d(ACC_MTKENI(iogrp),H2D_MTKENI(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','mtkeni',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','mtkeni',
. 'Mixed layer turbulent kinetic energy tendency '//
. 'per unit area related to near inertial motions',
. ' ','kg s-3')
c
call wrth2d(ACC_MTKEBF(iogrp),H2D_MTKEBF(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','mtkebf',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','mtkebf',
. 'Mixed layer turbulent kinetic energy tendency '//
. 'per unit area related to buoyancy forcing',
. ' ','kg s-3')
c
call wrth2d(ACC_MTKERS(iogrp),H2D_MTKERS(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','mtkers',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','mtkers',
. 'Mixed layer turbulent kinetic energy tendency '//
. 'per unit area related to eddy restratification',
. ' ','kg s-3')
c
call wrth2d(ACC_MTKEPE(iogrp),H2D_MTKEPE(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','mtkepe',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','mtkepe',
. 'Mixed layer turbulent kinetic energy tendency '//
. 'per unit area related to potential energy change',
. ' ','kg s-3')
c
call wrth2d(ACC_MTKEKE(iogrp),H2D_MTKEKE(iogrp),
- . rnacc*1.e-3/alpha0,0.,cmpflg,ip,'p','mtkeke',
+ . rnacc*M_cgs2mks/alpha0,0.,cmpflg,ip,'p','mtkeke',
. 'Mixed layer turbulent kinetic energy tendency '//
. 'per unit area related to kinetic energy change',
. ' ','kg s-3')
@@ -2409,23 +2421,23 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 1./onem,0.,cmpflg,ip,'p','maxmld','Maximum mixed layer depth',
. ' ','m')
c
- call wrth2d(ACC_MLTS(iogrp),H2D_MLTS(iogrp),rnacc*1e-2,
+ call wrth2d(ACC_MLTS(iogrp),H2D_MLTS(iogrp),rnacc*L_cgs2mks,
. 0.,cmpflg,ip,'p','mlts',
. 'Mixed layer thickness defined by sigma t',' ','m')
c
- call wrth2d(ACC_MLTSMN(iogrp),H2D_MLTSMN(iogrp),1e-2,
+ call wrth2d(ACC_MLTSMN(iogrp),H2D_MLTSMN(iogrp),L_cgs2mks,
. 0.,cmpflg,ip,'p','mltsmn',
. 'Minimum mixed layer thickness defined by sigma t',' ','m')
c
- call wrth2d(ACC_MLTSMX(iogrp),H2D_MLTSMX(iogrp),1e-2,
+ call wrth2d(ACC_MLTSMX(iogrp),H2D_MLTSMX(iogrp),L_cgs2mks,
. 0.,cmpflg,ip,'p','mltsmx',
. 'Maximum mixed layer thickness defined by sigma t',' ','m')
c
- call wrth2d(ACC_MLTSSQ(iogrp),H2D_MLTSSQ(iogrp),rnacc*1e-4,
+ call wrth2d(ACC_MLTSSQ(iogrp),H2D_MLTSSQ(iogrp),rnacc*A_cgs2mks,
. 0.,cmpflg,ip,'p','mltssq',
- . 'Mixed layer thickness squared defined by sigma t',' ','m')
+ . 'Mixed layer thickness squared defined by sigma t',' ','m2')
c
- call wrth2d(ACC_T20D(iogrp),H2D_T20D(iogrp),rnacc*1e-2,
+ call wrth2d(ACC_T20D(iogrp),H2D_T20D(iogrp),rnacc*L_cgs2mks,
. 0.,cmpflg,ip,'p','t20d','20C isoterm depth',' ','m')
c
call wrth2d(ACC_BRNPD(iogrp),H2D_BRNPD(iogrp),rnacc/onem,
@@ -2452,11 +2464,11 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. cmpflg,ip,'p','tbot','Bottom temperature',' ','degC')
c
c --- write 3d layer fields
- call wrtlyr(ACC_DP(iogrp),LYR_DP(iogrp),rnacc*.1,0.,
+ call wrtlyr(ACC_DP(iogrp),LYR_DP(iogrp),rnacc*P_cgs2mks,0.,
. cmpflg,ip,'p','dp','Layer pressure thickness',' ','Pa')
c
call wrtlyr(ACC_DZ(iogrp),LYR_DZ(iogrp),
- . rnacc*1e-2,0.,cmpflg,ip,'p','dz','Layer thickness',' ','m')
+ . rnacc*L_cgs2mks,0.,cmpflg,ip,'p','dz','Layer thickness',' ','m')
c
call wrtlyr(ACC_TEMP(iogrp),LYR_TEMP(iogrp),1.,0.,
. cmpflg,ip,'p','temp','Temperature','Ocean temperature',
@@ -2465,18 +2477,18 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
call wrtlyr(ACC_SALN(iogrp),LYR_SALN(iogrp),1.,0.,
. cmpflg,ip,'p','saln','Salinity','Ocean salinity','g kg-1')
c
- call wrtlyr(ACC_UVEL(iogrp),LYR_UVEL(iogrp),1e-2,
+ call wrtlyr(ACC_UVEL(iogrp),LYR_UVEL(iogrp),L_cgs2mks,
. 0.,cmpflg,iuu,'u','uvel','Velocity x-component',' ','m s-1')
c
- call wrtlyr(ACC_VVEL(iogrp),LYR_VVEL(iogrp),1e-2,
+ call wrtlyr(ACC_VVEL(iogrp),LYR_VVEL(iogrp),L_cgs2mks,
. 0.,cmpflg,ivv,'v','vvel','Velocity y-component',' ','m s-1')
c
call wrtlyr(ACC_UFLX(iogrp),LYR_UFLX(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,iuu,'u','uflx',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','uflx',
. 'Mass flux in x-direction',' ','kg s-1')
c
call wrtlyr(ACC_VFLX(iogrp),LYR_VFLX(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,ivv,'v','vflx',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vflx',
. 'Mass flux in y-direction',' ','kg s-1')
c
call wrtlyr(ACC_UTFLX(iogrp),LYR_UTFLX(iogrp),
@@ -2488,20 +2500,20 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'Heat flux in y-direction',' ','W')
c
call wrtlyr(ACC_USFLX(iogrp),LYR_USFLX(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,iuu,'u','usflx',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','usflx',
. 'Salt flux in x-direction',' ','kg s-1')
c
call wrtlyr(ACC_VSFLX(iogrp),LYR_VSFLX(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,ivv,'v','vsflx',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vsflx',
. 'Salt flux in y-direction',' ','kg s-1')
c
call wrtlyr(ACC_UMFLTD(iogrp),LYR_UMFLTD(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,iuu,'u','umfltd',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','umfltd',
. 'Mass flux due to thickness diffusion in x-direction',' ',
. 'kg s-1')
c
call wrtlyr(ACC_VMFLTD(iogrp),LYR_VMFLTD(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,ivv,'v','vmfltd',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vmfltd',
. 'Mass flux due to thickness diffusion in y-direction',' ',
. 'kg s-1')
c
@@ -2526,31 +2538,31 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'W')
c
call wrtlyr(ACC_USFLTD(iogrp),LYR_USFLTD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,iuu,'u','usfltd',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','usfltd',
. 'Salt flux due to thickness diffusion in x-direction',' ',
. 'kg s-1')
c
call wrtlyr(ACC_VSFLTD(iogrp),LYR_VSFLTD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,ivv,'v','vsfltd',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vsfltd',
. 'Salt flux due to thickness diffusion in y-direction',' ',
. 'kg s-1')
c
call wrtlyr(ACC_USFLLD(iogrp),LYR_USFLLD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,iuu,'u','usflld',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','usflld',
. 'Salt flux due to lateral diffusion in x-direction',' ',
. 'kg s-1')
c
call wrtlyr(ACC_VSFLLD(iogrp),LYR_VSFLLD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,ivv,'v','vsflld',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vsflld',
. 'Salt flux due to lateral diffusion in y-direction',' ',
. 'kg s-1')
c
call wrtlyr(ACC_WFLX(iogrp),LYR_WFLX(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,ip,'p','wflx',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,ip,'p','wflx',
. 'Vertical mass flux',' ','kg s-1')
c
call wrtlyr(ACC_WFLX2(iogrp),LYR_WFLX2(iogrp),
- . rnacc*(0.5e-3/(g*baclin))**2,0.,cmpflg,ip,'p','wflx2',
+ . rnacc*(0.5*M_cgs2mks/(g*baclin))**2,0.,cmpflg,ip,'p','wflx2',
. 'Vertical mass flux squared',' ','kg2 s-2')
c
call wrtlyr(ACC_BFSQ(iogrp),LYR_BFSQ(iogrp),1.,0.,
@@ -2558,7 +2570,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 's-1')
c
call wrtlyr(ACC_AVDSG(iogrp),LYR_PV(iogrp),
- . 1.e2*g,0.,cmpflg,ip,'p','pv','Potential vorticity',' ',
+ . L_mks2cgs*g,0.,cmpflg,ip,'p','pv','Potential vorticity',' ',
. 'm-1 s-1')
c
if (LYR_DIFINT(iogrp).eq.2) then
@@ -2566,7 +2578,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difint','Layer interface diffusivity',
. ' ','log10(m2 s-1)')
else
- call wrtlyr(ACC_DIFINT(iogrp),LYR_DIFINT(iogrp),1e-4,
+ call wrtlyr(ACC_DIFINT(iogrp),LYR_DIFINT(iogrp),A_cgs2mks,
. 0.,cmpflg,ip,'p','difint','Layer interface diffusivity',
. ' ','m2 s-1')
endif
@@ -2576,7 +2588,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difiso','Isopycnal diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlyr(ACC_DIFISO(iogrp),LYR_DIFISO(iogrp),1e-4,
+ call wrtlyr(ACC_DIFISO(iogrp),LYR_DIFISO(iogrp),A_cgs2mks,
. 0.,cmpflg,ip,'p','difiso','Isopycnal diffusivity',' ',
. 'm2 s-1')
endif
@@ -2586,7 +2598,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difdia','Vertical diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlyr(ACC_DIFDIA(iogrp),LYR_DIFDIA(iogrp),1e-4,
+ call wrtlyr(ACC_DIFDIA(iogrp),LYR_DIFDIA(iogrp),A_cgs2mks,
. 0.,cmpflg,ip,'p','difdia','Vertical diffusivity',' ',
. 'm2 s-1')
endif
@@ -2596,7 +2608,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difvmo','Vertical momentum diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlyr(ACC_DIFVMO(iogrp),LYR_DIFVMO(iogrp),1e-4,
+ call wrtlyr(ACC_DIFVMO(iogrp),LYR_DIFVMO(iogrp),A_cgs2mks,
. 0.,cmpflg,ip,'p','difvmo','Vertical momentum diffusivity',' ',
. 'm2 s-1')
endif
@@ -2606,7 +2618,7 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difvho','Vertical heat diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlyr(ACC_DIFVHO(iogrp),LYR_DIFVHO(iogrp),1e-4,
+ call wrtlyr(ACC_DIFVHO(iogrp),LYR_DIFVHO(iogrp),A_cgs2mks,
. 0.,cmpflg,ip,'p','difvho','Vertical heat diffusivity',' ',
. 'm2 s-1')
endif
@@ -2616,24 +2628,25 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difvso','Vertical salt diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlyr(ACC_DIFVSO(iogrp),LYR_DIFVSO(iogrp),1e-4,
+ call wrtlyr(ACC_DIFVSO(iogrp),LYR_DIFVSO(iogrp),A_cgs2mks,
. 0.,cmpflg,ip,'p','difvso','Vertical salt diffusivity',' ',
. 'm2 s-1')
endif
c
#if defined(TRC) && defined(TKE)
- call wrtlyr(ACC_TKE(iogrp),LYR_TKE(iogrp),1e-4,0.,
+ call wrtlyr(ACC_TKE(iogrp),LYR_TKE(iogrp),A_cgs2mks,0.,
. cmpflg,ip,'p','tke','TKE','Turbulent kinetic energy',
. 'm2 s-2')
c
- call wrtlyr(ACC_GLS_PSI(iogrp),LYR_GLS_PSI(iogrp),1.e-4,0.,
+ call wrtlyr(ACC_GLS_PSI(iogrp),LYR_GLS_PSI(iogrp),A_cgs2mks,0.,
. cmpflg,ip,'p','gls_psi','GLS_PSI','Generic length scale',
. 'm2 s-3')
c
#endif
c --- Write 3d depth fields
call wrtlvl(ACC_DZLVL(iogrp),LVL_DZ(iogrp),
- . rnacc*1e-2,0.,cmpflg,ip,'p','dzlvl','Layer thickness',' ','m')
+ . rnacc*L_cgs2mks,0.,cmpflg,ip,
+ . 'p','dzlvl','Layer thickness',' ','m')
c
call wrtlvl(ACC_TEMPLVL(iogrp),LVL_TEMP(iogrp),
. rnacc,0.,cmpflg,ip,'p','templvl','Temperature',
@@ -2644,19 +2657,19 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'Ocean salinity','g kg-1')
c
call wrtlvl(ACC_UVELLVL(iogrp),LVL_UVEL(iogrp),
- . rnacc*1e-2,0.,cmpflg,iuu,'u','uvellvl',
+ . rnacc*L_cgs2mks,0.,cmpflg,iuu,'u','uvellvl',
. 'Velocity x-component',' ','m s-1')
c
call wrtlvl(ACC_VVELLVL(iogrp),LVL_VVEL(iogrp),
- . rnacc*1e-2,0.,cmpflg,ivv,'v','vvellvl',
+ . rnacc*L_cgs2mks,0.,cmpflg,ivv,'v','vvellvl',
. 'Velocity y-component',' ','m s-1')
c
call wrtlvl(ACC_UFLXLVL(iogrp),LVL_UFLX(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,iuu,'u','uflxlvl',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','uflxlvl',
. 'Mass flux in x-direction',' ','kg s-1')
c
call wrtlvl(ACC_VFLXLVL(iogrp),LVL_VFLX(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,ivv,'v','vflxlvl',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vflxlvl',
. 'Mass flux in y-direction',' ','kg s-1')
c
call wrtlvl(ACC_UTFLXLVL(iogrp),LVL_UTFLX(iogrp),
@@ -2668,20 +2681,20 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'Heat flux in y-direction',' ','W')
c
call wrtlvl(ACC_USFLXLVL(iogrp),LVL_USFLX(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,iuu,'u','usflxlvl',
- . 'Salt flux in x-direction',' ','kg s-1')
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u',
+ . 'usflxlvl','Salt flux in x-direction',' ','kg s-1')
c
call wrtlvl(ACC_VSFLXLVL(iogrp),LVL_VSFLX(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,ivv,'v','vsflxlvl',
- . 'Salt flux in y-direction',' ','kg s-1')
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v',
+ . 'vsflxlvl','Salt flux in y-direction',' ','kg s-1')
c
call wrtlvl(ACC_UMFLTDLVL(iogrp),LVL_UMFLTD(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,iuu,'u','umfltdlvl',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,'u','umfltdlvl',
. 'Mass flux due to thickness diffusion in x-direction',' ',
. 'kg s-1')
c
call wrtlvl(ACC_VMFLTDLVL(iogrp),LVL_VMFLTD(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,ivv,'v','vmfltdlvl',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,'v','vmfltdlvl',
. 'Mass flux due to thickness diffusion in y-direction',' ',
. 'kg s-1')
c
@@ -2706,31 +2719,35 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 'W')
c
call wrtlvl(ACC_USFLTDLVL(iogrp),LVL_USFLTD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,iuu,'u','usfltdlvl',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,
+ . 'u','usfltdlvl',
. 'Salt flux due to thickness diffusion in x-direction',' ',
. 'kg s-1')
c
call wrtlvl(ACC_VSFLTDLVL(iogrp),LVL_VSFLTD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,ivv,'v','vsfltdlvl',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,
+ . 'v','vsfltdlvl',
. 'Salt flux due to thickness diffusion in y-direction',' ',
. 'kg s-1')
c
call wrtlvl(ACC_USFLLDLVL(iogrp),LVL_USFLLD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,iuu,'u','usflldlvl',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,iuu,
+ . 'u','usflldlvl',
. 'Salt flux due to lateral diffusion in x-direction',' ',
. 'kg s-1')
c
call wrtlvl(ACC_VSFLLDLVL(iogrp),LVL_VSFLLD(iogrp),
- . rnacc*0.5e-6/(g*baclin),0.,cmpflg,ivv,'v','vsflldlvl',
+ . rnacc*0.5*g2kg*M_cgs2mks/(g*baclin),0.,cmpflg,ivv,
+ . 'v','vsflldlvl',
. 'Salt flux due to lateral diffusion in y-direction',' ',
. 'kg s-1')
c
call wrtlvl(ACC_WFLXLVL(iogrp),LVL_WFLX(iogrp),
- . rnacc*0.5e-3/(g*baclin),0.,cmpflg,ip,'p','wflxlvl',
+ . rnacc*0.5*M_cgs2mks/(g*baclin),0.,cmpflg,ip,'p','wflxlvl',
. 'Vertical mass flux',' ','kg s-1')
c
call wrtlvl(ACC_WFLX2LVL(iogrp),LVL_WFLX2(iogrp),
- . rnacc*(0.5e-3/(g*baclin))**2,0.,cmpflg,ip,'p','wflx2lvl',
+ . rnacc*(0.5*M_cgs2mks/(g*baclin))**2,0.,cmpflg,ip,'p','wflx2lvl',
. 'Vertical mass flux squared',' ','kg2 s-2')
c
call wrtlvl(ACC_BFSQLVL(iogrp),LVL_BFSQ(iogrp),
@@ -2738,17 +2755,17 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. ' ','s-1')
c
call wrtlvl(ACC_PVLVL(iogrp),LVL_PV(iogrp),
- . rnacc*1.e2*g,0.,cmpflg,ip,'p','pvlvl','Potential vorticity',
- . ' ','m-1 s-1')
+ . rnacc*L_mks2cgs*g,0.,cmpflg,ip,
+ . 'p','pvlvl','Potential vorticity',' ','m-1 s-1')
c
if (LVL_DIFINT(iogrp).eq.2) then
call wrtlvl(ACC_DIFINTLVL(iogrp),LVL_DIFINT(iogrp),1.,
. 0.,cmpflg,ip,'p','difintlvl','Layer interface diffusivity',
. ' ','log10(m2 s-1)')
else
- call wrtlvl(ACC_DIFINTLVL(iogrp),LVL_DIFINT(iogrp),1e-4*rnacc,
- . 0.,cmpflg,ip,'p','difintlvl','Layer interface diffusivity',
- . ' ','m2 s-1')
+ call wrtlvl(ACC_DIFINTLVL(iogrp),LVL_DIFINT(iogrp),
+ . A_cgs2mks*rnacc,0.,cmpflg,ip,'p','difintlvl',
+ . 'Layer interface diffusivity',' ','m2 s-1')
endif
c
if (LVL_DIFISO(iogrp).eq.2) then
@@ -2756,9 +2773,9 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difisolvl','Isopycnal diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlvl(ACC_DIFISOLVL(iogrp),LVL_DIFISO(iogrp),1e-4*rnacc,
- . 0.,cmpflg,ip,'p','difisolvl','Isopycnal diffusivity',' ',
- . 'm2 s-1')
+ call wrtlvl(ACC_DIFISOLVL(iogrp),LVL_DIFISO(iogrp),
+ . A_cgs2mks*rnacc,0.,cmpflg,ip,'p','difisolvl',
+ . 'Isopycnal diffusivity',' ','m2 s-1')
endif
c
if (LVL_DIFDIA(iogrp).eq.2) then
@@ -2766,9 +2783,9 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difdialvl','Vertical diffusivity',' ',
. 'log10(m2 s-1)')
else
- call wrtlvl(ACC_DIFDIALVL(iogrp),LVL_DIFDIA(iogrp),1e-4*rnacc,
- . 0.,cmpflg,ip,'p','difdialvl','Vertical diffusivity',' ',
- . 'm2 s-1')
+ call wrtlvl(ACC_DIFDIALVL(iogrp),LVL_DIFDIA(iogrp),
+ . A_cgs2mks*rnacc,0.,cmpflg,ip,'p','difdialvl',
+ . 'Vertical diffusivity',' ','m2 s-1')
endif
c
if (LVL_DIFVMO(iogrp).eq.2) then
@@ -2776,9 +2793,9 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difvmolvl','Vertical momentum diffusivity',
. ' ','log10(m2 s-1)')
else
- call wrtlvl(ACC_DIFVMOLVL(iogrp),LVL_DIFVMO(iogrp),1e-4*rnacc,
- . 0.,cmpflg,ip,'p','difvmolvl','Vertical momentum diffusivity',
- . ' ','m2 s-1')
+ call wrtlvl(ACC_DIFVMOLVL(iogrp),LVL_DIFVMO(iogrp),
+ . A_cgs2mks*rnacc,0.,cmpflg,ip,'p','difvmolvl',
+ . 'Vertical momentum diffusivity',' ','m2 s-1')
endif
c
if (LVL_DIFVHO(iogrp).eq.2) then
@@ -2786,9 +2803,9 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difvholvl','Vertical heat diffusivity',
. ' ','log10(m2 s-1)')
else
- call wrtlvl(ACC_DIFVHOLVL(iogrp),LVL_DIFVHO(iogrp),1e-4*rnacc,
- . 0.,cmpflg,ip,'p','difvholvl','Vertical heat diffusivity',
- . ' ','m2 s-1')
+ call wrtlvl(ACC_DIFVHOLVL(iogrp),LVL_DIFVHO(iogrp),
+ . A_cgs2mks*rnacc,0.,cmpflg,ip,'p','difvholvl',
+ . 'Vertical heat diffusivity',' ','m2 s-1')
endif
c
if (LVL_DIFVSO(iogrp).eq.2) then
@@ -2796,19 +2813,19 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
. 0.,cmpflg,ip,'p','difvsolvl','Vertical salt diffusivity',
. ' ','log10(m2 s-1)')
else
- call wrtlvl(ACC_DIFVSOLVL(iogrp),LVL_DIFVSO(iogrp),1e-4*rnacc,
- . 0.,cmpflg,ip,'p','difvsolvl','Vertical salt diffusivity',
- . ' ','m2 s-1')
+ call wrtlvl(ACC_DIFVSOLVL(iogrp),LVL_DIFVSO(iogrp),
+ . A_cgs2mks*rnacc,0.,cmpflg,ip,'p','difvsolvl',
+ . 'Vertical salt diffusivity',' ','m2 s-1')
endif
c
#if defined(TRC) && defined(TKE)
- call wrtlvl(ACC_TKELVL(iogrp),LVL_TKE(iogrp),rnacc*1.e-4,
+ call wrtlvl(ACC_TKELVL(iogrp),LVL_TKE(iogrp),rnacc*A_cgs2mks,
. 0.,cmpflg,ip,'p','tkelvl','Turbulent kinetic energy',' ',
. 'm2 s-2')
c
- call wrtlvl(ACC_GLS_PSILVL(iogrp),LVL_GLS_PSI(iogrp),rnacc*1.e-4,
- . 0.,cmpflg,ip,'p','gls_psilvl','Generic length scale',' ',
- . 'm2 s-3')
+ call wrtlvl(ACC_GLS_PSILVL(iogrp),LVL_GLS_PSI(iogrp),
+ . rnacc*A_cgs2mks,0.,cmpflg,ip,'p','gls_psilvl',
+ . 'Generic length scale',' ','m2 s-3')
c
#endif
c
@@ -2919,8 +2936,8 @@ subroutine diaout(iogrp,m,n,mm,nn,k1m,k1n)
call inilyr(ACC_UTILLYR(1),'p',0.)
call acclyr(ACC_UTILLYR,dp(1-nbdy,1-nbdy,k1m),tmp3d,0,'p')
call wrtlyr(ACC_UTILLYR(1),
- . max(LYR_IDLAGE(iogrp),LYR_TRC(iogrp)),.1,0.,cmpflg,ip,'p',
- . 'dp_trc','Layer pressure thickness',' ','Pa')
+ . max(LYR_IDLAGE(iogrp),LYR_TRC(iogrp)),P_cgs2mks,0.,cmpflg,ip,
+ . 'p','dp_trc','Layer pressure thickness',' ','Pa')
endif
# ifdef IDLAGE
c
@@ -3134,14 +3151,14 @@ subroutine diasec(iogrp)
do i=max(1,ifu(j,l)),min(ii,ilu(j,l))
uflx_cum(i,j)=uflx_cum(i,j)+
. phylvl(i,j,k,ACC_UFLXLVL(iogrp))
- . *0.5e-3/(g*baclin*nacc_phy(iogrp))
+ . *0.5*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
enddo
enddo
do l=1,isv(j)
do i=max(1,ifv(j,l)),min(ii,ilv(j,l))
vflx_cum(i,j)=vflx_cum(i,j)+
. phylvl(i,j,k,ACC_VFLXLVL(iogrp))
- . *0.5e-3/(g*baclin*nacc_phy(iogrp))
+ . *0.5*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
enddo
enddo
c
@@ -3370,17 +3387,17 @@ subroutine diamer(iogrp)
if (ACC_MSFLX(iogrp).eq.0) exit
ACC_UIND=ACC_USFLX(iogrp)
ACC_VIND=ACC_VSFLX(iogrp)
- r=0.5e-6/(g*baclin*nacc_phy(iogrp))
+ r=0.5*g2kg*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
elseif (nfld.eq.5) then
if (ACC_MSFTD(iogrp).eq.0) exit
ACC_UIND=ACC_USFLTD(iogrp)
ACC_VIND=ACC_VSFLTD(iogrp)
- r=0.5e-6/(g*baclin*nacc_phy(iogrp))
+ r=0.5*g2kg*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
elseif (nfld.eq.6) then
if (ACC_MSFLD(iogrp).eq.0) exit
ACC_UIND=ACC_USFLLD(iogrp)
ACC_VIND=ACC_VSFLLD(iogrp)
- r=0.5e-6/(g*baclin*nacc_phy(iogrp))
+ r=0.5*g2kg*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
else
write(lp,*) 'field index out of range'
call xchalt('(diamer)')
@@ -3541,7 +3558,7 @@ subroutine diamer(iogrp)
enddo
c$OMP END PARALLEL DO
c
- r=0.5e-3/(g*baclin*nacc_phy(iogrp))
+ r=0.5*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
c
do nfld=1,2
c
@@ -3643,7 +3660,7 @@ subroutine diamer(iogrp)
enddo
endif
if (abs(mflx_mr(l,m)-mflx_last_mr(l,m)).lt.
- . 1.e5*epsil) then
+ . 1.e5*epsilp) then
mflx_last_mr(l,m)=mflx_mr(l,m)
mflx_mr(l,m)=nf90_fill_double
else
@@ -3706,7 +3723,7 @@ subroutine diamer(iogrp)
enddo
endif
c
- r=0.5e-3/(g*baclin*nacc_phy(iogrp))
+ r=0.5*M_cgs2mks/(g*baclin*nacc_phy(iogrp))
c
do nfld=1,2
c
@@ -3947,9 +3964,9 @@ subroutine diazlv(gridid,k,mm,nn,ind1,ind2,weights,weightsflx)
save ztop,zbot,dlevp,dlevu,dlevv,iniflg
c
c --- Define thresholds
- dzeps=1e1*epsil
- dpeps=1e5*epsil
- flxeps=1e5*epsil
+ dzeps=1e1*epsilp
+ dpeps=1e5*epsilp
+ flxeps=1e5*epsilp
c
c --- Sort out stuff related to time stepping
km=k+mm
diff --git a/phy/mod_difest.F b/phy/mod_difest.F
index 3831bf38..67f0f486 100644
--- a/phy/mod_difest.F
+++ b/phy/mod_difest.F
@@ -20,7 +20,8 @@
module mod_difest
c
use mod_types, only: r8
- use mod_constants, only: g, alpha0, pi, epsil, spval, onem, onecm
+ use mod_constants, only: g, alpha0, pi, epsilp, spval, onem,
+ . onecm, L_mks2cgs, M_mks2cgs, R_mks2cgs
use mod_time, only: delt1
use mod_xc
use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml,
@@ -101,6 +102,12 @@ module mod_difest
type(CVMix_kpp_params_type) :: KPP_params
c type(CVMix_kpp_params_type), pointer :: CVmix_kpp_params_in
c type(CVMix_kpp_params_type) :: CVmix_kpp_params_in
+c
+ real, parameter ::
+ . iL_mks2cgs = 1./L_mks2cgs,
+ . iM_mks2cgs = 1./M_mks2cgs,
+ . A_mks2cgs = L_mks2cgs**2,
+ . A_cgs2mks = 1./(L_mks2cgs*L_mks2cgs)
c
c --- parameters:
c --- iidtyp - type of interface and isopycnal diffusivities. If
@@ -162,22 +169,33 @@ module mod_difest
c --- non-isopycnic layers [g/cm/s**2].
c --- dpnbav - thickness of region near the bottom used to estimate
c --- bottom Brunt-Vaisala frequency [g/cm/s**2].
+c --- cpsemin - Lower bound of zonal eddy phase speed minus zonal
+c --- barotropic velocity [cm/s].
+c --- urmsemin- Lower bound of absolute value of RMS eddy velocity
+c --- [cm/s].
integer iidtyp,bdmldp,tdmflg,iwdflg
real dptmin,dpbmin,drhomn,thkdff,temdff,nu0,nus0,nug0,drho0,nuls0,
. iwdfac,dmxeff,tdmq,tdmls0,tdmls1,tdclat,tddlat,tkepls,niwls,
. cori30,bvf0,nubmin,dpgc,dpgrav,dpdiav,dpddav,dpnbav,ustmin,
- . kappa,bfeps,sleps,zetas,as,cs,minOBLdepth
- parameter (iidtyp=2,bdmldp=1,tdmflg=1,iwdflg=1,dptmin=98060.,
- . dpbmin=980.6,drhomn=6.e-6,thkdff=.5,temdff=.35,nu0=.1,
- . nus0=50.,nug0=2500.,drho0=6.e-6,nuls0=500.,iwdfac=.06,
- . dmxeff=.2,tdmq=1./3.,tdmls0=500.*98060.,
- . tdmls1=100.*98060.,tdclat=74.5,tddlat=3.,
- . tkepls=20.*98060.,niwls=300.*98060.,cori30=7.2722e-5,
- . bvf0=5.24e-3,nubmin=.01,dpgc=300.*98060.,
- . dpgrav=100.*98060.,dpdiav=100.*98060.,
- . dpddav=10.*98060.,dpnbav=250.*98060.,ustmin=.1,
- . kappa=.4,bfeps=1.e-12,sleps=.1,zetas=-1.,as=-28.86,
- . cs=98.96,minOBLdepth=1.0)
+ . kappa,bfeps,sleps,zetas,as,cs,minOBLdepth,
+ . cpsemin,urmsemin
+ parameter (iidtyp=2,bdmldp=1,tdmflg=1,iwdflg=1,dptmin=onem,
+ . dpbmin=onecm,drhomn=6.e-3*R_mks2cgs,
+ . thkdff=5.e-3*L_mks2cgs,temdff=3.5e-3*L_mks2cgs,
+ . nu0=1.e-5*A_mks2cgs,
+ . nus0=5.e-3*A_mks2cgs,
+ . nug0=2.5e-1*A_mks2cgs,
+ . drho0=6.e-3*R_mks2cgs,
+ . nuls0=5.e-2*A_mks2cgs,iwdfac=.06,
+ . dmxeff=.2,tdmq=1./3.,tdmls0=500.*onem,
+ . tdmls1=100.*onem,tdclat=74.5,tddlat=3.,
+ . tkepls=20.*onem,niwls=300.*onem,cori30=7.2722e-5,
+ . bvf0=5.24e-3,nubmin=1.e-6*A_mks2cgs,
+ . dpgc=300.*onem,dpgrav=100.*onem,dpdiav=100.*onem,
+ . dpddav=10.*onem,dpnbav=250.*onem,ustmin=.001*L_mks2cgs,
+ . kappa=.4,bfeps=1.e-16*A_mks2cgs,sleps=.1,zetas=-1.,
+ . cpsemin=-0.2*L_mks2cgs,urmsemin=0.05*L_mks2cgs,
+ . as=-28.86,cs=98.96,minOBLdepth=1.0)
c
public :: OBLdepth, inivar_difest, init_difest, difest_isobml,
. difest_lateral_hybrid, difest_vertical_hybrid, hOBL
@@ -236,8 +254,8 @@ subroutine init_difest
c --- ------ convection routine based on N2 not rho
c --- ------ if lBruntVaisala is TRUE, otherwise based on rho
c --- ------ convert nuls0 to m2/s
- call CVMix_init_conv(convect_diff=20.0*nuls0*1e-4,
- . convect_visc=20.0*nuls0*1e-4,
+ call CVMix_init_conv(convect_diff=20.0*nuls0*A_cgs2mks,
+ . convect_visc=20.0*nuls0*A_cgs2mks,
. lBruntVaisala=.true.,
. BVsqr_convect=0.0)
call CVMix_put(CVMix_glb_params,'max_nlev',kk)
@@ -245,7 +263,7 @@ subroutine init_difest
call CVMix_put(CVMix_glb_params,'FreshWaterDensity',1000.0)
call CVMix_put(CVMix_glb_params,'SaltWaterDensity',1025.0)
call cvmix_init_shear(mix_scheme='KPP',
- . KPP_nu_zero=nus0*1e-4,
+ . KPP_nu_zero=nus0*A_cgs2mks,
. KPP_Ri_zero=ri0,
. KPP_exp=3.0)
! CVmix_kpp_params_in => CVmix_kpp_params_user
@@ -501,7 +519,7 @@ subroutine difest_common_iso(m,n,mm,nn,k1m,k1n)
elseif (k.lt.kmax(i,j)) then
q=max(0.,rho(p(i,j,k+1),temp(i,j,kn+1),saln(i,j,kn+1))
. -rho(p(i,j,k+1),temp(i,j,kn ),saln(i,j,kn )))
- drhol(i,j,k)=2.*tup(i)*q/max(1.e-14,tup(i)+q)
+ drhol(i,j,k)=2.*tup(i)*q/max(1.e-11*R_mks2cgs,tup(i)+q)
tup(i)=q
else
drhol(i,j,k)=tup(i)
@@ -518,7 +536,7 @@ subroutine difest_common_iso(m,n,mm,nn,k1m,k1n)
c --- ------- Local gradient Richardson number.
rig(i,j,k)=alpha0*alpha0
. *max(drhomn,drhol(i,j,k))*dp(i,j,kn)
- . /max(1.e-9,du2l(i,j,k))
+ . /max(1.e-13*A_mks2cgs,du2l(i,j,k))
c
endif
enddo
@@ -640,7 +658,8 @@ subroutine difest_common_hyb(m,n,mm,nn,k1m,k1n)
. +(mskv(i,j,k)*dv2(i,j,k)+mskv(i,j+1,k)*dv2(i,j+1,k))
. /max(1,mskv(i,j,k)+mskv(i,j+1,k))
dz=.5*(dp(i,j,kn-1)+dp(i,j,kn))*alpha0/g
- rig(i,j,k)=max(0.,bfsqi(i,j,k)*dz*dz)/max(1.e-9,q)
+ rig(i,j,k)=max(0.,bfsqi(i,j,k)*dz*dz)
+ . /max(1.e-13*A_mks2cgs,q)
else
rig(i,j,k)=rig(i,j,k-1)
endif
@@ -964,16 +983,16 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
Kt_kpp = 0.0
Ks_kpp = 0.0
do k=1,kk+1
- Kv_kpp(k) = Kvisc_m(i,j,k)*1e-4
- Kt_kpp(k) = Kdiff_t(i,j,k)*1e-4
- Ks_kpp(k) = Kdiff_s(i,j,k)*1e-4
+ Kv_kpp(k) = Kvisc_m(i,j,k)*A_cgs2mks
+ Kt_kpp(k) = Kdiff_t(i,j,k)*A_cgs2mks
+ Ks_kpp(k) = Kdiff_s(i,j,k)*A_cgs2mks
enddo
depth_int(1) = p(i,j,1)/onem
iFaceHeight(1) = -depth_int(1)
! convert cm/s to m/s
- surfFricVel = ustar(i,j) * 1e-2
+ surfFricVel = ustar(i,j) * iL_mks2cgs
! convert cm2/s3 to m2/s3
- surfBuoyFlux = - buoyfl(i,j,1) * 1e-4
+ surfBuoyFlux = - buoyfl(i,j,1) * A_cgs2mks
do k=1,kk
kn = k + nn
kn1 = max(nn+1,kn-1)
@@ -1028,7 +1047,7 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
surfU = surfHu / hTot
surfV = surfHv / hTot
surfRho = rho(p(i,j,k),surfTemp,surfSalt)
- if (p(i,j,kk+1)-p(i,j,k) < epsil) then
+ if (p(i,j,kk+1)-p(i,j,k) < epsilp) then
deltaRho(k) = deltaRho(k-1)
else
deltaRho(k) = rho_1d(k) - surfRho
@@ -1056,12 +1075,12 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
rig_i(k)=rig(i,j,k)
surfBuoyFlux2(k) = ( buoyfl(i,j,k+1)
- . - buoyfl(i,j,1 )) * 1e-4
+ . - buoyfl(i,j,1 )) * A_cgs2mks
c
enddo ! k
if(dps.gt.0.) bvfbot=bvfbot/dps
! convert cm2/s2 to m2/s2
- deltaU2 = deltaU2*1e-4
+ deltaU2 = deltaU2*A_cgs2mks
! bottom values for the Ri, N2, and N
rig_i(kk+1) = rig_i(kk)
@@ -1079,8 +1098,8 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
elseif (bdmtyp.eq.2) then
c --- --------- Type 2: Background diffusivity is a constant
! convert cm2/s2 to m2/s2
- Kv_col(:) = bdmc2*1e-4
- Kd_col(:) = bdmc2*1e-4
+ Kv_col(:) = bdmc2*A_cgs2mks
+ Kd_col(:) = bdmc2*A_cgs2mks
else
Kv_col(:) = 0.
Kd_col(:) = 0.
@@ -1099,7 +1118,7 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
. efficiency=dmxeff, local_mixing_frac=tdmq)
call CVMix_compute_Simmons_invariant(nlev=kk,
- . energy_flux=twedon(i,j)*bvfbot*1e-3,
+ . energy_flux=twedon(i,j)*bvfbot*iM_mks2cgs,
. rho=CVMix_glb_params%FreshWaterDensity,
. SimmonsCoeff = Simmons_coeff, VertDep = vert_dep,
. zw = iFaceHeight, zt = cellHeight,
@@ -1146,7 +1165,7 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
! Calculate Bulk Richardson number from eq (21) of LMD94
BulkRi_1d = CVmix_kpp_compute_bulk_Richardson(
. zt_cntr = cellHeight, ! Depth of cell center [m]
- . delta_buoy_cntr=g*alpha0*deltaRho*1e-2, ! Bulk buoyancy difference, Br-B(z) [m s-2]
+ . delta_buoy_cntr=g*alpha0*deltaRho*iL_mks2cgs, ! Bulk buoyancy difference, Br-B(z) [m s-2]
. delta_Vsqr_cntr=deltaU2, ! Square of resolved velocity difference [m2 s-2]
. Vt_sqr_cntr=VT2(:), ! Unresolved shear [m2 s-2]
. ws_cntr=Ws_1d, ! Turbulent velocity scale profile [m s-1]
@@ -1201,7 +1220,7 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
! Buoyancy flux acting on the OBL
surfBuoyFlux = ( buoyfl(i,j,kOBL+1)
- . - buoyfl(i,j,1 )) * 1e-4
+ . - buoyfl(i,j,1 )) * A_cgs2mks
! Compute KPP using CVMix
call CVMix_coeffs_kpp(Kv_kpp(:), ! (inout) Total viscosity [m2 s-1]
@@ -1224,9 +1243,9 @@ subroutine difest_vertical_hyb(m,n,mm,nn,k1m,k1n)
c ---- ccc -------
! convert m2/s to cm2/s
- Kv_kpp = Kv_kpp*1e4
- Kt_kpp = Kt_kpp*1e4
- Ks_kpp = Ks_kpp*1e4
+ Kv_kpp = Kv_kpp*A_mks2cgs
+ Kt_kpp = Kt_kpp*A_mks2cgs
+ Ks_kpp = Ks_kpp*A_mks2cgs
Kv_kpp=max(nubmin,Kv_kpp)
Kt_kpp=max(nubmin,Kt_kpp)
Ks_kpp=max(nubmin,Ks_kpp)
@@ -1295,7 +1314,7 @@ subroutine difest_lateral_hyb(m,n,mm,nn,k1m,k1n)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
kfil(i,j)=kk+1
do k=kk,2,-1
- if (p(i,j,k).gt.mlts(i,j)*onecm) kfil(i,j)=k
+ if (p(i,j,k).gt.mlts(i,j)*(onem*iL_mks2cgs)) kfil(i,j)=k
enddo
enddo
enddo
@@ -1322,7 +1341,7 @@ subroutine difest_lateral_hyb(m,n,mm,nn,k1m,k1n)
kn=k+nn
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- if (p(i,j,kk+1)-p(i,j,k+1).lt.epsil) then
+ if (p(i,j,kk+1)-p(i,j,k+1).lt.epsilp) then
plo=p(i,j,kk+1)
else
plo=.5*(p(i,j,k)+p(i,j,k+1))
@@ -1519,7 +1538,7 @@ subroutine difest_lateral_hyb(m,n,mm,nn,k1m,k1n)
. kmax(i,j)-kfil(i,j).ge.1) then
c
c --- --------- Rhines scale.
- rhisc=egr(i,k)/max(1.e-24,betafp(i,j))
+ rhisc=egr(i,k)/max(1.e-22*iL_mks2cgs,betafp(i,j))
c
c --- --------- Eddy length scale.
els=max(eglsmn,min(bcrrd(i),rhisc))
@@ -1582,13 +1601,13 @@ subroutine difest_lateral_hyb(m,n,mm,nn,k1m,k1n)
c --- ------- RMS eddy velocity estimated from K = Gamma*u_rms*L, where
c --- ------- a mixing efficiency of Gamma = 0.35 is used (Klocker and
c --- ------- Abernathey, 2014).
- rhisc=egrs(i)/max(1.e-24,betafp(i,j))
+ rhisc=egrs(i)/max(1.e-22*iL_mks2cgs,betafp(i,j))
els=max(eglsmn,min(bcrrd(i),rhisc))
urmse(i)=2.86*egc*egrs(i)*els
c
c --- ------- Zonal eddy phase speed minus zonal barotropic velocity
c --- ------- with a lower bound of -20 cm s-1.
- cpse(i)=max(-20.,-betafp(i,j)*bcrrd(i)**2)
+ cpse(i)=max(cpsemin,-betafp(i,j)*bcrrd(i)**2)
c
endif
c
@@ -1606,7 +1625,8 @@ subroutine difest_lateral_hyb(m,n,mm,nn,k1m,k1n)
c --- --------- zonal velocity minus eddy phase speed and absolute value
c --- --------- of RMS eddy velocity is set to -20 cm s-1 and 5 cm s-1,
c --- --------- respectively.
- esfac=1./(1.+4.*(umnsc/max(5.,abs(urmse(i))))**2)
+ esfac=1./
+ . (1.+4.*(umnsc/max(urmsemin,abs(urmse(i))))**2)
c
else
esfac=1.
@@ -1637,7 +1657,8 @@ subroutine difest_lateral_hyb(m,n,mm,nn,k1m,k1n)
. -cpse(i)
c
c --- ----------- Eddy mixing suppresion factor.
- esfac=1./(1.+4.*(umnsc/max(5.,abs(urmse(i))))**2)
+ esfac=1./
+ . (1.+4.*(umnsc/max(urmsemin,abs(urmse(i))))**2)
c
else
esfac=1.
@@ -1718,7 +1739,7 @@ subroutine difest_lateral_iso(m,n,mm,nn,k1m,k1n)
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
if (k.ge.kfpla(i,j,n)) then
- if (p(i,j,kk+1)-p(i,j,k+1).lt.epsil) then
+ if (p(i,j,kk+1)-p(i,j,k+1).lt.epsilp) then
plo=p(i,j,kk+1)
else
plo=.5*(p(i,j,k)+p(i,j,k+1))
@@ -1915,7 +1936,7 @@ subroutine difest_lateral_iso(m,n,mm,nn,k1m,k1n)
. kmax(i,j)-kfil(i,j).ge.1) then
c
c --- --------- Rhines scale.
- rhisc=egr(i,k)/max(1.e-24,betafp(i,j))
+ rhisc=egr(i,k)/max(1.e-22*iL_mks2cgs,betafp(i,j))
c
c --- --------- Eddy length scale.
els=max(eglsmn,min(bcrrd(i),rhisc))
@@ -1951,13 +1972,13 @@ subroutine difest_lateral_iso(m,n,mm,nn,k1m,k1n)
c --- ------- RMS eddy velocity estimated from K = Gamma*u_rms*L, where
c --- ------- a mixing efficiency of Gamma = 0.35 is used (Klocker and
c --- ------- Abernathey, 2014).
- rhisc=egrs(i)/max(1.e-24,betafp(i,j))
+ rhisc=egrs(i)/max(1.e-22*iL_mks2cgs,betafp(i,j))
els=max(eglsmn,min(bcrrd(i),rhisc))
urmse(i)=2.86*egc*egrs(i)*els
c
c --- ------- Zonal eddy phase speed minus zonal barotropic velocity
c --- ------- with a lower bound of -20 cm s-1.
- cpse(i)=max(-20.,-betafp(i,j)*bcrrd(i)**2)
+ cpse(i)=max(cpsemin,-betafp(i,j)*bcrrd(i)**2)
c
endif
c
@@ -2012,7 +2033,8 @@ subroutine difest_lateral_iso(m,n,mm,nn,k1m,k1n)
c --- --------- zonal velocity minus eddy phase speed and absolute value
c --- --------- of RMS eddy velocity is set to -20 cm s-1 and 5 cm s-1,
c --- --------- respectively.
- esfac=1./(1.+4.*(umnsc/max(5.,abs(urmse(i))))**2)
+ esfac=1./
+ . (1.+4.*(umnsc/max(urmsemin,abs(urmse(i))))**2)
c
else
esfac=1.
@@ -2047,7 +2069,8 @@ subroutine difest_lateral_iso(m,n,mm,nn,k1m,k1n)
. -cpse(i)
c
c --- ----------- Eddy mixing suppresion factor.
- esfac=1./(1.+4.*(umnsc/max(5.,abs(urmse(i))))**2)
+ esfac=
+ . 1./(1.+4.*(umnsc/max(urmsemin,abs(urmse(i))))**2)
c
else
esfac=1.
@@ -2141,7 +2164,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
c
c --- ------- Brunt-Vaisala frequency squared
bvfsq(i,k)=g*g*max(drhomn,drhol(i,j,k))
- . /max(epsil,dp(i,j,kn))
+ . /max(epsilp,dp(i,j,kn))
c
c --- ------- Brunt-Vaisala frequency
bvf(i,k)=sqrt(bvfsq(i,k))
@@ -2152,7 +2175,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
h=max(onem,dp(i,j,kn))*alpha0/g
c h=max(onem*1e-8,dp(i,j,kn))*alpha0/g
c h=max(onemm,dp(i,j,kn))*alpha0/g
- Shear2(i,j,k)=max(1.e-9,du2l(i,j,k))/(h*h)
+ Shear2(i,j,k)=max(1.e-13*A_mks2cgs,du2l(i,j,k))/(h*h)
Prod(i,j,k)=difdia(i,j,k)*Pr_t*Shear2(i,j,k)
else
Buoy(i,j,k)=0.
@@ -2307,7 +2330,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
c
c --- ------- Penetration of surface TKE below mixed layer.
if (tkepf.gt.0.) then
- if (dp(i,j,kn).lt.epsil) then
+ if (dp(i,j,kn).lt.epsilp) then
q=exp(-p(i,j,k)/tkepls)
else
q=tkepls*(exp(-p(i,j,k )/tkepls)
@@ -2319,7 +2342,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
c
c --- ------- Set TKE and GLS to prescribed minimum values in surface
c --- ------- mixed layers and thin layers
- if (dp(i,j,kn).lt.epsil) then
+ if (dp(i,j,kn).lt.epsilp) then
trc(i,j,kn,itrtke)=tke_min
trc(i,j,kn,itrgls)=gls_psi_min
endif
@@ -2336,7 +2359,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
trc(i,j,kn,itrgls)=max(gls_psi_min,
. (gls_cmu0**(gls_p-2.*gls_m))
. *(ust**(2.*gls_m))
- . *(kappa*1.e2)**gls_n)
+ . *(kappa*L_mks2cgs)**gls_n)
# endif
endif
c
@@ -2408,7 +2431,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
if (tdmflg.eq.1) then
q=.5*(tanh(4.*(abs(plat(i,j))-tdclat)/tddlat-2.)+1.)
q=(1.-q)*tdmls0+q*tdmls1
- if (dp(i,j,kn).lt.epsil) then
+ if (dp(i,j,kn).lt.epsilp) then
vsf=exp(p(i,j,k)/q)/(q*(exp(p(i,j,kk+1)/q)-1.))
else
vsf=(exp(p(i,j,k+1)/q)-exp(p(i,j,k)/q))
@@ -2470,7 +2493,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
if (k.lt.kfil(i,j)) then
if (k.gt.2.and.kfil(i,j).le.kk.and.
- . p(i,j,min(kk,kfil(i,j)))-p(i,j,3).gt.epsil) then
+ . p(i,j,min(kk,kfil(i,j)))-p(i,j,3).gt.epsilp) then
q=.5*(p(i,j,k+1)+p(i,j,k))
difdia(i,j,k)=((q-p(i,j,3))*dfddsl(i)
. +(p(i,j,kfil(i,j))-q)*dfddsu(i))
@@ -2491,7 +2514,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
if (k.le.kmax(i,j).and.kmax(i,j)-kfil(i,j).ge.1) then
q=niwls
- if (k.eq.2.or.dp(i,j,kn).lt.epsil) then
+ if (k.eq.2.or.dp(i,j,kn).lt.epsilp) then
vsf=exp((p(i,j,3)-p(i,j,k+1))/q)
. /(q*(1.-exp((p(i,j,3)-p(i,j,kk+1))/q)))
else
@@ -2519,7 +2542,7 @@ subroutine difest_vertical_iso(m,n,mm,nn,k1m,k1n)
. -buoyfl(i,j,1)))
c
c --- --- Mixed layer thickness
- h=(p(i,j,3)-p(i,j,1))/onecm
+ h=(p(i,j,3)-p(i,j,1))/(onem*iL_mks2cgs)
c
c --- --- Dimensionless vertical coordinate in the boundary layer
sg=(p(i,j,2)-p(i,j,1))/(p(i,j,3)-p(i,j,1))
diff --git a/phy/mod_diffusion.F90 b/phy/mod_diffusion.F90
index 638a5f6d..4f82d958 100644
--- a/phy/mod_diffusion.F90
+++ b/phy/mod_diffusion.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2020-2022 Mats Bentsen
+! Copyright (C) 2020-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,7 +24,7 @@ module mod_diffusion
use mod_types, only: r8
use mod_config, only: inst_suffix
- use mod_constants, only: spval, epsil
+ use mod_constants, only: spval, epsilk
use mod_xc
implicit none
@@ -320,9 +320,9 @@ subroutine inivar_diffusion
enddo
do k = 1, kk+1
do i = 1 - nbdy, ii + nbdy
- Kvisc_m(i, j, k) = epsil
- Kdiff_t(i, j, k) = epsil
- Kdiff_s(i, j, k) = epsil
+ Kvisc_m(i, j, k) = epsilk
+ Kdiff_t(i, j, k) = epsilk
+ Kdiff_s(i, j, k) = epsilk
enddo
enddo
enddo
diff --git a/phy/mod_eddtra.F90 b/phy/mod_eddtra.F90
index 87df3bda..1490cf92 100644
--- a/phy/mod_eddtra.F90
+++ b/phy/mod_eddtra.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2022 Mats Bentsen
+! Copyright (C) 2015-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,7 +24,8 @@ module mod_eddtra
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: g, alpha0, epsil, onecm, onemm
+ use mod_constants, only: g, alpha0, rho0, epsilp, onem, onecm, onemm, &
+ L_mks2cgs
use mod_time, only: delt1
use mod_xc
use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml, cntiso_hybrid
@@ -39,6 +40,9 @@ module mod_eddtra
implicit none
+ real(r8), parameter :: &
+ iL_mks2cgs = 1./L_mks2cgs
+
private
public :: eddtra
@@ -149,12 +153,10 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
real(r8), dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: ptu, ptv
real(r8), dimension(kdm+1) :: upsilon, mfl
real(r8), dimension(kdm) :: dlm, dlp
- real(r8) :: rho0, q, et2mf, kappa, fhi, flo
+ real(r8) :: q, et2mf, kappa, fhi, flo
integer :: i, j, k, l, km, kn, kintr, kmax, kmin, niter, kdir
logical :: changed
- rho0 = 1._r8/alpha0
-
call xctilr(difint, 1, kk, 2, 2, halo_ps)
call xctilr(pbu, 1, 2, 2, 2, halo_us)
call xctilr(pbv, 1, 2, 2, 2, halo_vs)
@@ -205,7 +207,8 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
kmax = 1
do k = 3, kk
kn = k + nn
- if (dp(i - 1, j, kn) > epsil .or. dp(i, j, kn) > epsil) kmax = k
+ if (dp(i - 1, j, kn) > epsilp .or. dp(i, j, kn) > epsilp) &
+ kmax = k
enddo
! ------------------------------------------------------------------
@@ -248,7 +251,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
temp(i - 1, j, kn), saln(i - 1, j, kn)) < &
rho(p(i , j, 3), &
temp(i , j, km), saln(i , j, km)) .or. &
- dp(i - 1, j, kn) < epsil)
+ dp(i - 1, j, kn) < epsilp)
kintr = kintr + 1
if (kintr == kmax + 1) exit
kn = kintr + nn
@@ -290,7 +293,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
temp(i , j, kn), saln(i , j, kn)) < &
rho(p(i - 1, j, 3), &
temp(i - 1, j, km), saln(i - 1, j, km)) .or. &
- dp(i , j, kn) < epsil)
+ dp(i , j, kn) < epsilp)
kintr = kintr + 1
if (kintr == kmax + 1) exit
kn = kintr + nn
@@ -412,10 +415,10 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
mfl(kmin + 1) = min(fhi, max(flo, mfl(kmin + 1)))
do k = kmin + 1, kmax - 1
if (mfl(k + 1) - mfl(k) > &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)) then
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)) then
mfl(k + 1) = mfl(k) + fface*dlm(k)*scp2(i - 1, j)
elseif (mfl(k + 1) - mfl(k) < &
- - ffac*max(epsil, dlp(k))*scp2(i , j)) then
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)) then
mfl(k + 1) = mfl(k) - fface*dlp(k)*scp2(i , j)
else
exit
@@ -462,11 +465,11 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
! difference between lower and upper interface is beyond the
! floating point accuracy limitation.
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scu2(i, j), &
+ 1.e-14_r8*max(epsilp*scu2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
if (mfl(k + 1) - mfl(k) > &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)) then
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i - 1, j, k). Limit the
! dominating interface flux.
@@ -488,7 +491,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
endif
changed = .true.
elseif (mfl(k + 1) - mfl(k) < &
- - ffac*max(epsil, dlp(k))*scp2(i , j)) then
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i, j, k). Limit the
! dominating interface flux.
@@ -522,7 +525,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
k = kmin
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scu2(i, j), &
+ 1.e-14_r8*max(epsilp*scu2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
umfltd(i, j, 2 + mm) = mfl(k + 1) - mfl(k)
umfltd(i, j, 1 + mm) = umfltd(i, j, 2 + mm) &
@@ -537,25 +540,25 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
do k = kintr, kmax
km = k + mm
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scu2(i, j), &
+ 1.e-14_r8*max(epsilp*scu2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
umfltd(i, j, km) = mfl(k + 1) - mfl(k)
else
umfltd(i, j, km) = 0._r8
endif
if (umfltd(i, j, km) > &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)) then
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)) then
write(lp,*) 'eddtra_gm_isopyc_bulkml u >', &
i + i0, j + j0, k, umfltd(i, j, km), &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)
call xchalt('(eddtra_gm_isopyc_bulkml)')
stop '(eddtra_gm_isopyc_bulkml)'
endif
if (umfltd(i, j, km) < &
- - ffac*max(epsil, dlp(k))*scp2(i , j)) then
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)) then
write(lp,*) 'eddtra_gm_isopyc_bulkml u <', &
i + i0, j + j0, k, umfltd(i, j, km), &
- - ffac*max(epsil, dlp(k))*scp2(i , j)
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)
call xchalt('(eddtra_gm_isopyc_bulkml)')
stop '(eddtra_gm_isopyc_bulkml)'
endif
@@ -591,7 +594,8 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
kmax = 1
do k = 3, kk
kn = k + nn
- if (dp(i, j - 1, kn) > epsil .or. dp(i, j, kn) > epsil) kmax = k
+ if (dp(i, j - 1, kn) > epsilp .or. dp(i, j, kn) > epsilp) &
+ kmax = k
enddo
! ------------------------------------------------------------------
@@ -634,7 +638,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
temp(i, j - 1, kn), saln(i, j - 1, kn)) < &
rho(p(i, j , 3), &
temp(i, j , km), saln(i, j , km)) .or. &
- dp(i, j - 1, kn) < epsil)
+ dp(i, j - 1, kn) < epsilp)
kintr = kintr + 1
if (kintr == kmax + 1) exit
kn = kintr + nn
@@ -676,7 +680,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
temp(i, j , kn), saln(i, j , kn)) < &
rho(p(i, j - 1, 3), &
temp(i, j - 1, km), saln(i, j - 1, km)) .or. &
- dp(i, j , kn) < epsil)
+ dp(i, j , kn) < epsilp)
kintr = kintr + 1
if (kintr == kmax + 1) exit
kn = kintr + nn
@@ -798,10 +802,10 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
mfl(kmin + 1) = min(fhi, max(flo, mfl(kmin + 1)))
do k = kmin + 1, kmax - 1
if (mfl(k + 1) - mfl(k) > &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)) then
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)) then
mfl(k + 1) = mfl(k) + fface*dlm(k)*scp2(i, j - 1)
elseif (mfl(k + 1) - mfl(k) < &
- - ffac*max(epsil, dlp(k))*scp2(i, j )) then
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )) then
mfl(k + 1) = mfl(k) - fface*dlp(k)*scp2(i, j )
else
exit
@@ -848,11 +852,11 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
! difference between lower and upper interface is beyond the
! floating point accuracy limitation.
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scv2(i, j), &
+ 1.e-14_r8*max(epsilp*scv2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
if (mfl(k + 1) - mfl(k) > &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)) then
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i, j - 1, k). Limit the
! dominating interface flux.
@@ -874,7 +878,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
endif
changed = .true.
elseif (mfl(k + 1) - mfl(k) < &
- - ffac*max(epsil, dlp(k))*scp2(i, j )) then
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i, j, k). Limit the
! dominating interface flux.
@@ -908,7 +912,7 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
k = kmin
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scv2(i, j), &
+ 1.e-14_r8*max(epsilp*scv2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
vmfltd(i, j, 2 + mm) = mfl(k + 1) - mfl(k)
vmfltd(i, j, 1 + mm) = vmfltd(i, j, 2 + mm) &
@@ -923,25 +927,25 @@ subroutine eddtra_gm_isopyc_bulkml(m, n, mm, nn, k1m, k1n)
do k = kintr, kmax
km = k + mm
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scv2(i, j), &
+ 1.e-14_r8*max(epsilp*scv2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
vmfltd(i, j, km) = mfl(k + 1) - mfl(k)
else
vmfltd(i, j, km) = 0._r8
endif
if (vmfltd(i, j, km) > &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)) then
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)) then
write(lp,*) 'eddtra_gm_isopyc_bulkml v >', &
i + i0, j + j0, k, vmfltd(i, j, km), &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)
call xchalt('(eddtra_gm_isopyc_bulkml)')
stop '(eddtra_gm_isopyc_bulkml)'
endif
if (vmfltd(i, j, km) < &
- - ffac*max(epsil, dlp(k))*scp2(i, j )) then
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )) then
write(lp,*) 'eddtra_gm_isopyc_bulkml v <', &
i + i0, j + j0, k, vmfltd(i, j, km), &
- - ffac*max(epsil, dlp(k))*scp2(i, j )
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )
call xchalt('(eddtra_gm_isopyc_bulkml)')
stop '(eddtra_gm_isopyc_bulkml)'
endif
@@ -971,12 +975,10 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
real(r8), dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: ptu, ptv
real(r8), dimension(kdm+1) :: mfl
real(r8), dimension(kdm) :: puv, dlm, dlp
- real(r8) :: rho0, q, et2mf, mlp, kappa
+ real(r8) :: q, et2mf, mlp, kappa
integer :: i, j, k, l, km, kn, kmax, kml, niter, kdir
logical :: changed
- rho0 = 1._r8/alpha0
-
call xctilr(difint, 1, kk, 2, 2, halo_ps)
call xctilr(pbu, 1, 2, 2, 2, halo_us)
call xctilr(pbv, 1, 2, 2, 2, halo_vs)
@@ -1028,12 +1030,13 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
do k = 2, kk
kn = k + nn
puv(k) = puv(k - 1) + dpu(i, j, kn - 1)
- if (dp(i - 1, j, kn) > epsil .or. dp(i, j, kn) > epsil) kmax = k
+ if (dp(i - 1, j, kn) > epsilp .or. dp(i, j, kn) > epsilp) &
+ kmax = k
enddo
! Compute the eddy induced mass flux at layer interfaces below the
! mixed layer.
- mlp = .5_r8*(mlts(i - 1, j) + mlts(i, j))*onecm
+ mlp = .5_r8*(mlts(i - 1, j) + mlts(i, j))*(onem*iL_mks2cgs)
kml = kmax + 1
mfl(kmax + 1) = 0._r8
do k = kmax, 2, -1
@@ -1106,11 +1109,11 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
! difference between lower and upper interface is beyond the
! floating point accuracy limitation.
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scu2(i, j), &
+ 1.e-14_r8*max(epsilp*scu2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
if (mfl(k + 1) - mfl(k) > &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)) then
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i - 1, j, k). Limit the
! dominating interface flux.
@@ -1132,7 +1135,7 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
endif
changed = .true.
elseif (mfl(k + 1) - mfl(k) < &
- - ffac*max(epsil, dlp(k))*scp2(i , j)) then
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i, j, k). Limit the
! dominating interface flux.
@@ -1167,25 +1170,25 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
do k = 1, kmax
km = k + mm
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scu2(i, j), &
+ 1.e-14_r8*max(epsilp*scu2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
umfltd(i, j, km) = mfl(k + 1) - mfl(k)
else
umfltd(i, j, km) = 0._r8
endif
if (umfltd(i, j, km) > &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)) then
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)) then
write(lp,*) 'eddtra_gm_cntiso_hybrid u >', &
i + i0, j + j0, k, umfltd(i, j, km), &
- ffac*max(epsil, dlm(k))*scp2(i - 1, j)
+ ffac*max(epsilp, dlm(k))*scp2(i - 1, j)
call xchalt('(eddtra_gm_cntiso_hybrid)')
stop '(eddtra_gm_cntiso_hybrid)'
endif
if (umfltd(i, j, km) < &
- - ffac*max(epsil, dlp(k))*scp2(i , j)) then
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)) then
write(lp,*) 'eddtra_gm_cntiso_hybrid u <', &
i + i0, j + j0, k, umfltd(i, j, km), &
- - ffac*max(epsil, dlp(k))*scp2(i , j)
+ - ffac*max(epsilp, dlp(k))*scp2(i , j)
call xchalt('(eddtra_gm_cntiso_hybrid)')
stop '(eddtra_gm_cntiso_hybrid)'
endif
@@ -1222,12 +1225,13 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
do k = 2, kk
kn = k + nn
puv(k) = puv(k - 1) + dpv(i, j, kn - 1)
- if (dp(i, j - 1, kn) > epsil .or. dp(i, j, kn) > epsil) kmax = k
+ if (dp(i, j - 1, kn) > epsilp .or. dp(i, j, kn) > epsilp) &
+ kmax = k
enddo
! Compute the eddy induced mass flux at layer interfaces below the
! mixed layer.
- mlp = .5_r8*(mlts(i, j - 1) + mlts(i, j))*onecm
+ mlp = .5_r8*(mlts(i, j - 1) + mlts(i, j))*(onem*iL_mks2cgs)
kml = kmax + 1
mfl(kmax + 1) = 0._r8
do k = kmax, 2, -1
@@ -1300,11 +1304,11 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
! difference between lower and upper interface is beyond the
! floating point accuracy limitation.
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scv2(i, j), &
+ 1.e-14_r8*max(epsilp*scv2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
if (mfl(k + 1) - mfl(k) > &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)) then
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i, j - 1, k). Limit the
! dominating interface flux.
@@ -1326,7 +1330,7 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
endif
changed = .true.
elseif (mfl(k + 1) - mfl(k) < &
- - ffac*max(epsil, dlp(k))*scp2(i, j )) then
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )) then
! In this case, the mass fluxes are removing too much
! mass from the grid cell at (i, j, k). Limit the
! dominating interface flux.
@@ -1361,25 +1365,25 @@ subroutine eddtra_gm_cntiso_hybrid(m, n, mm, nn, k1m, k1n)
do k = 1, kmax
km = k + mm
if (abs(mfl(k + 1) - mfl(k)) > &
- 1.e-14_r8*max(epsil*scv2(i, j), &
+ 1.e-14_r8*max(epsilp*scv2(i, j), &
abs(mfl(k + 1) + mfl(k)))) then
vmfltd(i, j, km) = mfl(k + 1) - mfl(k)
else
vmfltd(i, j, km) = 0._r8
endif
if (vmfltd(i, j, km) > &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)) then
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)) then
write(lp,*) 'eddtra_gm_cntiso_hybrid v >', &
i + i0, j + j0, k, vmfltd(i, j, km), &
- ffac*max(epsil, dlm(k))*scp2(i, j - 1)
+ ffac*max(epsilp, dlm(k))*scp2(i, j - 1)
call xchalt('(eddtra_gm_cntiso_hybrid)')
stop '(eddtra_gm_cntiso_hybrid)'
endif
if (vmfltd(i, j, km) < &
- - ffac*max(epsil, dlp(k))*scp2(i, j )) then
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )) then
write(lp,*) 'eddtra_gm_cntiso_hybrid v <', &
i + i0, j + j0, k, vmfltd(i, j, km), &
- - ffac*max(epsil, dlp(k))*scp2(i, j )
+ - ffac*max(epsilp, dlp(k))*scp2(i, j )
call xchalt('(eddtra_gm_cntiso_hybrid)')
stop '(eddtra_gm_cntiso_hybrid)'
endif
diff --git a/phy/mod_eos.F90 b/phy/mod_eos.F90
index 49bbef46..7226abff 100644
--- a/phy/mod_eos.F90
+++ b/phy/mod_eos.F90
@@ -24,6 +24,7 @@ module mod_eos
! ------------------------------------------------------------------------------
use mod_types, only: r8
+ use mod_constants, only: alpha0
use mod_config, only: expcnf
use mod_xc, only: mnproc, lp, xcstop
@@ -32,6 +33,27 @@ module mod_eos
private
! Coefficients for the functional fit of in situ density.
+#ifdef MKS
+ real(r8), parameter :: &
+ a11 = 9.9985372432159340e+02_r8, &
+ a12 = 1.0380621928183473e+01_r8, &
+ a13 = 1.7073577195684715e+00_r8, &
+ a14 = -3.6570490496333680e-02_r8, &
+ a15 = -7.3677944503527477e-03_r8, &
+ a16 = -3.5529175999643348e-03_r8, &
+ b11 = 1.7083494994335439e-06_r8, &
+ b12 = 7.1567921402953455e-09_r8, &
+ b13 = 1.2821026080049485e-09_r8, &
+ a21 = 1.0_r8 , &
+ a22 = 1.0316374535350838e-02_r8, &
+ a23 = 8.9521792365142522e-04_r8, &
+ a24 = -2.8438341552142710e-05_r8, &
+ a25 = -1.1887778959461776e-05_r8, &
+ a26 = -4.0163964812921489e-06_r8, &
+ b21 = 1.1995545126831476e-09_r8, &
+ b22 = 5.5234008384648383e-12_r8, &
+ b23 = 8.4310335919950873e-13_r8
+#else
real(r8), parameter :: &
a11 = 9.9985372432159340e-01_r8, &
a12 = 1.0380621928183473e-02_r8, &
@@ -51,6 +73,7 @@ module mod_eos
b21 = 1.1995545126831476e-10_r8, &
b22 = 5.5234008384648383e-13_r8, &
b23 = 8.4310335919950873e-14_r8
+#endif
! Reference pressure [g cm-1 s-2].
real(r8) :: pref
@@ -106,12 +129,12 @@ subroutine inieos
ap24 = a24
ap25 = a25
ap26 = a26
- ap11 = a11 + b11*pref - ap21
- ap12 = a12 + b12*pref - ap22
- ap13 = a13 + b13*pref - ap23
- ap14 = a14 - ap24
- ap15 = a15 - ap25
- ap16 = a16 - ap26
+ ap11 = a11 + b11*pref - ap21/alpha0
+ ap12 = a12 + b12*pref - ap22/alpha0
+ ap13 = a13 + b13*pref - ap23/alpha0
+ ap14 = a14 - ap24/alpha0
+ ap15 = a15 - ap25/alpha0
+ ap16 = a16 - ap26/alpha0
ap210 = a21
ap220 = a22
@@ -119,12 +142,12 @@ subroutine inieos
ap240 = a24
ap250 = a25
ap260 = a26
- ap110 = a11 - ap210
- ap120 = a12 - ap220
- ap130 = a13 - ap230
- ap140 = a14 - ap240
- ap150 = a15 - ap250
- ap160 = a16 - ap260
+ ap110 = a11 - ap210/alpha0
+ ap120 = a12 - ap220/alpha0
+ ap130 = a13 - ap230/alpha0
+ ap140 = a14 - ap240/alpha0
+ ap150 = a15 - ap250/alpha0
+ ap160 = a16 - ap260/alpha0
! Coefficients for freezing temperature.
select case (trim(expcnf))
diff --git a/phy/mod_inicon.F b/phy/mod_inicon.F
index fdfebd50..70e0c371 100644
--- a/phy/mod_inicon.F
+++ b/phy/mod_inicon.F
@@ -27,7 +27,8 @@ module mod_inicon
c
use mod_types, only: r8
use mod_config, only: expcnf
- use mod_constants, only: g, epsil, onem
+ use mod_constants, only: g, epsilp, onem,
+ . L_mks2cgs, M_mks2cgs, P_mks2cgs
use mod_time, only: nstep, delt1, dlt
use mod_xc
use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml,
@@ -98,7 +99,7 @@ function getpl(th,s,phiu,phil,pup) result(plo)
c --- improve the accuracy of the pressure interface by an
c --- iterative procedure
q=1._r8
- do while (abs(q).gt.1.e-4_r8)
+ do while (abs(q).gt.1.e-5_r8*P_mks2cgs)
call delphi(pup,plo,th,s,dphi,alpu,alpl)
q=(phil-phiu-dphi)/alpl
plo=plo-q
@@ -121,6 +122,9 @@ subroutine ictsz_file
real dsig,a0,a1,a2
integer, dimension(3) :: start,count
integer i,j,kdmic,k,l,status,ncid,dimid,varid,kb
+ real iM_mks2cgs
+c
+ iM_mks2cgs = 1.0 / M_mks2cgs
c
if (mnproc.eq.1) then
write (lp,'(2a)') ' reading initial condition from ',
@@ -344,7 +348,7 @@ subroutine ictsz_file
do j=1,jj
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
-c z(i,j,1)=z(i,j,1)*1.e2
+c z(i,j,1)=z(i,j,1)*L_mks2cgs
z(i,j,1)=0.
enddo
enddo
@@ -355,7 +359,8 @@ subroutine ictsz_file
do k=1,kdmic
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- z(i,j,k+1)=min(depths(i,j)*1.e2,z(i,j,k)+dz(i,j,k)*1.e2)
+ z(i,j,k+1)=min(depths(i,j)*L_mks2cgs,
+ . z(i,j,k)+dz(i,j,k)*L_mks2cgs)
enddo
enddo
enddo
@@ -369,20 +374,20 @@ subroutine ictsz_file
do k=2,kk
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- if (z(i,j,kk+1)-z(i,j,k).lt.1.e-4)
- . z(i,j,k)=depths(i,j)*1.e2
+ if (z(i,j,kk+1)-z(i,j,k).lt.1.e-6*L_mks2cgs)
+ . z(i,j,k)=depths(i,j)*L_mks2cgs
enddo
enddo
enddo
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- z(i,j,kk+1)=depths(i,j)*1.e2
+ z(i,j,kk+1)=depths(i,j)*L_mks2cgs
enddo
enddo
do k=1,kk
do l=1,isp(j)
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
- sigmar(i,j,k)=sigmar(i,j,k)*1.e-3
+ sigmar(i,j,k)=sigmar(i,j,k)*iM_mks2cgs
enddo
enddo
enddo
@@ -876,7 +881,7 @@ subroutine inicon
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
k=3
dps=0.
- do while (dp(i,j,k).lt.epsil)
+ do while (dp(i,j,k).lt.epsilp)
dps=dps+dp(i,j,k)
dp(i,j,k)=0.
dp(i,j,k+kk)=0.
@@ -920,7 +925,8 @@ subroutine inicon
j=jtest
write (lp,103) nstep,i0+i,j0+j,
. ' init.profile temp saln dens thkns dpth',
- . (k,temp(i,j,k),saln(i,j,k),1000.*sig(temp(i,j,k),saln(i,j,k)),
+ . (k,temp(i,j,k),saln(i,j,k),
+ . M_mks2cgs*sig(temp(i,j,k),saln(i,j,k)),
. dp(i,j,k)/onem,p(i,j,k+1)/onem,k=1,kk)
103 format (i9,2i5,a/(28x,i3,3f8.2,2f8.1))
endif
diff --git a/phy/mod_momtum.F b/phy/mod_momtum.F
index ad98f07a..9eea25bf 100644
--- a/phy/mod_momtum.F
+++ b/phy/mod_momtum.F
@@ -1,6 +1,6 @@
! ------------------------------------------------------------------------------
! Copyright (C) 2000 HYCOM Consortium and contributors
-! Copyright (C) 2001-2020 Mats Bentsen, Lars Inge Enstad, Mehmet Ilicak
+! Copyright (C) 2001-2022 Mats Bentsen, Lars Inge Enstad, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -26,7 +26,8 @@ module mod_momtum
c --- ------------------------------------------------------------------
c
use mod_types, only: r8
- use mod_constants, only: g, alpha0, epsil, spval, onem, onemm
+ use mod_constants, only: g, alpha0, epsilp, epsilpl, spval,
+ . onem, onemm
use mod_time, only: delt1, dlt
use mod_xc
use mod_grid, only: scqx, scqy, scpx, scpy, scux, scuy,
@@ -276,14 +277,14 @@ subroutine momtum(m,n,mm,nn,k1m,k1n)
c
do i=max(0,ifp(j,l)),min(ii,ilp(j,l))
ubot=(ubflxs_p(i ,j,n)
- . /max(epsil,pbu(i ,j,n)*scuy(i ,j))
+ . /max(epsilpl,pbu(i ,j,n)*scuy(i ,j))
. +ubflxs_p(i+1,j,n)
- . /max(epsil,pbu(i+1,j,n)*scuy(i+1,j)))*tsfac
+ . /max(epsilpl,pbu(i+1,j,n)*scuy(i+1,j)))*tsfac
. +util1(i,j)/thkbop
vbot=(vbflxs_p(i,j ,n)
- . /max(epsil,pbv(i,j ,n)*scvx(i,j ))
+ . /max(epsilpl,pbv(i,j ,n)*scvx(i,j ))
. +vbflxs_p(i,j+1,n)
- . /max(epsil,pbv(i,j+1,n)*scvx(i,j+1)))*tsfac
+ . /max(epsilpl,pbv(i,j+1,n)*scvx(i,j+1)))*tsfac
. +util2(i,j)/thkbop
ubbl=.5*sqrt(ubot*ubot+vbot*vbot)
q=cb*(ubbl+cbar)
@@ -445,9 +446,9 @@ subroutine momtum(m,n,mm,nn,k1m,k1n)
do l=1,isu(j)
do i=max(0,ifu(j,l)),min(ii+2,ilu(j,l))
wgtja(i,j)=max(0.,min(1.,(pu(i,j,k+1)-pbu(i,j-1,m))
- . /max(pu(i,j,k+1)-pu(i,j,k),epsil)))
+ . /max(pu(i,j,k+1)-pu(i,j,k),epsilp)))
wgtjb(i,j)=max(0.,min(1.,(pu(i,j,k+1)-pbu(i,j+1,m))
- . /max(pu(i,j,k+1)-pu(i,j,k),epsil)))
+ . /max(pu(i,j,k+1)-pu(i,j,k),epsilp)))
uja(i,j)=(1.-wgtja(i,j))*utotn(i,j-1)
. +wgtja(i,j)*slip*utotn(i,j)
ujb(i,j)=(1.-wgtjb(i,j))*utotn(i,j+1)
@@ -464,9 +465,9 @@ subroutine momtum(m,n,mm,nn,k1m,k1n)
do l=1,isv(j)
do i=max(-1,ifv(j,l)),min(ii+2,ilv(j,l))
wgtia(i,j)=max(0.,min(1.,(pv(i,j,k+1)-pbv(i-1,j,m))
- . /max(pv(i,j,k+1)-pv(i,j,k),epsil)))
+ . /max(pv(i,j,k+1)-pv(i,j,k),epsilp)))
wgtib(i,j)=max(0.,min(1.,(pv(i,j,k+1)-pbv(i+1,j,m))
- . /max(pv(i,j,k+1)-pv(i,j,k),epsil)))
+ . /max(pv(i,j,k+1)-pv(i,j,k),epsilp)))
via(i,j)=(1.-wgtia(i,j))*vtotn(i-1,j)
. +wgtia(i,j)*slip*vtotn(i,j)
vib(i,j)=(1.-wgtib(i,j))*vtotn(i+1,j)
diff --git a/phy/mod_mxlayr.F b/phy/mod_mxlayr.F
index c81abb6a..fa5ebffb 100644
--- a/phy/mod_mxlayr.F
+++ b/phy/mod_mxlayr.F
@@ -25,8 +25,9 @@ module mod_mxlayr
c --- ------------------------------------------------------------------
c
use mod_types, only: r8
- use mod_constants, only: g, spcifh, alpha0, epsil, spval, onem,
- . tencm, onecm, onemm, onemu
+ use mod_constants, only: g, spcifh, alpha0, epsilp, spval, onem,
+ . tencm, onecm, onemm, onemu,
+ . L_mks2cgs, R_mks2cgs
use mod_time, only: delt1
use mod_xc
use mod_vcoord, only: sigmar
@@ -89,6 +90,11 @@ module mod_mxlayr
. mtkeke, ! Mixed layer TKE tendency related to kin.
! energy change [cm3 s-3].
. pbrnda ! Brine plume pressure depth [g cm-1 s-2].
+c
+ real(r8), parameter ::
+ . iL_mks2cgs = 1./L_mks2cgs,
+ . A_cgs2mks = 1./(L_mks2cgs*L_mks2cgs),
+ . V_mks2cgs = L_mks2cgs**3
c
public :: rm0,rm5,ce,mlrttp,mltmin,
. mtkeus,mtkeni,mtkebf,mtkers,mtkepe,mtkeke,pbrnda,
@@ -163,7 +169,8 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c --- of TKE balance [].
real kappa,mu,ustmin,mldjmp
integer maxitr
- parameter (kappa=.4,mu=2.,ustmin=.1,mldjmp=1.e-6,maxitr=20)
+ parameter (kappa=.4,mu=2.,ustmin=.001*L_mks2cgs,
+ . mldjmp=1.e-3*R_mks2cgs,maxitr=20)
c
c --- Parameters for the parameterization of restratification by mixed
c --- layer eddies by Fox-Kemper et al. (2008):
@@ -175,12 +182,12 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c --- ci - constant that appears when integrating the shape
c --- function over the mixed layer depth [].
real rtau,cori20,rlf,ci,slbg0
- parameter (rtau=1./86400.,cori20=4.9745e-5,rlf=1./5.e5,
- . ci=44./63.,slbg0=0.)
+ parameter (rtau=1./86400.,cori20=4.9745e-5,
+ . rlf=1./(5.e3*L_mks2cgs),ci=44./63.,slbg0=0.)
c
c --- Parameters for brine plume parameterization:
c --- bpdrho - density contrast between surface and brine plume depth
-c --- [g/cm/s**2].
+c --- [g/cm**3].
c --- bpmndp - minimum distribution thickness of salt from sea-ice
c --- freezing [g/cm/s**2].
c --- bpmxdp - maximum distribution depth below the mixed layer base
@@ -190,8 +197,8 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c --- dsgmnr - minimum ratio of linearized density jump to target
c --- density jump across a layer interface [].
real bpdrho,bpmndp,bpmxdp,bpdpmn,dsgmnr
- parameter (bpdrho=.4e-3,bpmndp=10.*98060.,bpmxdp=500.*98060.,
- . bpdpmn=1.*98060.,dsgmnr=.1)
+ parameter (bpdrho=.4*R_mks2cgs,bpmndp=10.*onem,
+ . bpmxdp=500.*onem,bpdpmn=1.*onem,dsgmnr=.1)
c
c --- ------------------------------------------------------------------
c --- Resolve type of mixed layer restratification time scale.
@@ -418,7 +425,8 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
tkew=mtkeus(i,j)+mtkeni(i,j)+mtkebf(i,j)+mtkers(i,j)
if (.not.(nitr.eq.1.and.pres(3)*lbi.gt.1.)) then
dtke=(tkew-tkeo)/dpmxl
- if (abs(dtke)<(abs(tkew)+1.e-16)/(pres(3)-pres(1))) then
+ if (abs(dtke)<(abs(tkew)+1.e-22*V_mks2cgs)
+ . /(pres(3)-pres(1))) then
if (tkew.lt.0.) then
dpmxl=.5*(pres(1)-pmxl)
else
@@ -437,9 +445,9 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
write (lp,*) 'dpth=',pres(3)/onem,';'
write (lp,*) 'pmxl=',pmxl/onem,';'
write (lp,*) 'corio=',coriop(i,j),';'
- write (lp,*) 'ustar=',ustar(i,j)*1.e-2,';'
- write (lp,*) 'bfltot=',bfltot*1.e-4,';'
- write (lp,*) 'bflpsw=',bflpsw*1.e-4,';'
+ write (lp,*) 'ustar=',ustar(i,j)*iL_mks2cgs,';'
+ write (lp,*) 'bfltot=',bfltot*A_cgs2mks,';'
+ write (lp,*) 'bflpsw=',bflpsw*A_cgs2mks,';'
write (lp,*) 'bg2=',util1(i,j),';'
write (lp,*) 'ce=',ce*sqrt(scp2(i,j))*rlf,';'
write (lp,*)
@@ -512,7 +520,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c
kmax=1
do k=2,kk
- if (delp(k).gt.epsil) kmax=k
+ if (delp(k).gt.epsilp) kmax=k
enddo
kfmax=0
c
@@ -589,7 +597,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
bc(k)=0.
endif
kfmax=k
- if (bdpsum.le.epsil) then
+ if (bdpsum.le.epsilp) then
if (dpfsl.gt.onemu) then
bpmldp=min(bpmndp,dpfsl+delp(2))
q=brnflx(i,j)*delt1*g/bpmldp
@@ -730,7 +738,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
endif
else
if (delp(k).gt.onemu.and.dens(k).gt.densr(k).and.
- . sigfsl.lt.densr(k)-1.e-9) then
+ . sigfsl.lt.densr(k)-(1.e-6*R_mks2cgs)) then
dps=min(dpfsl,delp(k)*(dens(k)-densr(k))
. /(densr(k)-sigfsl))
q=1./(dps+delp(k))
@@ -857,7 +865,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
do
if (k.gt.kk) then
exit
- elseif (delp(k).lt.epsil) then
+ elseif (delp(k).lt.epsilp) then
k=k+1
else
pmxl=pres(k+1)
@@ -929,7 +937,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
endif
if (.not.chngd) then
if (abs(dtke).lt.
- . (abs(tkew)+1.e-16)/delp(k)) then
+ . (abs(tkew)+1.e-22*V_mks2cgs)/delp(k)) then
if (tkew.lt.0.) then
dpmxl=.5*(pres(k)-pmxl)
else
@@ -952,9 +960,9 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
write (lp,*) 'dpth=',pres(3)/onem,';'
write (lp,*) 'pmxl=',pmxl/onem,';'
write (lp,*) 'corio=',coriop(i,j),';'
- write (lp,*) 'ustar=',ustar(i,j)*1.e-2,';'
- write (lp,*) 'bfltot=',bfltot*1.e-4,';'
- write (lp,*) 'bflpsw=',bflpsw*1.e-4,';'
+ write (lp,*) 'ustar=',ustar(i,j)*iL_mks2cgs,';'
+ write (lp,*) 'bfltot=',bfltot*A_cgs2mks,';'
+ write (lp,*) 'bflpsw=',bflpsw*A_cgs2mks,';'
write (lp,*) 'bg2=',util1(i,j),';'
write (lp,*) 'ce=',ce*sqrt(scp2(i,j))*rlf,';'
write (lp,*) 'pres(3)=',pres(3)/onem,';'
@@ -972,7 +980,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c call xchalt('(mxlayr)')
c stop '(mxlayr)'
endif
- if (pmxl.lt.pres(k+1)-epsil.and.nitr.lt.maxitr) then
+ if (pmxl.lt.pres(k+1)-epsilp.and.nitr.lt.maxitr) then
tdps=tdps+ttem(k)*(pmxl-pres(k))
sdps=sdps+ssal(k)*(pmxl-pres(k))
#ifdef TRC
@@ -1072,7 +1080,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c
kmax=1
do k=2,kk
- if (delp(k).gt.epsil) kmax=k
+ if (delp(k).gt.epsilp) kmax=k
enddo
kfmax=0
c
@@ -1142,7 +1150,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
bc(k)=0.
endif
kfmax=k
- if (bdpsum.le.epsil) then
+ if (bdpsum.le.epsilp) then
ssal(2)=ssal(2)-brnflx(i,j)*delt1*g/delp(2)
else
if (bdpsum.lt.bpmndp) then
@@ -1206,7 +1214,7 @@ subroutine mxlayr(m,n,mm,nn,k1m,k1n)
c --- --- Define first physical layer.
k=3
dps=0.
- do while (delp(k).lt.epsil)
+ do while (delp(k).lt.epsilp)
dps=dps+delp(k)
delp(k)=0.
k=k+1
diff --git a/phy/mod_ndiff.F90 b/phy/mod_ndiff.F90
index 5f875313..d31b0f34 100644
--- a/phy/mod_ndiff.F90
+++ b/phy/mod_ndiff.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2022 Mats Bentsen
+! Copyright (C) 2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -23,7 +23,7 @@ module mod_ndiff
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: g, alpha0, epsil, onemm
+ use mod_constants, only: g, alpha0, epsilp, onemm, P_mks2cgs, R_mks2cgs
use mod_time, only: delt1
use mod_xc
use mod_grid, only: scuy, scvx, scp2, scuxi, scvyi
@@ -42,8 +42,8 @@ module mod_ndiff
private
real(r8), parameter :: &
- rhoeps = 1.e-8_r8, &
- dpeps = 1.e-4_r8
+ rhoeps = 1.e-5_r8*R_mks2cgs, &
+ dpeps = 1.e-5_r8*P_mks2cgs
integer, parameter :: &
p_ord = 4, &
it = 1, &
@@ -846,7 +846,7 @@ subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
flxconv_rs(kd_p,is,i_p,j_rs_p) - sflx
p_ni_up = .5_r8*(p_ni_m(nip) + p_ni_p(nip))
p_ni_lo = .5_r8*(p_ni_m(nic) + p_ni_p(nic))
- dp_ni_i = 1._r8/max(epsil, p_ni_lo - p_ni_up)
+ dp_ni_i = 1._r8/max(epsilp, p_ni_lo - p_ni_up)
do while (kuv <= kk)
kuvm = kuv + mm
if (puv(i_p,j_p,kuv+1) < p_ni_lo) then
@@ -918,7 +918,7 @@ subroutine ndiff_flx(p_srcdi_m, t_srcdi_m, tpc_src_m, &
ks = ks + 1
enddo
q = (p_nslp_src(ks) - p_nslp_dst) &
- /max(p_nslp_src(ks) - p_nslp_src(ks-1), epsil)
+ /max(p_nslp_src(ks) - p_nslp_src(ks-1), epsilp)
nslpxy(i_p,j_p,kd) = q*nslp_src(ks-1) + (1._r8 - q)*nslp_src(ks)
kd = kd + 1
if (kd > kk) exit
diff --git a/phy/mod_pbcor.F b/phy/mod_pbcor.F
index 3a350abb..b0d84df6 100644
--- a/phy/mod_pbcor.F
+++ b/phy/mod_pbcor.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2005-2020 Mats Bentsen, Mehmet Ilicak
+! Copyright (C) 2005-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -27,7 +27,7 @@ module mod_pbcor
c --- ------------------------------------------------------------------
c
use mod_types, only: r8
- use mod_constants, only: epsil
+ use mod_constants, only: epsilp, P_mks2cgs
use mod_time, only: dlt
use mod_xc
use mod_grid, only: scp2i
@@ -54,10 +54,10 @@ module mod_pbcor
c
c --- Parameters:
real(r8), parameter ::
- . dpeps1 = 1.e-4_r8, ! Small layer pressure thickness
- ! [g cm-1 s-2].
- . dpeps2 = 1.e-6_r8 ! Small layer pressure thickness
- ! [g cm-1 s-2].
+ . dpeps1 = 1.e-5_r8*P_mks2cgs, ! Small layer pressure thickness
+ ! [g cm-1 s-2].
+ . dpeps2 = 1.e-7_r8*P_mks2cgs ! Small layer pressure thickness
+ ! [g cm-1 s-2].
c
public :: bmcmth, pbcor1, pbcor2
c
@@ -459,7 +459,7 @@ subroutine pbcor2(m,n,mm,nn,k1m,k1n)
km=k+mm
do l=1,isp(j)
do i=max(0,ifp(j,l)),min(ii+1,ilp(j,l))
- dp(i,j,km)=max(0.,dp(i,j,km))+epsil
+ dp(i,j,km)=max(0.,dp(i,j,km))+epsilp
p(i,j,k+1)=p(i,j,k)+dp(i,j,km)
enddo
enddo
@@ -719,7 +719,7 @@ subroutine pbcor2(m,n,mm,nn,k1m,k1n)
enddo
#endif
sigma(i,j,km)=sig(temp(i,j,km),saln(i,j,km))
- dp(i,j,km)=dp(i,j,km)-epsil
+ dp(i,j,km)=dp(i,j,km)-epsilp
if (dp(i,j,km).lt.dpeps2) dp(i,j,km)=0.
enddo
enddo
diff --git a/phy/mod_pgforc.F b/phy/mod_pgforc.F
index 22eed8c6..7fdac606 100644
--- a/phy/mod_pgforc.F
+++ b/phy/mod_pgforc.F
@@ -25,7 +25,7 @@ module mod_pgforc
c --- ------------------------------------------------------------------
c
use mod_types, only: r8
- use mod_constants, only: g, epsil, spval
+ use mod_constants, only: g, epsilp, spval
use mod_xc
use mod_state, only: dp, dpu, dpv, temp, saln, p, pu, pv, phi,
. pb_p, pbu_p, pbv_p, sealv
@@ -140,9 +140,6 @@ subroutine pgforc(m,n,mm,nn,k1m,k1n)
c --- ------------------------------------------------------------------
c --- compute the pressure gradient force
c --- ------------------------------------------------------------------
-c
- use mod_constants, only: g, epsil
- use mod_xc
c
implicit none
c
@@ -206,7 +203,7 @@ subroutine pgforc(m,n,mm,nn,k1m,k1n)
kn=k+nn
do l=1,isp(j)
do i=max(0,ifp(j,l)),min(ii,ilp(j,l))
- if (dp(i,j,kn).lt.epsil) then
+ if (dp(i,j,kn).lt.epsilp) then
phi (i,j,k)=phi (i,j,k+1)
phip(i,j,k)=phip(i,j,k+1)
else
diff --git a/phy/mod_remap.F b/phy/mod_remap.F
index 44180dd1..10e19318 100644
--- a/phy/mod_remap.F
+++ b/phy/mod_remap.F
@@ -26,6 +26,7 @@ module mod_remap
c
use mod_types, only: r8
use mod_xc
+ use mod_constants, only: P_mks2cgs
#ifdef TRC
use mod_tracers, only: ntr, itrtke, itrgls
#endif
@@ -36,8 +37,8 @@ module mod_remap
c
c --- Parameters:
real(r8), parameter ::
- . dpeps = 1.e-11_r8 ! Small layer pressure thickness (equivalent
- ! to approximately 10-16 m) [g cm-1 s-2].
+ . dpeps = 1.e-12_r8*P_mks2cgs ! Small layer pressure thickness (equivalent
+ ! to approximately 10-16 m) [g cm-1 s-2].
#if defined(TRC) && defined(ATRC)
real(r8), parameter ::
. treps = 1.e-14_r8 ! Small tracer concentration.
diff --git a/phy/mod_swabs.F b/phy/mod_swabs.F
index 6b1b6b38..f8c57064 100644
--- a/phy/mod_swabs.F
+++ b/phy/mod_swabs.F
@@ -113,7 +113,7 @@ module mod_swabs
. ma94z2=(/ 7.925,-6.644, 3.662,-1.815, -.218, .502/)
c
c --- Other parameters:
-c---- swamxd: Maximum depth of shortwave radiation penetration.
+c---- swamxd: Maximum depth of shortwave radiation penetration [m].
real, parameter ::
. swamxd = 200.
c
diff --git a/phy/mod_tidaldissip.F90 b/phy/mod_tidaldissip.F90
index 2005289d..d3a61236 100644
--- a/phy/mod_tidaldissip.F90
+++ b/phy/mod_tidaldissip.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2020 Mats Bentsen
+! Copyright (C) 2015-2020 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -24,7 +24,7 @@ module mod_tidaldissip
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: spval
+ use mod_constants, only: spval, M_mks2cgs
use mod_xc
use mod_checksum, only: csdiag, chksummsk
use netcdf
@@ -157,7 +157,7 @@ subroutine read_tidaldissip
do j = 1, jj
do l = 1, isp(j)
do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
- twedon(i, j) = twedon(i, j)*1.e3_r8
+ twedon(i, j) = twedon(i, j)*M_mks2cgs
enddo
enddo
enddo
diff --git a/phy/mod_time.F90 b/phy/mod_time.F90
index 5e32c7d3..f4f0442b 100644
--- a/phy/mod_time.F90
+++ b/phy/mod_time.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2020-2021 Mats Bentsen, Mehmet Ilicak, Aleksi Nummelin
+! Copyright (C) 2020-2022 Mats Bentsen, Mehmet Ilicak, Aleksi Nummelin
!
! This file is part of BLOM.
!
@@ -24,7 +24,7 @@ module mod_time
use mod_types, only: r8
use mod_config, only: expcnf
- use mod_constants, only: epsil
+ use mod_constants, only: epsilt
use mod_calendar, only: date_type, daynum_diff, date_offset, &
calendar_noerr, calendar_errstr
use mod_xc, only: lp, mnproc, xcstop
@@ -118,7 +118,7 @@ subroutine init_timevars
! Get number of baroclinic time steps per day and verify that an integer
! number of steps fits in a day.
nstep_in_day = nint(86400._r8/baclin)
- if (abs(86400._r8/baclin - nstep_in_day) > epsil) then
+ if (abs(86400._r8/baclin - nstep_in_day) > epsilt) then
if (mnproc == 1) then
write (lp, *) &
'init_timevars: '// &
diff --git a/phy/mod_tke.F90 b/phy/mod_tke.F90
index fc5655ae..73dc2c95 100644
--- a/phy/mod_tke.F90
+++ b/phy/mod_tke.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2013-2020 Mehmet Ilicak, Mats Bentsen
+! Copyright (C) 2013-2022 Mehmet Ilicak, Mats Bentsen
!
! This file is part of BLOM.
!
@@ -24,6 +24,7 @@ module mod_tke
! ------------------------------------------------------------------------------
use mod_types, only: r8
+ use mod_constants, only: spval
use mod_xc
use mod_diffusion, only: difdia
use mod_forcing, only: ustarb
@@ -34,10 +35,8 @@ module mod_tke
real(r8), parameter :: &
gls_cmu0 = .527_r8, & ! cmu0
- Pr_t = 1._r8, & ! Turbulent Prandtl number [].
- tke_min = 7.6e-4_r8, & ! Minimum TKE value [?].
+ Pr_t = 1._r8, & ! Turbulent Prandtl number [non-dimensional].
zos = .0002_r8, & !
- gls_psi_min = 1.e-10_r8, & ! Minimum GLS value [?].
gls_p = 3._r8, & !
gls_m = 1.5_r8, & !
gls_n = -1._r8, & !
@@ -56,8 +55,19 @@ module mod_tke
gls_Gh0 = .0329_r8, & !
gls_Ghmin = -.28_r8, & !
gls_Ghcri = .03_r8, & !
- vonKar = .4_r8, & !
- Ls_unlmt_min = 1.e-6_r8 !
+ vonKar = .4_r8 !
+
+#ifdef MKS
+ real(r8), parameter :: &
+ tke_min = 7.6e-8_r8, & ! Minimum TKE value [m2/s2].
+ gls_psi_min = 1.e-14_r8, & ! Minimum GLS value [m2/s3].
+ Ls_unlmt_min = 1.e-8_r8 ! [m]
+#else
+ real(r8), parameter :: &
+ tke_min = 7.6e-4_r8, & ! Minimum TKE value [cm2/s2].
+ gls_psi_min = 1.e-10_r8, & ! Minimum GLS value [cm2/s3].
+ Ls_unlmt_min = 1.e-6_r8 ! [cm]
+#endif
real(r8), dimension(1 - nbdy:idm + nbdy, 1 - nbdy:jdm + nbdy, kdm) :: &
Prod, & ! Shear production [?].
@@ -93,6 +103,18 @@ subroutine initke
! Initialize fields holding turbulent kinetic energy, generic length
! scale, and other fields used in the turbulence closure.
+ !$omp parallel do private(i, k)
+ do j = 1 - nbdy, jj + nbdy
+ do i = 1 - nbdy, ii + nbdy
+ do k = 1, kk
+ Prod(i ,j ,k) = spval
+ Buoy(i ,j ,k) = spval
+ Shear2(i ,j ,k) = spval
+ L_scale(i ,j ,k) = spval
+ enddo
+ enddo
+ enddo
+ !$omp end parallel do
!$omp parallel do private(k, l, i)
do j = 1 - nbdy, jj + nbdy
do k = 1, 2*kdm
diff --git a/phy/mod_tmsmt.F b/phy/mod_tmsmt.F
index 5aec8d56..cf7ef21a 100644
--- a/phy/mod_tmsmt.F
+++ b/phy/mod_tmsmt.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2005-2022 Mats Bentsen
+! Copyright (C) 2005-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -27,7 +27,7 @@ module mod_tmsmt
c --- ------------------------------------------------------------------
c
use mod_types, only: r8
- use mod_constants, only: epsil, spval
+ use mod_constants, only: epsilp, spval
use mod_xc
use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml
use mod_state, only: dp, dpu, dpv, temp, saln, p, pb
@@ -265,7 +265,7 @@ subroutine tmsmt2(m,n,mm,nn,k1m,k1n)
c
c --- time smoothing of layer thickness, temperature and salinity
c
- use mod_constants, only: epsil
+ use mod_constants, only: epsilp
use mod_xc
c
implicit none
@@ -318,23 +318,23 @@ subroutine tmsmt2(m,n,mm,nn,k1m,k1n)
pnew=max(0.,dp(i,j,kn)*pbfacn(i))
dp(i,j,km)=wts1*pmid+wts2*(pold+pnew)
dpold(i,j,km)=dp(i,j,km)
- pold=pold+epsil
- pmid=pmid+epsil
- pnew=pnew+epsil
+ pold=pold+epsilp
+ pmid=pmid+epsilp
+ pnew=pnew+epsilp
temp(i,j,km)=(wts1*pmid*temp(i,j,km)
. +wts2*(pold*told(i,j,k)+pnew*temp(i,j,kn)))
- . /(dp(i,j,km)+epsil)
+ . /(dp(i,j,km)+epsilp)
told(i,j,k)=temp(i,j,km)
saln(i,j,km)=(wts1*pmid*saln(i,j,km)
. +wts2*(pold*sold(i,j,k)+pnew*saln(i,j,kn)))
- . /(dp(i,j,km)+epsil)
+ . /(dp(i,j,km)+epsilp)
sold(i,j,k)=saln(i,j,km)
#ifdef TRC
do nt=1,ntr
trc(i,j,km,nt)=(wts1*pmid*trc(i,j,km,nt)
. +wts2*(pold*trcold(i,j,k,nt)
. +pnew*trc(i,j,kn,nt)))
- . /(dp(i,j,km)+epsil)
+ . /(dp(i,j,km)+epsilp)
trcold(i,j,k,nt)=trc(i,j,km,nt)
enddo
#endif
diff --git a/phy/mod_vcoord.F90 b/phy/mod_vcoord.F90
index 849380b9..f63e6cbd 100644
--- a/phy/mod_vcoord.F90
+++ b/phy/mod_vcoord.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2021-2022 Mats Bentsen
+! Copyright (C) 2021-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -25,7 +25,7 @@ module mod_vcoord
use mod_types, only: r8
use mod_config, only: inst_suffix
- use mod_constants, only: g, epsil, spval, onem
+ use mod_constants, only: g, epsilp, spval, onem
use mod_xc
use mod_eos, only: sig, dsigdt, dsigds
use mod_state, only: u, v, dp, dpu, dpv, temp, saln, sigma, p, pu, pv
@@ -296,7 +296,7 @@ subroutine cntiso_regrid_direct_jslice(p_src, p_dst, i_lb, i_ub, j, j_rs, nn)
kl = kk
ku = kl - 1
do while (ku > 0)
- thin_layers = p_src(kl+1,i) - p_src(ku,i) < epsil
+ thin_layers = p_src(kl+1,i) - p_src(ku,i) < epsilp
if (thin_layers .or. &
sigma_1d(kl) - sigma_1d(ku) &
< .5_r8*beta*(p_src(kl+1,i) - p_src(ku,i))) then
@@ -311,7 +311,7 @@ subroutine cntiso_regrid_direct_jslice(p_src, p_dst, i_lb, i_ub, j, j_rs, nn)
ku = ku - 1
sdpsum = sdpsum &
+ sigma_1d(ku)*(p_src(ku+1,i) - p_src(ku,i))
- thin_layers = p_src(kl+1,i) - p_src(ku,i) < epsil
+ thin_layers = p_src(kl+1,i) - p_src(ku,i) < epsilp
if (.not. thin_layers) &
smean = sdpsum/(p_src(kl+1,i) - p_src(ku,i))
layer_added = .true.
@@ -331,7 +331,7 @@ subroutine cntiso_regrid_direct_jslice(p_src, p_dst, i_lb, i_ub, j, j_rs, nn)
kl = kl + 1
sdpsum = sdpsum &
+ sigma_1d(kl)*(p_src(kl+1,i) - p_src(kl,i))
- thin_layers = p_src(kl+1,i) - p_src(ku,i) < epsil
+ thin_layers = p_src(kl+1,i) - p_src(ku,i) < epsilp
if (.not. thin_layers) &
smean = sdpsum/(p_src(kl+1,i) - p_src(ku,i))
layer_added = .true.
diff --git a/phy/mod_vdiff.F90 b/phy/mod_vdiff.F90
index 3d1798fd..61778a72 100644
--- a/phy/mod_vdiff.F90
+++ b/phy/mod_vdiff.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2021-2022 Mats Bentsen
+! Copyright (C) 2021-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -23,7 +23,7 @@ module mod_vdiff
! ------------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: g, spcifh, alpha0
+ use mod_constants, only: g, spcifh, alpha0, onem
use mod_time, only: delt1
use mod_xc
use mod_eos, only: sig
@@ -40,7 +40,7 @@ module mod_vdiff
private
real(r8), parameter :: &
- dpmin_vdiff = 0.1_r8*98060._r8
+ dpmin_vdiff = 0.1_r8*onem
public :: cntiso_hybrid_vdifft, cntiso_hybrid_vdiffm
diff --git a/phy/numerical_bounds.F90 b/phy/numerical_bounds.F90
index f4579056..a38c5b55 100644
--- a/phy/numerical_bounds.F90
+++ b/phy/numerical_bounds.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2020 Mats Bentsen
+! Copyright (C) 2020-2022 Mats Bentsen, Mehmet Ilicak
!
! This file is part of BLOM.
!
@@ -23,7 +23,7 @@ subroutine numerical_bounds
! ---------------------------------------------------------------------------
use mod_types, only: r8
- use mod_constants, only: g, spval
+ use mod_constants, only: g, spval, L_mks2cgs
use mod_time, only: baclin
use mod_xc
use mod_grid, only: scqx, scqy, scpx, scpy, scuy, scvx, scp2, depths
@@ -61,8 +61,8 @@ subroutine numerical_bounds
do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
btdtmx = min(btdtmx, &
scpx(i, j)*scpy(i, j) &
- /sqrt(g*depths(i, j)*100._r8*( scpx(i, j)*scpx(i, j) &
- + scpy(i, j)*scpy(i, j))))
+ /sqrt(g*depths(i, j)*L_mks2cgs*( scpx(i, j)*scpx(i, j) &
+ + scpy(i, j)*scpy(i, j))))
enddo
enddo
enddo
diff --git a/phy/rdlim.F b/phy/rdlim.F
index 44949352..f4d00f76 100644
--- a/phy/rdlim.F
+++ b/phy/rdlim.F
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2008-2021 Mats Bentsen, Mehmet Ilicak, Ingo Bethke,
+! Copyright (C) 2008-2022 Mats Bentsen, Mehmet Ilicak, Ingo Bethke,
! Ping-Gin Chiu, Aleksi Nummelin
!
! This file is part of BLOM.
@@ -25,6 +25,7 @@ subroutine rdlim
c --- ------------------------------------------------------------------
c
use mod_config, only: expcnf, runid, inst_suffix
+ use mod_constants, only: epsilt
use mod_calendar, only: date_type, daynum_diff, calendar_errstr,
. calendar_noerr, operator(==), operator(<),
. operator(/=)
@@ -732,7 +733,7 @@ subroutine rdlim
c
c --- - verify integer number of baroclinic time steps per coupling
c --- - interval
- if (mod(ocn_cpl_dt_cesm+epsil,baclin).gt.2.*epsil) then
+ if (mod(ocn_cpl_dt_cesm+epsilt,baclin).gt.2.*epsilt) then
if (mnproc.eq.1) then
write (lp,*) 'rdlim: must have an integer number of '//
. 'baroclinic time steps in a coupling'
diff --git a/single_column/mod_single_column.F90 b/single_column/mod_single_column.F90
index 58eb6826..6017c0d9 100644
--- a/single_column/mod_single_column.F90
+++ b/single_column/mod_single_column.F90
@@ -1,5 +1,5 @@
! ------------------------------------------------------------------------------
-! Copyright (C) 2021 Mehmet Ilicak, Mats Bentsen
+! Copyright (C) 2021-2022 Mehmet Ilicak, Mats Bentsen
!
! This file is part of BLOM.
!
@@ -24,6 +24,7 @@ module mod_single_column
! ----------------------------------------------------------------------
use mod_types, only: r8
+ use mod_constants, only: L_mks2cgs
use mod_xc
use mod_vcoord, only: sigmar
use mod_grid, only: qclon, qclat, pclon, pclat, uclon, uclat, vclon, vclat, &
@@ -64,13 +65,13 @@ subroutine geoenv_single_column
uclat = 0._r8
vclon = 0._r8
vclat = 0._r8
- scqx = 1100000.0_r8
- scqy = 1100000.0_r8
- scpx = 1100000.0_r8
- scpy = 1100000.0_r8
- scux = 1100000.0_r8
- scuy = 1100000.0_r8
- scvx = 1100000.0_r8
+ scqx = 11000.0_r8*L_mks2cgs
+ scqy = 11000.0_r8*L_mks2cgs
+ scpx = 11000.0_r8*L_mks2cgs
+ scpy = 11000.0_r8*L_mks2cgs
+ scux = 11000.0_r8*L_mks2cgs
+ scuy = 11000.0_r8*L_mks2cgs
+ scvx = 11000.0_r8*L_mks2cgs
scvy = scuy
scq2 = scqx*scqy
scp2 = scpx*scpy