Merge pull request #255 from matsbn/feature-hybrid_enhancements2

Hybrid coordinate enhancements 2

bld/tnx0.125v4/kdm.cntiso_hybrid

@@ -0,0 +1 @@
+56

cime_config/buildnml

@@ -136,6 +136,7 @@ set DITFLX   = .false.
 set DISFLX   = .false.
 set SRXBAL   = .false.
 set SCFILE   = "'unset'"
+set WAVSRC   = "'none'"
 set SMTFRC   = .true.
 set SPRFAC   = .false.
 set ATM_PATH = "'unset'"
@@ -207,6 +208,8 @@ else
   set BDMC2    = 1.e-5
 endif
 set TKEPF    = .006
+set SMOBLD   = .true.
+set LNGMTP   = "'none'"
 if ($BLOM_VCOORD == isopyc_bulkml) then
   set BDMLDP   = .true.
   set LTEDTP   = "'layer'"
@@ -231,7 +234,7 @@ if ($BLOM_N_DEPOSITION == TRUE) then
   else if( $BLOM_NDEP_SCENARIO == 1850 && $OCN_GRID == tnx0.25v4) then
     set NDEPFNAME = ndep_1850_CMIP6_tnx0.25v4_20190912.nc
   else if( $BLOM_NDEP_SCENARIO == 1850 && $OCN_GRID == tnx0.125v4) then
-    set NDEPFNAME = ndep_1850_CMIP6_tnx0.125v4_20190912.nc
+    set NDEPFNAME = ndep_1850_CMIP6_tnx0.125v4_20221013.nc
   else if( $BLOM_NDEP_SCENARIO == 2000 && $OCN_GRID == tnx2v1) then
     set NDEPFNAME = ndep_2000_CMIP6_tnx2v1_20200826.nc
   else if( $BLOM_NDEP_SCENARIO == 2000 && $OCN_GRID == tnx1v4) then
@@ -311,6 +314,8 @@ set H2D_HMLTFZ   =  '0,   4,   0'
 set H2D_HSNW     =  '0,   0,   0'
 set H2D_IAGE     =  '0,   0,   0'
 set H2D_IDKEDT   =  '0,   4,   0'
+set H2D_LAMULT   =  '0,   4,   0'
+set H2D_LASL     =  '0,   4,   0'
 set H2D_LIP      =  '0,   4,   0'
 set H2D_MAXMLD   =  '4,   4,   0'
 set H2D_MLD      =  '0,   4,   0'
@@ -355,8 +360,10 @@ set H2D_UB       =  '0,   4,   0'
 set H2D_UICE     =  '0,   0,   0'
 set H2D_USTAR    =  '0,   4,   0'
 set H2D_USTAR3   =  '0,   4,   0'
+set H2D_USTOKES  =  '0,   0,   0'
 set H2D_VB       =  '0,   4,   0'
 set H2D_VICE     =  '0,   0,   0'
+set H2D_VSTOKES  =  '0,   0,   0'
 set H2D_ZTX      =  '0,   4,   0'
 set LYR_BFSQ     =  '0,   4,   0'
 set LYR_DIFDIA   =  '0,   4,   0'
@@ -776,24 +783,25 @@ else if ($OCN_GRID == tnx0.125v4) then
   set CWBDTS = .75e-4
   set CWBDLS = 25
   if ($BLOM_UNIT == cgs) then
-    set MDV2HI = .5
+    set MDV2HI = .1
     set MDV2LO = .1
     set MDV4HI = 0.
     set MDV4LO = 0.
     set MDC2HI = 300.e4
     set MDC2LO = 100.e4
   else
-    set MDV2HI = .005
+    set MDV2HI = .001
     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
-  set VSC4LO = 0.0
+  set VSC2HI = 0.
+  set VSC2LO = 0.
+  set VSC4HI = .06
+  set VSC4LO = .06
+  set LTEDTP = "'layer'"
 else
   echo "$0 ERROR: Cannot deal with GRID = $OCN_GRID "
   exit -1
@@ -885,19 +893,19 @@ else if ($OCN_GRID == tnx0.25v4) then
     set NDEPFILE  = "''"
   endif
 else if ($OCN_GRID == tnx0.125v4) then
