Changeset 13400 for NEMO/branches/2019

Timestamp:

2020-08-14T12:07:55+02:00 (4 years ago)

Author:

agn

Message:

New code to reduce Stokes drift and allow ice sheltering

ln_zdfosm_ice_shelter if true=> ice sheltering
rn_zdfosm_adjust_sd sets factor to reduce SD by if not using new van Roekel method
nn_osm_SD_reduce if 0 => use rn_zdfosm_adjust_sd

if 1 => reduce SD to hbl/10 average
if 2 => ewduce SD to matched exp @ hbl/10

Location:

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0

Files:

• cfgs/SHARED/namelist_ref (modified) (2 diffs)
• src/OCE/ZDF/zdfosm.F90 (modified) (12 diffs)

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/cfgs/SHARED/namelist_ref

@@ -1091 +1091 @@
    ln_use_osm_la = .false.     !  Use   rn_osm_la
    rn_osm_la     = 0.3         !  Turbulent Langmuir number
-   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m)
+   rn_zdfosm_adjust_sd = 1.0   ! Stokes drift reduction factor
+   rn_osm_hblfrac = 0.1        ! specify top part of hbl for nn_osm_wave = 3 or 4
    nn_ave = 0                  ! choice of horizontal averaging on avt, avmu, avmv
    ln_dia_osm = .true.         ! output OSMOSIS-OBL variables
@@ -1099 +1100 @@
    rn_difri  =  0.005          ! max Ri# diffusivity at Ri_g = 0 (m^2/s)
    ln_convmix  = .true.        ! Use convective instability mixing below BL
-   rn_difconv = 1.             ! diffusivity when unstable below BL  (m2/s)
+   rn_difconv = 1. !0.01 !1.             ! diffusivity when unstable below BL  (m2/s)
+   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m)
    nn_osm_wave = 0             ! Method used to calculate Stokes drift
       !                        !  = 2: Use ECMWF wave fields
       !                        !  = 1: Pierson Moskowitz wave spectrum
       !                        !  = 0: Constant La# = 0.3
+   nn_osm_SD_reduce = 0        ! Method used to get active Stokes drift from surface value
+      !                        !  = 0: No reduction
+                               !  = 1: use SD avged over top 10% hbl
+                               !  = 2:use surface value of SD fit to slope at rn_osm_hblfrac*hbl below surface
+   ln_zdfosm_ice_shelter = .true.  ! reduce surface SD and depth scale under ice
    ln_osm_mle = .false.        !  Use integrated FK-OSM model
 /

NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90

@@ -103 +103 @@
    REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number
    REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift
+   REAL(wp) ::   rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by
+   REAL(wp) ::   rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl
+   LOGICAL  ::   ln_zdfosm_ice_shelter      ! flag to activate ice sheltering
    REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs
    INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt
    INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
+   INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value
    LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
-
 
    LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing
@@ -323 +326 @@
       REAL(wp) :: zzeta_v = 0.46
       REAL(wp) :: zabsstke
+      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness
+      REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc
 
       ! For debugging
@@ -425 +430 @@
               zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
               zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
+              ! Linearly
               zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 )
-              ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
+              dstokes(ji,jj) = rn_osm_dstokes
            END DO
         END DO
@@ -434 +440 @@
            DO ji = 2, jpim1
               ! Use wind speed wndm included in sbc_oce module
-              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
-              dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
+              zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
+              dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
            END DO
         END DO
@@ -441 +447 @@
      CASE(2)
         zfac =  2.0_wp * rpi / 16.0_wp
