Changeset 13400


Ignore:
Timestamp:
2020-08-14T12:07:55+02:00 (7 months 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:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/cfgs/SHARED/namelist_ref

    r13398 r13400  
    10911091   ln_use_osm_la = .false.     !  Use   rn_osm_la 
    10921092   rn_osm_la     = 0.3         !  Turbulent Langmuir number 
    1093    rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
     1093   rn_zdfosm_adjust_sd = 1.0   ! Stokes drift reduction factor 
     1094   rn_osm_hblfrac = 0.1        ! specify top part of hbl for nn_osm_wave = 3 or 4 
    10941095   nn_ave = 0                  ! choice of horizontal averaging on avt, avmu, avmv 
    10951096   ln_dia_osm = .true.         ! output OSMOSIS-OBL variables 
     
    10991100   rn_difri  =  0.005          ! max Ri# diffusivity at Ri_g = 0 (m^2/s) 
    11001101   ln_convmix  = .true.        ! Use convective instability mixing below BL 
    1101    rn_difconv = 1.             ! diffusivity when unstable below BL  (m2/s) 
     1102   rn_difconv = 1. !0.01 !1.             ! diffusivity when unstable below BL  (m2/s) 
     1103   rn_osm_dstokes     = 5.     !  Depth scale of Stokes drift (m) 
    11021104   nn_osm_wave = 0             ! Method used to calculate Stokes drift 
    11031105      !                        !  = 2: Use ECMWF wave fields 
    11041106      !                        !  = 1: Pierson Moskowitz wave spectrum 
    11051107      !                        !  = 0: Constant La# = 0.3 
     1108   nn_osm_SD_reduce = 0        ! Method used to get active Stokes drift from surface value 
     1109      !                        !  = 0: No reduction 
     1110                               !  = 1: use SD avged over top 10% hbl 
     1111                               !  = 2:use surface value of SD fit to slope at rn_osm_hblfrac*hbl below surface 
     1112   ln_zdfosm_ice_shelter = .true.  ! reduce surface SD and depth scale under ice 
    11061113   ln_osm_mle = .false.        !  Use integrated FK-OSM model 
    11071114/ 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90

    r13399 r13400  
    103103   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number 
    104104   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift 
     105   REAL(wp) ::   rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by 
     106   REAL(wp) ::   rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl 
     107   LOGICAL  ::   ln_zdfosm_ice_shelter      ! flag to activate ice sheltering 
    105108   REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs 
    106109   INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt 
    107110   INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave 
     111   INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value 
    108112   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la 
    109  
    110113 
    111114   LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing 
     
    323326      REAL(wp) :: zzeta_v = 0.46 
    324327      REAL(wp) :: zabsstke 
     328      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness 
     329      REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc 
    325330 
    326331      ! For debugging 
     
    425430              zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
    426431              zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 
     432              ! Linearly 
    427433              zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 ) 
    428               ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init 
     434              dstokes(ji,jj) = rn_osm_dstokes 
    429435           END DO 
    430436        END DO 
     
    434440           DO ji = 2, jpim1 
    435441              ! Use wind speed wndm included in sbc_oce module 
    436               zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
    437               dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
     442              zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     443              dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
    438444           END DO 
    439445        END DO 
     
    441447     CASE(2) 
    442448        zfac =  2.0_wp * rpi / 16.0_wp 
     449 
    443450        DO jj = 2, jpjm1 
    444451           DO ji = 2, jpim1 
    445               ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 
    446               !    The coefficient 0.8 gives La=0.3  in this situation. 
    447               ! It could represent the effects of the spread of wave directions 
    448               ! around the mean wind. The effect of this adjustment needs to be tested. 
    449               zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
    450               zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
    451               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 ! 
     452              IF (hsw(ji,jj) > 1.e-4) THEN 
     453                 ! Use  wave fields 
     454                 zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2) 
     455                 zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8) 
     456                 dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) 
     457              ELSE 
     458                 ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run) 
     459                 ! .. so default to Pierson-Moskowitz 
     460                 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) 
     461                 dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1) 
     462              END IF 
     463            END DO 
     464        END DO 
     465     END SELECT 
     466 
     467     IF (ln_zdfosm_ice_shelter) THEN 
     468        ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice 
     469        DO jj = 2, jpjm1 
     470           DO ji = 2, jpim1 
     471              zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj)) 
     472              dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj)) 
     473           END DO 
     474        END DO 
     475     END IF 
     476 
     477     SELECT CASE (nn_osm_SD_reduce) 
     478     ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van  Roekel (2012) or Grant (2020). 
     479     CASE(0) 
     480        ! The Langmur number from the ECMWF model (or from PM)  appears to give La<0.3 for wind-driven seas. 
     481        !    The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3  in this situation. 
     482        ! It could represent the effects of the spread of wave directions 
     483        ! around the mean wind. The effect of this adjustment needs to be tested. 
     484        IF(nn_osm_wave > 0) THEN 
     485           zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1) 
     486        END IF 
     487     CASE(1) 
     488        ! van  Roekel (2012): consider average SD over top 10% of boundary layer 
     489        ! assumes approximate depth profile of SD from Breivik (2016) 
     490        zsqrtpi = SQRT(rpi) 
     491        z_two_thirds = 2.0_wp / 3.0_wp 
     492 
     493        DO jj = 2, jpjm1 
     494           DO ji = 2, jpim1 
     495              zthickness = rn_osm_hblfrac*hbl(ji,jj) 
     496              z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 
     497              zsqrt_depth = SQRT(z2k_times_thickness) 
     498              zexp_depth  = EXP(-z2k_times_thickness) 
     499              zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth  & 
     500                   &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) & 
     501                   &              + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness 
     502 
     503           END DO 
     504        END DO 
     505     CASE(2) 
     506        ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer 
     507        ! assumes approximate depth profile of SD from Breivik (2016) 
     508        zsqrtpi = SQRT(rpi) 
     509 
     510        DO jj = 2, jpjm1 
     511           DO ji = 2, jpim1 
     512              zthickness = rn_osm_hblfrac*hbl(ji,jj) 
     513              z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp ) 
     514 
     515              IF(z2k_times_thickness < 50._wp) THEN 
     516                 zsqrt_depth = SQRT(z2k_times_thickness) 
     517                 zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness) 
     518              ELSE 
     519                 ! asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness 
     520                 ! See Abramowitz and Stegun, Eq. 7.1.23 
     521                 ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness)  + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3) 
     522                 zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp 
     523              END IF 
     524              zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp) 
     525              dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj) 
     526              zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc) 
    452527           END DO 
    453528        END DO 
     
    488563     ! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must 
    489564     ! previously exist for hbl also. 
     565 
     566     ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway 
     567     ! ########################################################################## 
    490568      hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) ) 
    491569      ibld(:,:) = 4 
     
    499577         END DO 
    500578      END DO 
     579     ! ########################################################################## 
    501580 
    502581      DO jj = 2, jpjm1 
     
    11651244            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
    11661245         ! Stokes drift read in from sbcwave  (=2). 
    1167          CASE(2) 
     1246         CASE(2:3) 
    11681247            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )               ! x surface Stokes drift 
    11691248            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )               ! y surface Stokes drift 
     
    15201599             zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 ) 
    15211600             zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 ) 
    1522              zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
    1523                   &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1601             IF (nn_osm_SD_reduce > 0 ) THEN 
     1602                ! Effective Stokes drift already reduced from surface value 
     1603                zr_stokes = 1.0_wp 
     1604             ELSE 
     1605                ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd, 
     1606                ! requires further reduction where BL is deep 
     1607                zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) & 
     1608                     &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) ) 
     1609             END IF 
    15241610 
    15251611             zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) & 
     
    20002086     !! 
    20012087     NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave & 
    2002           & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 & 
     2088          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd & 
    20032089          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave & 
    2004           & ,ln_osm_mle 
     2090          & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, ln_zdfosm_ice_shelter 
    20052091! Namelist for Fox-Kemper parametrization. 
    20062092      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,& 
     
    20262112        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle 
    20272113        WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la 
     2114        WRITE(numout,*) '     Stokes drift reduction factor                 rn_zdfosm_adjust_sd   = ', rn_zdfosm_adjust_sd 
    20282115        WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0 
    20292116        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes 
     
    20372124        CASE(2) 
    20382125           WRITE(numout,*) '     calculated from ECMWF wave fields' 
     2126         END SELECT 
     2127        WRITE(numout,*) '     Stokes drift reduction                        nn_osm_SD_reduce', nn_osm_SD_reduce 
     2128        WRITE(numout,*) '     fraction of hbl to average SD over/fit' 
     2129        WRITE(numout,*) '     exponential with nn_osm_SD_reduce = 1 or 2    rn_osm_hblfrac =  ', rn_osm_hblfrac 
     2130        SELECT CASE (nn_osm_SD_reduce) 
     2131        CASE(0) 
     2132           WRITE(numout,*) '     No reduction' 
     2133        CASE(1) 
     2134           WRITE(numout,*) '     Average SD over upper rn_osm_hblfrac of BL' 
     2135        CASE(2) 
     2136           WRITE(numout,*) '     Fit exponential to slope rn_osm_hblfrac of BL' 
    20392137        END SELECT 
     2138        WRITE(numout,*) '     reduce surface SD and depth scale under ice   ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter 
    20402139        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm 
    20412140        WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix 
Note: See TracChangeset for help on using the changeset viewer.