-  set GRFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/grid_tnx0.125v4_20200722.nc'"
-  set ICFILE = "'$DIN_LOC_ROOT/ocn/blom/inicon/inicon_tnx0.125v4_20200728.nc'"
-  set TDFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/tidal_dissipation_tnx0.125v4_20200722.nc'"
-  set MER_ORFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/ocean_regions_tnx0.125v4_20200722.nc'"
-  set MER_MIFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/mertra_index_tnx0.125v4_20200722.dat'"
-  set SEC_SIFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/section_index_tnx0.125v4_20200722.dat'"
-  set CCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/chlorophyll_concentration_tnx0.125v4_20200722.nc'"
-  set SCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/sss_clim_core_tnx0.125v4_20200722.nc'"
-  set FEDEPFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/dustdep_mhw2006_tnx0.125v4_20200722.nc'"
+  set GRFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/grid_tnx0.125v4_20221013.nc'"
+  set ICFILE = "'$DIN_LOC_ROOT/ocn/blom/inicon/inicon_tnx0.125v4_20230318.nc'"
+  set TDFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/tidal_dissipation_tnx0.125v4_20221013.nc'"
+  set MER_ORFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/ocean_regions_tnx0.125v4_20221013.nc'"
+  set MER_MIFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/mertra_index_tnx0.125v4_20221013.dat'"
+  set SEC_SIFILE = "'$DIN_LOC_ROOT/ocn/blom/grid/section_index_tnx0.125v4_20221013.dat'"
+  set CCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/chlorophyll_concentration_tnx0.125v4_20221013.nc'"
+  set SCFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/sss_clim_core_tnx0.125v4_20221013.nc'"
+  set FEDEPFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/dustdep_mhw2006_tnx0.125v4_20221013.nc'"
   set SWACLIMFILE = "''"
   set SEDPORFILE = "''"
   if ($BLOM_RIVER_NUTRIENTS == TRUE) then
-    set RIVINFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/river_nutrients_GNEWS2000_tnx0.125v4_20170821.nc'"
+    set RIVINFILE = "'$DIN_LOC_ROOT/ocn/blom/bndcon/river_nutrients_GNEWS2000_tnx0.125v4_20221013.nc'"
   else
     set RIVINFILE = "''"
   endif
@@ -1007,6 +1015,7 @@ cat >! $RUNDIR/ocn_in$inststr << EOF
 ! DISFLX   : Diagnose freshwater flux flag (l)
 ! SRXBAL   : Balance the SSS relaxation (l)
 ! SCFILE   : Name of file containing SSS climatology used for relaxation (a)
+! WAVSRC   : Source of wave fields. Valid source: 'none', 'param', 'extern' (a)
 ! SMTFRC   : Smooth CESM forcing (l)
 ! SPRFAC   : Send precipitation/runoff factor to CESM coupler (l)
 ! ATM_PATH : Path to forcing fields in case of EXPCNF 'ben02clim' or
@@ -1074,6 +1083,7 @@ cat >! $RUNDIR/ocn_in$inststr << EOF
   DISFLX   = $DISFLX
   SRXBAL   = $SRXBAL
   SCFILE   = $SCFILE
+  WAVSRC   = $WAVSRC
   SMTFRC   = $SMTFRC
   SPRFAC   = $SPRFAC
   ATM_PATH = $ATM_PATH
@@ -1157,6 +1167,10 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
 !            Gregg et al. (2003) (l)
 ! TKEPF    : Fraction of surface TKE that penetrates beneath mixed layer
 !            () (f)
+! SMOBLD   : If true, apply lateral smoothing of CVMix estimated
+!            boundary layer depth (l)
+! LNGMTP   : Type of CVMix Langmuir turbulence parameterization. Valid
+!            types: 'none', 'vr12-ma', 'lf17'
 ! LTEDTP   : Type of lateral tracer eddy diffusion: Valid methods:
 !            'layer', 'neutral'.
 &DIFFUSION
@@ -1182,6 +1196,8 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   BDMC2    = $BDMC2
   BDMLDP   = $BDMLDP
   TKEPF    = $TKEPF
+  SMOBLD   = $SMOBLD
+  LNGMTP   = $LNGMTP
   LTEDTP   = $LTEDTP
 /
 EOF