+
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
-              ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
-              !    The coefficient 0.8 gives La=0.3  in this situation.
-              ! It could represent the effects of the spread of wave directions
-              ! around the mean wind. The effect of this adjustment needs to be tested.
-              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)
-              zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8)
-              dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
+              IF (hsw(ji,jj) > 1.e-4) THEN
+                 ! Use  wave fields
+                 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)
+                 zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8)
+                 dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1)
+              ELSE
+                 ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run)
+                 ! .. so default to Pierson-Moskowitz
+                 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
+                 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
+              END IF
+            END DO
+        END DO
+     END SELECT
+
+     IF (ln_zdfosm_ice_shelter) THEN
+        ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj))
+              dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj))
+           END DO
+        END DO
+     END IF
+
+     SELECT CASE (nn_osm_SD_reduce)
+     ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van  Roekel (2012) or Grant (2020).
+     CASE(0)
+        ! The Langmur number from the ECMWF model (or from PM)  appears to give La<0.3 for wind-driven seas.
+        !    The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3  in this situation.
+        ! It could represent the effects of the spread of wave directions
+        ! around the mean wind. The effect of this adjustment needs to be tested.
+        IF(nn_osm_wave > 0) THEN
+           zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1)
+        END IF
+     CASE(1)
+        ! van  Roekel (2012): consider average SD over top 10% of boundary layer
+        ! assumes approximate depth profile of SD from Breivik (2016)
+        zsqrtpi = SQRT(rpi)
+        z_two_thirds = 2.0_wp / 3.0_wp
+
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              zthickness = rn_osm_hblfrac*hbl(ji,jj)
+              z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
+              zsqrt_depth = SQRT(z2k_times_thickness)
+              zexp_depth  = EXP(-z2k_times_thickness)
+              zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth  &
+                   &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) &
+                   &              + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness
+
+           END DO
+        END DO
+     CASE(2)
+        ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer
+        ! assumes approximate depth profile of SD from Breivik (2016)
+        zsqrtpi = SQRT(rpi)
+
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              zthickness = rn_osm_hblfrac*hbl(ji,jj)
+              z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
+
+              IF(z2k_times_thickness < 50._wp) THEN
+                 zsqrt_depth = SQRT(z2k_times_thickness)
+                 zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness)
+              ELSE
+                 ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness
+                 ! See Abramowitz and Stegun, Eq. 7.1.23
+                 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness)  + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3)
+                 zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp
+              END IF
+              zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp)
+              dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj)
+              zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc)
            END DO
         END DO
@@ -488 +563 @@
      ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must
      ! previously exist for hbl also.
+
+     ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway
+     ! ##########################################################################
       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
       ibld(:,:) = 4
@@ -499 +577 @@
          END DO
       END DO
+     ! ##########################################################################
 
       DO jj = 2, jpjm1
@@ -1165 +1244 @@
             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke )
          ! Stokes drift read in from sbcwave  (=2).
-         CASE(2)
+         CASE(2:3)
             IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift
             IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift
@@ -1520 +1599 @@
              zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
              zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
-             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
-                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
+             IF (nn_osm_SD_reduce > 0 ) THEN
+                ! Effective Stokes drift already reduced from surface value
+                zr_stokes = 1.0_wp
+             ELSE
+                ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
+                ! requires further reduction where BL is deep
+                zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
+                     &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
+             END IF
 
              zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
@@ -2000 +2086 @@
      !!
      NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
-          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
+          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd &
           & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
-          & ,ln_osm_mle
+          & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, ln_zdfosm_ice_shelter
 ! Namelist for Fox-Kemper parametrization.
       NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
@@ -2026 +2112 @@
         WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle
         WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
+        WRITE(numout,*) '     Stokes drift reduction factor                 rn_zdfosm_adjust_sd   = ', rn_zdfosm_adjust_sd
         WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
         WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes
@@ -2037 +2124 @@
         CASE(2)
            WRITE(numout,*) '     calculated from ECMWF wave fields'
+         END SELECT
+        WRITE(numout,*) '     Stokes drift reduction                        nn_osm_SD_reduce', nn_osm_SD_reduce
+        WRITE(numout,*) '     fraction of hbl to average SD over/fit'
+        WRITE(numout,*) '     exponential with nn_osm_SD_reduce = 1 or 2    rn_osm_hblfrac =  ', rn_osm_hblfrac
+        SELECT CASE (nn_osm_SD_reduce)
+        CASE(0)
+           WRITE(numout,*) '     No reduction'
+        CASE(1)
+           WRITE(numout,*) '     Average SD over upper rn_osm_hblfrac of BL'
+        CASE(2)
+           WRITE(numout,*) '     Fit exponential to slope rn_osm_hblfrac of BL'
         END SELECT
+        WRITE(numout,*) '     reduce surface SD and depth scale under ice   ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter
         WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
         WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix