Merge pull request #217 from matsbn/SI_units

Updates and corrections for CGS/MKS unit systems

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.
 !
@@ -183,14 +183,13 @@ module mod_ben02
      .  atm_cswa_era       ! short-wave radiation adjustment factor
                            ! (NCEP)
 c
-#if defined(CGS)
-      data atm_ice_csmt_ncep,atm_rnf_csmt_ncep /2.e14,1.e13/,
+#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.e13/,
+     .     atm_ice_csmt_era,atm_rnf_csmt_era /0.0,1.e9/,
      .     atm_crnf_era,atm_cswa_era /0.7234,0.9721/
-#endif
-#if defined(MKS)
-      data atm_ice_csmt_ncep,atm_rnf_csmt_ncep /2.e10,1.e9/,
+#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/
@@ -2098,13 +2097,12 @@ 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
-      real iL_mks2cgssq
-      iL_mks2cgssq = 1.0 / (L_mks2cgs*L_mks2cgs)
+      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),
@@ -2786,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)*iL_mks2cgssq ! [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
@@ -2830,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)*iL_mks2cgssq ! [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))

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.
 !

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,8 +21,8 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
 c
 c --- NERSC version of thermf. 
 c
-      use mod_constants, only: spcifh, t0deg, epsilt, onem
-      use mod_constants, only: L_mks2cgs, M_mks2cgs, g2kg, alpha0
+      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
@@ -67,12 +67,12 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
       real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: vrtsfl
 c
       integer i,j,k,l,m1,m2,m3,m4,m5
-      real*8 dt,cpsw,rnf_fac,sag_fac,y,
+      real dt,cpsw,rnf_fac,sag_fac,y,
      .     dpotl,hotl,totl,sotl,dpmxl,hmxl,tmxl,smxl,tice_f,hice_min,
      .     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) ::
@@ -83,9 +83,8 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
 c
       real intp1d
       external intp1d
-      real iL_mks2cgssq
 c
-      iL_mks2cgssq = 1.0 / (L_mks2cgs**2)
+      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
@@ -419,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))*iL_mks2cgssq
+          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)*iL_mks2cgssq
+          sswflx(i,j)=-qsww*(1.-fice0)*A_cgs2mks
 c
 #ifdef TRC
 c --- ------------------------------------------------------------------
@@ -478,7 +477,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
      .                  ricclm(i,j,l5mi),xmi)
             sstc=(1.-rice)*max(sstc,tice_f)+rice*tice_f
             trxflx=spcifh*L_mks2cgs*min(hmxl,trxdpt)/(trxday*86400.)
-     .                        *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
+     .             *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
             surrlx(i,j)=-trxflx
           else
             trxflx=0.
@@ -505,7 +504,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
      .                  sssclm(i,j,l3mi),sssclm(i,j,l4mi),
      .                  sssclm(i,j,l5mi),xmi)
             srxflx=L_mks2cgs*min(hmxl,srxdpt)/(srxday*86400.)
-     .                 *min(srxlim,max(-srxlim,sssc-smxl))/alpha0
+     .             *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)
@@ -572,9 +571,9 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
         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))
-     .                 *(1e3*(M_mks2cgs/L_mks2cgs**2))
+     .                 *(kg2g*(M_mks2cgs/L_mks2cgs**2))
           brnflx(i,j)=-brnflx(i,j)
-     .                 *(1e3*(M_mks2cgs/L_mks2cgs**2))
+     .                 *(kg2g*(M_mks2cgs/L_mks2cgs**2))
         enddo
         enddo
       enddo

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.
 !

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,8 +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, epsilt, onem
-      use mod_constants, only: L_mks2cgs, M_mks2cgs, g2kg, alpha0
+      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
@@ -62,19 +62,18 @@ 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) ::
      .  ttrsf,ttrav
       real tottrsf,tottrav,trflxc
 #endif
 c
-      real iL_mks2cgssq
       real intp1d
       external intp1d
 c
-      iL_mks2cgssq = 1.0 / (L_mks2cgs*L_mks2cgs)
+      A_cgs2mks=1./(L_mks2cgs**2)
 c      
 c --- Set parameters for time interpolation when applying diagnosed heat
 c --- and salt relaxation fluxes
@@ -155,7 +154,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
 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)*(L_mks2cgs**2)
+     .                /(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*(L_mks2cgs**2)
@@ -164,11 +163,11 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
           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))*iL_mks2cgssq
+          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)*iL_mks2cgssq
+          sswflx(i,j)=-swa(i,j)*A_cgs2mks
 c
 #ifdef TRC
 c --- ------------------------------------------------------------------
@@ -187,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
@@ -206,6 +205,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
             ttrav(nt,i,j)=trc(i,j,k1n,nt)*scp2(i,j)
           enddo
 #endif
+c
 c --- ------------------------------------------------------------------
 c --- --- Relaxation fluxes
 c --- ------------------------------------------------------------------
@@ -222,7 +222,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
      .                  ricclm(i,j,l5mi),xmi)
             sstc=(1.-rice)*max(sstc,tice_f)+rice*tice_f
             trxflx=spcifh*L_mks2cgs*min(hmxl,trxdpt)/(trxday*86400.)
-     .                        *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
+     .             *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
             surrlx(i,j)=-trxflx
           else
             trxflx=0.
@@ -249,7 +249,7 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
      .                  sssclm(i,j,l3mi),sssclm(i,j,l4mi),
      .                  sssclm(i,j,l5mi),xmi)
             srxflx=L_mks2cgs*min(hmxl,srxdpt)/(srxday*86400.)
-     .                 *min(srxlim,max(-srxlim,sssc-smxl))/alpha0
+     .             *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)
@@ -304,9 +304,9 @@ subroutine thermf_cesm(m,n,mm,nn,k1m,k1n)
         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))
-     .                 *(1e3*(M_mks2cgs/L_mks2cgs**2))
+     .                 *(kg2g*(M_mks2cgs/L_mks2cgs**2))
           brnflx(i,j)=-brnflx(i,j)
-     .                 *(1e3*(M_mks2cgs/L_mks2cgs**2))
+     .                 *(kg2g*(M_mks2cgs/L_mks2cgs**2))
         enddo
         enddo
       enddo

cime_config/buildcpp

@@ -87,7 +87,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")
-    cgsmks = case.get_value("BLOM_UNIT")
+    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))
 
@@ -147,14 +147,10 @@ def buildcpp(case):
             else:
                 expect(False, "tracer module {} is not recognized".format(module))
 
-    if cgsmks:
-        for option in cgsmks.split():
-            if option == "cgs":
-                blom_cppdefs = blom_cppdefs + " -DCGS"
-            elif option == "mks":
-                blom_cppdefs = blom_cppdefs + " -DMKS"
-            else:
-                expect(False, "SI_UNIT module {} is not recognized".format(option))
+    if blom_unit == "mks":
+        blom_cppdefs = blom_cppdefs + " -DMKS"
+    else:
+        expect(blom_unit == "cgs", "Unit system {} is not recognized".format(option))
 
     blom_cppdefs = "-DMPI" + blom_cppdefs