@@ -1309,6 +1325,8 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
 !   HSNW     - snow depth [m]
 !   IAGE     - ice age [d]
 !   IDKEDT   - mixed layer inertial kinetic energy tendency [kg s-3]
+!   LAMULT   - Langmuir enhancement factor []
+!   LASL     - Surface layer averaged Langmuir number []
 !   LIP      - liquid precipitation [kg m-2 s-1]
 !   MAXMLD   - maximum mixed layer depth [m]
 !   MLD      - mixed layer depth [m]
@@ -1353,8 +1371,10 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
 !   UICE     - ice velocity x-component [m s-1]
 !   USTAR    - friction velocity [m s-1]
 !   USTAR3   - friction velocity cubed [m3 s-3]
+!   USTOKES  - surface Stokes drift x-componen [m s-1]
 !   VB       - barotropic velocity y-component [m s-1]
 !   VICE     - ice velocity y-component [m s-1]
+!   VSTOKES  - surface Stokes drift y-componen [m s-1]
 !   ZTX      - wind stress x-component [N m-2]
 !   BFSQ     - buoyancy frequency squared [s-1]
 !   DIFDIA   - vertical diffusivity [log10(m2 s-1)|m2 s-1]
@@ -1440,6 +1460,8 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   H2D_HSNW     = $H2D_HSNW
   H2D_IAGE     = $H2D_IAGE
   H2D_IDKEDT   = $H2D_IDKEDT
+  H2D_LAMULT   = $H2D_LAMULT
+  H2D_LASL     = $H2D_LASL
   H2D_LIP      = $H2D_LIP
   H2D_MAXMLD   = $H2D_MAXMLD
   H2D_MLD      = $H2D_MLD
@@ -1484,8 +1506,10 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   H2D_UICE     = $H2D_UICE
   H2D_USTAR    = $H2D_USTAR
   H2D_USTAR3   = $H2D_USTAR3
+  H2D_USTOKES  = $H2D_USTOKES
   H2D_VB       = $H2D_VB
   H2D_VICE     = $H2D_VICE
+  H2D_VSTOKES  = $H2D_VSTOKES
   H2D_ZTX      = $H2D_ZTX
   LYR_BFSQ     = $LYR_BFSQ
   LYR_DIFDIA   = $LYR_DIFDIA

drivers/nuopc/mod_nuopc_methods.F90

@@ -30,7 +30,8 @@ module mod_nuopc_methods
    use mod_grid, only: scuy, scvx, scp2, scuxi, scvyi, plon, plat, &
                        cosang, sinang
    use mod_state, only: u, v, dp, temp, saln, pbu, pbv, ubflxs, vbflxs, sealv
-   use mod_forcing, only: sprfac, prfac, flxco2, flxdms, flxbrf
+   use mod_forcing, only: wavsrc_opt, wavsrc_extern, sprfac, prfac, flxco2, &
+                          flxdms, flxbrf
    use mod_difest, only: obldepth
    use mod_vcoord, only: vcoord_type_tag, isopyc_bulkml, cntiso_hybrid
    use mod_cesm, only: frzpot, mltpot, &
@@ -687,67 +688,71 @@ subroutine blom_importflds(fldlist_num, fldlist)
       call fill_global(mval, fval, halo_ps, abswnd_da(1-nbdy,1-nbdy,l2ci))
       call fill_global(mval, fval, halo_ps, ficem_da(1-nbdy,1-nbdy,l2ci))
 
-      call getfldindex(fldlist_num, fldlist, 'Sw_lamult', index_Sw_lamult)
-      call getfldindex(fldlist_num, fldlist, 'Sw_ustokes', index_Sw_ustokes)
-      call getfldindex(fldlist_num, fldlist, 'Sw_vstokes', index_Sw_vstokes)
-      call getfldindex(fldlist_num, fldlist, 'Sw_hstokes', index_Sw_hstokes)
+      if (wavsrc_opt == wavsrc_extern) then
 
-   !$omp parallel do private(i, n, utmp, vtmp)
-      do j = 1, jjcpl
-         do i = 1, ii
-            if     (ip(i,j) == 0) then
-               util1(i,j) = mval
-               util2(i,j) = mval
-               lamult_da(i,j,l2ci) = mval
-               lasl_da(i,j,l2ci) = mval
-            elseif (cplmsk(i,j) == 0) then
-               util1(i,j) = fval
-               util2(i,j) = fval
-               lamult_da(i,j,l2ci) = fval
-               lasl_da(i,j,l2ci) = fval
-            else
-               n = (j - 1)*ii + i
+         call getfldindex(fldlist_num, fldlist, 'Sw_lamult', index_Sw_lamult)
+         call getfldindex(fldlist_num, fldlist, 'Sw_ustokes', index_Sw_ustokes)
+         call getfldindex(fldlist_num, fldlist, 'Sw_vstokes', index_Sw_vstokes)
+         call getfldindex(fldlist_num, fldlist, 'Sw_hstokes', index_Sw_hstokes)
 
-               utmp = fldlist(index_Sw_ustokes)%dataptr(n)
-               vtmp = fldlist(index_Sw_vstokes)%dataptr(n)
-               util1(i,j) =   utmp*cosang(i,j) + vtmp*sinang(i,j)
-               util2(i,j) = - utmp*sinang(i,j) + vtmp*cosang(i,j)
+      !$omp parallel do private(i, n, utmp, vtmp)
+         do j = 1, jjcpl
+            do i = 1, ii
+               if     (ip(i,j) == 0) then
+                  util1(i,j) = mval
+                  util2(i,j) = mval
+                  lamult_da(i,j,l2ci) = mval
+                  lasl_da(i,j,l2ci) = mval
+               elseif (cplmsk(i,j) == 0) then
+                  util1(i,j) = fval
+                  util2(i,j) = fval
+                  lamult_da(i,j,l2ci) = fval
+                  lasl_da(i,j,l2ci) = fval
+               else
+                  n = (j - 1)*ii + i
 
-               ! Langmuir enhancement factor [].
-               lamult_da(i,j,l2ci) = fldlist(index_Sw_lamult)%dataptr(n)
+                  utmp = fldlist(index_Sw_ustokes)%dataptr(n)
+                  vtmp = fldlist(index_Sw_vstokes)%dataptr(n)
+                  util1(i,j) =   utmp*cosang(i,j) + vtmp*sinang(i,j)
+                  util2(i,j) = - utmp*sinang(i,j) + vtmp*cosang(i,j)
 
-               ! Surface layer averaged Langmuir number [].
-               lasl_da(i,j,l2ci) = fldlist(index_Sw_hstokes)%dataptr(n)
+                  ! Langmuir enhancement factor [].
+                  lamult_da(i,j,l2ci) = fldlist(index_Sw_lamult)%dataptr(n)
 
-            endif
+                  ! Surface layer averaged Langmuir number [].
+                  lasl_da(i,j,l2ci) = fldlist(index_Sw_hstokes)%dataptr(n)
+
+               endif
+            enddo
          enddo
-      enddo
-   !$omp end parallel do
+      !$omp end parallel do
 
-      call fill_global(mval, fval, halo_pv, util1)
-      call fill_global(mval, fval, halo_pv, util2)
-      call fill_global(mval, fval, halo_ps, lamult_da(1-nbdy,1-nbdy,l2ci))
-      call fill_global(mval, fval, halo_ps, lasl_da(1-nbdy,1-nbdy,l2ci))
+         call fill_global(mval, fval, halo_pv, util1)
+         call fill_global(mval, fval, halo_pv, util2)
+         call fill_global(mval, fval, halo_ps, lamult_da(1-nbdy,1-nbdy,l2ci))
+         call fill_global(mval, fval, halo_ps, lasl_da(1-nbdy,1-nbdy,l2ci))
 
-      call xctilr(util1, 1,1, 1,1, halo_pv)
-      call xctilr(util2, 1,1, 1,1, halo_pv)
+         call xctilr(util1, 1,1, 1,1, halo_pv)
+         call xctilr(util2, 1,1, 1,1, halo_pv)
 
-   !$omp parallel do private(l, i)
-      do j = 1, jj
-         do l = 1, isu(j)
-         do i = max(1,ifu(j,l)), min(ii,ilu(j,l))
-            ! x-component of surface Stokes drift [m s-1].
-            ustokes_da(i,j,l2ci) = .5_r8*(util1(i-1,j) + util1(i,j))
-         enddo
-         enddo
-         do l = 1,isv(j)
-         do i = max(1,ifv(j,l)), min(ii,ilv(j,l))
-            ! y-component of surface Stokes drift [m s-1].
-            vstokes_da(i,j,l2ci) = .5_r8*(util2(i,j-1) + util2(i,j))
-         enddo
+      !$omp parallel do private(l, i)
+         do j = 1, jj
+            do l = 1, isu(j)
+            do i = max(1,ifu(j,l)), min(ii,ilu(j,l))
+               ! x-component of surface Stokes drift [m s-1].
+               ustokes_da(i,j,l2ci) = .5_r8*(util1(i-1,j) + util1(i,j))
+            enddo
+            enddo
+            do l = 1,isv(j)
+            do i = max(1,ifv(j,l)), min(ii,ilv(j,l))
+               ! y-component of surface Stokes drift [m s-1].
+               vstokes_da(i,j,l2ci) = .5_r8*(util2(i,j-1) + util2(i,j))
+            enddo
+            enddo
          enddo
-      enddo
-   !$omp end parallel do
+      !$omp end parallel do
+
+      endif
 
 #ifdef PROGCO2
       call getfldindex(fldlist_num, fldlist, 'Sa_co2prog', index_Sa_co2prog)

meson.build

@@ -29,7 +29,7 @@ elif fcc.get_id() == 'intel'
                           language: 'fortran')
   endif
   if get_option('buildtype') == 'release'
-    add_project_arguments('-fp-model precise'.split(),language: 'fortran')
+    add_project_arguments('-O2 -fp-model precise'.split(),language: 'fortran')
   endif
   if get_option('processors') == 1 and get_option('grid') == 'channel'
     add_project_arguments('-mcmodel=medium', language: 'fortran')

phy/blom_init.F

@@ -30,7 +30,7 @@ subroutine blom_init
       use mod_xc
       use mod_pointtest, only: init_ptest
       use mod_inicon, only: icfile
-      use mod_state, only: dp, dpu, dpv, uflx, vflx, p, phi
+      use mod_state, only: dp, dpu, dpv, uflx, vflx, p, pu, pv, phi
       use mod_barotp, only: pvtrop
       use mod_pgforc, only: pgfxm, pgfym, xixp, xixm, xiyp, xiym
       use mod_niw, only: uml, vml, umlres, vmlres
@@ -49,7 +49,7 @@ subroutine blom_init
 c
       implicit none
 c
-      integer istat,ncid,varid,i,j,k,l,m,n,mm,nn,k1m,k1n,mt,mmt,km
+      integer istat,ncid,varid,i,j,k,l,m,n,mm,nn,k1m,k1n,mt,mmt,kn,km
       real q
       logical icrest,fexist
 c
@@ -274,12 +274,31 @@ subroutine blom_init
         call xctilr(dpu, 1,2*kk, 3,3, halo_us)
         call xctilr(dpv, 1,2*kk, 3,3, halo_vs)
 c
-c$OMP PARALLEL DO PRIVATE(k,l,i)
+c$OMP PARALLEL DO PRIVATE(k,km,l,i)
         do j=-2,jj+2
           do k=1,kk
+            km=k+mm
             do l=1,isp(j)
             do i=max(-2,ifp(j,l)),min(ii+2,ilp(j,l))
-              p(i,j,k+1)=p(i,j,k)+dp(i,j,k+mm)
+              p(i,j,k+1)=p(i,j,k)+dp(i,j,km)
+            enddo
+            enddo
+          enddo
+        enddo
+c$OMP END PARALLEL DO
+c
+c$OMP PARALLEL DO PRIVATE(k,kn,l,i)
+        do j=1,jj
+          do k=1,kk
+            kn=k+nn
+            do l=1,isu(j)
+            do i=max(1,ifu(j,l)),min(ii,ilu(j,l))
+              pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,kn)
+            enddo
+            enddo
+            do l=1,isv(j)
+            do i=max(1,ifv(j,l)),min(ii,ilv(j,l))
+              pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,kn)
             enddo
             enddo
           enddo