Changeset 12532


Ignore:
Timestamp:
2020-03-11T10:26:18+01:00 (8 months ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0.1_biharmonic_GM: science changes.

Location:
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/cfgs/SHARED/namelist_ref

    r12083 r12532  
    923923   ln_dynldf_lap = .false.     !    laplacian operator 
    924924   ln_dynldf_blp = .false.     !  bilaplacian operator 
     925   ln_dynldf_bgm = .false.     !  bilaplician Gent-McWilliams  
    925926   !                       !  Direction of action  : 
    926927   ln_dynldf_lev = .false.     !  iso-level 
     
    947948      !                       !  iso-neutral laplacian operator (ln_dynldf_iso=T) : 
    948949      rn_ahm_b      = 0.0         !  background eddy viscosity  [m2/s] 
     950   nn_bhm_ijk_t  = 0           !  space/time variation of bilaplacian GM coefficient : 
     951      !                             !  = 11  uniform coefficient in interior, zero top and bottom 
     952      !                             !  = 12  linear tapering with dz to surface, zero at top and bottom 
     953      !                             !  = 13  dx^2.dz^2 scaling, zero at top and bottom 
     954      rn_bhm_b      = 0.0         ! bilaplacian GM coefficient 
     955      rn_bgmzcrit   = 0.0         ! critical depth for depth variation of bilaplacian GM coefficient 
    949956/ 
    950957!----------------------------------------------------------------------- 
  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90

    r11715 r12532  
    4242   !!---------------------------------------------------------------------- 
    4343   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    44    !! $Id$ 
     44   !! $Id: diaar5.F90 10425 2018-12-19 21:54:16Z smasson $ 
    4545   !! Software governed by the CeCILL license (see ./LICENSE) 
    4646   !!---------------------------------------------------------------------- 
     
    136136         !                                         ! steric sea surface height 
    137137         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density 
    138          zrhop(:,:,jpk) = 0._wp 
     138         !CW: commenting out the following line, as seems weird to zero the array before output & only way to get rhop output 
     139         !zrhop(:,:,jpk) = 0._wp 
    139140         CALL iom_put( 'rhop', zrhop ) 
    140141         ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90

    r11715 r12532  
    77   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
    88   !!                 !                                  add velocity dependent coefficient and optional read in file 
     9   !!            4.0  ! 2020-02  (C. Wilson, ...) add bi-Laplacian GM implementation via momentum      
    910   !!---------------------------------------------------------------------- 
    1011 
     
    3839   !!---------------------------------------------------------------------- 
    3940   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    40    !! $Id$ 
     41   !! $Id: dynldf.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ 
    4142   !! Software governed by the CeCILL license (see ./LICENSE) 
    4243   !!---------------------------------------------------------------------- 
     
    6667      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
    6768      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian 
     69      CASE ( np_bgm )  ;       CALL dyn_ldf_bgm( kt, ub, vb, ua, va)  ! bi-Laplacian GM parameterisation 
    6870      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 
    6971      ! 
     
    104106         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator' 
    105107         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background' 
     108         CASE( np_bgm )      ;   WRITE(numout,*) '   ==>>>   bi-Laplacian GM parameterisation via momentum equation' 
    106109         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator' 
    107110         END SELECT 
  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90

    r11715 r12532  
    55   !!====================================================================== 
    66   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
     7   !!           4.0  ! 2020-02  (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum      
    78   !!---------------------------------------------------------------------- 
    89 
     
    2526   PUBLIC dyn_ldf_lap  ! called by dynldf.F90 
    2627   PUBLIC dyn_ldf_blp  ! called by dynldf.F90 
     28   PUBLIC dyn_ldf_bgm  ! called by dynldf.F90 
     29   PRIVATE dyn_ldf_lap_no_ahm !needed by dyn_ldf_bgm 
    2730 
    2831   !! * Substitutions 
     
    3033   !!---------------------------------------------------------------------- 
    3134   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    32    !! $Id$ 
     35   !! $Id: dynldf_lap_blp.F90 10425 2018-12-19 21:54:16Z smasson $  
    3336   !! Software governed by the CeCILL license (see ./LICENSE) 
    3437   !!---------------------------------------------------------------------- 
     
    104107   END SUBROUTINE dyn_ldf_lap 
    105108 
     109!CW: new subroutine for the Laplacian of velocity, copied from dyn_ldf_lap above, but without eddy viscosity ahmf, ahmt 
     110 
     111   SUBROUTINE dyn_ldf_lap_no_ahm( kt, pub, pvb, pulap, pvlap, kpass ) 
     112      !!---------------------------------------------------------------------- 
     113      !!                     ***  ROUTINE dyn_ldf_lap_no_ahm  *** 
     114      !!                        
     115      !! ** Purpose :   Companion subroutine to dyn_ldf_bgm to calculate  
     116      !!      the before horizontal momentum Laplacian  
     117      !!      trend and add it to the general trend of momentum equation.   
     118      !!      Note: mimics dyn_ldf_lap but without eddy viscosity ahmf, ahmt,  
     119      !!      because biharmonic GM eddy diffusivity is applied in dyn_bgm. 
     120      !! 
     121      !! ** Method  :   The Laplacian operator apply on horizontal velocity is  
     122      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
     123      !! 
     124      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. 
     125      !!---------------------------------------------------------------------- 
     126      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     127      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     128      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s] 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) ::   pulap, pvlap ! velocity trend   [m/s2] 
     130      ! 
     131      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     132      REAL(wp) ::   zsign        ! local scalars 
     133      REAL(wp) ::   zua, zva     ! local scalars 
     134      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv 
     135      !!---------------------------------------------------------------------- 
     136      ! 
     137      IF( kt == nit000 .AND. lwp ) THEN 
     138         WRITE(numout,*) 
     139         WRITE(numout,*) 'dyn_ldf_lap_no_ahm : companion iso-level harmonic (laplacian) operator to dyn_bgm, pass=', kpass 
     140         WRITE(numout,*) '~~~~~~~ ' 
     141      ENDIF 
     142      ! 
     143      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign 
     144      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0) 
     145      ENDIF 
     146      ! 
     147      !                                                ! =============== 
     148      DO jk = 1, jpkm1                                 ! Horizontal slab 
     149         !                                             ! =============== 
     150         DO jj = 2, jpj 
     151            DO ji = fs_2, jpi   ! vector opt. 
     152               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
     153!!gm open question here : e3f  at before or now ?    probably now... 
     154 
     155!!CW: note that the removed ahmf had already been multiplied by fmask, so we multiply by fmask explicitly here, 
     156!!    with the i,j,k of fmask aligning with those of ahmf, as defined in ldfdyn.F90. 
     157               zcur(ji-1,jj-1) = fmask(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      & 
     158                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  & 
     159                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) 
     160               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
     161!!CW: note that the removed ahmt had already been multiplied by tmask, so we multiply by tmask explicitly here, 
     162!!    with the i,j,k of tmask aligning with those of ahmt, as defined in ldfdyn.F90 
     163               zdiv(ji,jj)     = tmask(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         & 
     164                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  & 
     165                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  ) 
     166            END DO   
     167         END DO   
     168         ! 
     169!CW: multiply pulap, pvlap by umask, vmask 
     170         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div ) 
     171            DO ji = fs_2, fs_jpim1   ! vector opt. 
     172               pulap(ji,jj,jk) = umask(ji,jj,jk) * zsign * (                                            & 
     173                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   & 
     174                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
     175                  ! 
     176               pvlap(ji,jj,jk) = vmask(ji,jj,jk) * zsign * (                                             & 
     177                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   & 
     178                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
     179            END DO 
     180         END DO 
     181         !                                             ! =============== 
     182      END DO                                           !   End of slab 
     183      !                                                ! =============== 
     184      ! 
     185   END SUBROUTINE dyn_ldf_lap_no_ahm 
     186 
     187 
     188 
     189 
     190!CW: new subroutine for biharmonic GM, copied from SUBROUTINE dyn_ldf_blp below 
     191 
     192   SUBROUTINE dyn_ldf_bgm( kt, pub, pvb, pua, pva ) 
     193      !!---------------------------------------------------------------------- 
     194      !!                 ***  ROUTINE dyn_bgm  *** 
     195      !!                     
     196      !! ** Purpose :   Compute the before lateral momentum trend due to the bi-Laplacian GM parameterisation  
     197      !!              and add it to the general trend of momentum equation. 
     198      !! 
     199      !! ** Method  :   The bi-Laplacian implementation of GM is via a -d/dz(diffusivity x d/dz(Laplacian of velocity)) 
     200      !!      operator applied at the 'now' time level.  The existing code for the Laplacian contains the 'before' time also in zdiv. 
     201      !!      It is computed by a call to dyn_ldf_lap routine and vertical differentiation applied twice. 
     202      !! 
     203      !! ** Action :   pua, pva increased with the before bi-Laplacian GM momentum trend calculated from pub, pvb. 
     204      !!---------------------------------------------------------------------- 
     205      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     206      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields 
     207      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend 
     208      !    
     209      INTEGER  ::   iku, ikv                         ! local integers 
     210      !CW 
     211      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     212      ! 
     213      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point 
     214      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz   ! -1*bhm * d/dz(del^2 u) at u- and v-point 
     215      !!---------------------------------------------------------------------- 
     216      ! 
     217      IF( kt == nit000 )  THEN 
     218         IF(lwp) WRITE(numout,*) 
     219         IF(lwp) WRITE(numout,*) 'dyn_bgm : bi-Laplacian GM operator via momentum ' 
     220         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     221      ENDIF 
     222      ! 
     223   
     224      CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 )    
     225 
     226      ! =============== 
     227      !CW: calculate -bhm * d/dz(del^2 u)   
     228      DO jk = 2, jpkm1  
     229         DO jj = 2, jpjm1 
     230            DO ji = fs_2, jpim1   ! vector opt. 
     231                zulapdz(ji,jj,jk) = -bhmu(ji,jj,jk)*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) / e3uw_n(ji,jj,jk) 
     232                zvlapdz(ji,jj,jk) = -bhmv(ji,jj,jk)*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) / e3vw_n(ji,jj,jk) 
     233            ENDDO 
     234         ENDDO 
     235      ENDDO 
     236 
     237 
     238!CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0 
     239!     Surface 
     240      zulapdz(:,:,1)   = 0._wp   ;   zvlapdz(:,:,1)   = 0._wp  
     241!     Flat bottom case 
     242      zulapdz(:,:,jpk) = 0._wp   ;   zvlapdz(:,:,jpk) = 0._wp 
     243!     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 
     244        DO jj = 2, jpjm1 
     245            DO ji = 2, jpim1 
     246                iku = mbku(ji,jj)+1 
     247                ikv = mbkv(ji,jj)+1 
     248                zulapdz(:,:,iku) = 0._wp 
     249                zvlapdz(:,:,ikv) = 0._wp 
     250            ENDDO 
     251        ENDDO 
     252 
     253 
     254 
     255 
     256!!    calculate d/dz(-bhm * d/dz(del^2 u)) 
     257                                                       ! =============== 
     258      DO jk = 1, jpkm1                                 ! Horizontal slab 
     259         !                                             ! =============== 
     260         DO jj = 2, jpjm1 
     261            DO ji = fs_2, jpim1   ! vector opt. 
     262            pua(ji,jj,jk) = pua(ji,jj,jk) + (zulapdz(ji,jj,jk) - zulapdz(ji,jj,jk+1)) / e3u_n(ji,jj,jk) 
     263            pva(ji,jj,jk) = pva(ji,jj,jk) + (zvlapdz(ji,jj,jk) - zvlapdz(ji,jj,jk+1)) / e3v_n(ji,jj,jk) 
     264            ENDDO 
     265         ENDDO 
     266      ENDDO 
     267 
     268      ! ----- 
     269 
     270    
     271   END SUBROUTINE dyn_ldf_bgm 
     272    
    106273 
    107274   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) 
  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90

    r11715 r12532  
    88   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
    99   !!                 !                                  add velocity dependent coefficient and optional read in file 
     10   !!            4.0  ! 2020-02  (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum      
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    1516   !!---------------------------------------------------------------------- 
    1617   USE oce             ! ocean dynamics and tracers    
    17    USE dom_oce         ! ocean space and time domain  
     18   USE dom_oce         ! ocean space and time domain 
    1819   USE phycst          ! physical constants 
    1920   USE ldfslp          ! lateral diffusion: slopes of mixing orientation 
     
    3536   LOGICAL , PUBLIC ::   ln_dynldf_OFF   !: No operator (i.e. no explicit diffusion) 
    3637   LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator 
     38   LOGICAL , PUBLIC ::   ln_dynldf_bgm   !: bi-Laplacian GM implementation via momentum 
    3739   LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator 
    3840   LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction 
     
    4244   !                                        !  time invariant coefficients:  aht = 1/2  Ud*Ld   (lap case)  
    4345   !                                           !                             bht = 1/12 Ud*Ld^3 (blp case) 
     46   INTEGER , PUBLIC ::   nn_bhm_ijk_t    !: choice of time & space variations of the lateral bi-Laplacian GM coef. 
    4447   REAL(wp), PUBLIC ::   rn_Uv                 !: lateral viscous velocity  [m/s] 
    4548   REAL(wp), PUBLIC ::   rn_Lv                 !: lateral viscous length    [m] 
     
    5053   !                                        ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 
    5154   REAL(wp), PUBLIC ::   rn_ahm_b              !: lateral laplacian background eddy viscosity  [m2/s] 
    52  
     55   REAL(wp), PUBLIC ::   rn_bhm_b              !: lateral bi-Laplacian GM background eddy viscosity  [m4/s] 
     56   REAL(wp), PUBLIC ::   rn_bgmzcrit           !: critical depth (>=0) for linear tapering of bi-Lap. GM coef. to zero at surface [m]  
    5357   !                                    !!* Parameter to control the type of lateral viscous operator 
    5458   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator 
     
    5761   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator 
    5862   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator 
    59    ! 
     63   INTEGER, PARAMETER, PUBLIC ::   np_bgm    = 30                       !: iso-level operator, bi-Laplacian GM 
    6064   INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
    6165   LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef. 
    6266 
    6367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 
     68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bhmu, bhmv   !: eddy viscosity coef. for bi-Laplacian GM at u- and v-points [m4/s] 
    6469   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only) 
    6570   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only) 
     
    7681   !!---------------------------------------------------------------------- 
    7782   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    78    !! $Id$ 
     83   !! $Id: ldfdyn.F90 11653 2019-10-04 12:34:02Z davestorkey $  
    7984   !! Software governed by the CeCILL license (see ./LICENSE) 
    8085   !!---------------------------------------------------------------------- 
     
    9196      !!             ln_dynldf_lap = T     laplacian operator 
    9297      !!             ln_dynldf_blp = T   bilaplacian operator 
     98      !!             ln_dynldf_bgm = T   bi-Laplacian GM operator 
    9399      !!              - the parameter nn_ahm_ijk_t: 
    94100      !!    nn_ahm_ijk_t  =  0 => = constant 
     
    103109      !!                                                           (   L^2|D|      laplacian operator 
    104110      !!                                                           or  L^4|D|/8  bilaplacian operator ) 
     111      !!                - the parameter nn_bhm_ijk_t: 
     112      !!    nn_bhm_ijk_t  = 11 => = F(z) :     = constant, with BCs set to zero at top and bottom 
     113      !!    nn_bhm_ijk_t  = 12 => = F(z) :     = constant in interior, linear taper to zero at surface, 
     114      !!                                         for depths < rn_bgmzcrit, with BCs set to zero at top and bottom 
     115      !!    nn_bhm_ijk_t  = 13 => = F(z) :     = bhm0 X dx^2 X dz^2, with BCs set to zero at top and bottom 
    105116      !!---------------------------------------------------------------------- 
    106117      INTEGER  ::   ji, jj, jk                     ! dummy loop indices 
    107118      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer 
     119      INTEGER  ::   iku, ikv, ikt                  ! local integer 
    108120      REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar 
     121      REAL(wp) ::   bhm0, lhscrit, lhscritmax      ! local scalar 
    109122      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s) 
    110123      !! 
    111124      NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator 
     125         &                 ln_dynldf_bgm,                                 &   ! type of operator 
    112126         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator 
    113127         &                 nn_ahm_ijk_t , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient 
     128         &                 nn_bhm_ijk_t, rn_bhm_b, rn_bgmzcrit,           &   ! lateral bi-Lap. GM coefficient 
    114129         &                 rn_csmc      , rn_minfac    , rn_maxfac            ! Smagorinsky settings 
    115130      !!---------------------------------------------------------------------- 
     
    134149         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap 
    135150         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp 
     151         WRITE(numout,*) '         bi-Laplacian GM operator             ln_dynldf_bgm = ', ln_dynldf_bgm 
     152         WRITE(numout,*) '         critical depth for taper (if used)   rn_bgmzcrit   = ', rn_bgmzcrit 
    136153         ! 
    137154         WRITE(numout,*) '      direction of action :' 
     
    145162         WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m' 
    146163         WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s' 
    147          ! 
     164         WRITE(numout,*) '         bi-Laplacian GM background viscosity    rn_bhm_b   = ', rn_bhm_b, ' m4/s' 
    148165         WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :' 
    149166         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc 
     
    159176      nldf_dyn = np_ERROR 
    160177      ioptio = 0 
     178      !CW exclude combination of lap+blp (as in existing code), but allow other combinations 
    161179      IF( ln_dynldf_OFF ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF 
    162       IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
    163       IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
    164       IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
    165       ! 
     180      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF 
     181      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF 
     182      IF( ln_dynldf_bgm ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
     183      ! 
     184      IF( (ioptio < 1).OR.(ioptio > 6) )   CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp/bgm) but not (lap+blp)' ) 
     185      ! 
     186 
    166187      IF(.NOT.ln_dynldf_OFF ) THEN     !==  direction ==>> type of operator  ==! 
    167188         ioptio = 0 
     
    209230         ENDIF 
    210231         ! 
     232         IF( ln_dynldf_bgm ) THEN         ! bi-Laplacian GM operator 
     233            IF( ln_zco ) THEN                   ! z-coordinate 
     234               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level = horizontal (no rotation) 
     235               IF( ln_dynldf_hor )   nldf_dyn = np_bgm   ! iso-level = horizontal (no rotation) 
     236               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation) 
     237            ENDIF 
     238            IF( ln_zps ) THEN                   ! z-coordinate with partial step 
     239               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level              (no rotation) 
     240               IF( ln_dynldf_hor )   nldf_dyn = np_bgm   ! iso-level              (no rotation) 
     241               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation) 
     242            ENDIF 
     243            IF( ln_sco ) THEN                   ! s-coordinate 
     244               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level              (no rotation) 
     245               IF( ln_dynldf_hor )   ierr = 3            ! horizontal             (   rotation) 
     246               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation) 
     247            ENDIF 
     248         ENDIF 
     249         !          
     250 
    211251         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 
    212252         ! 
     253         IF( ierr == 3 )   CALL ctl_stop( 'rotated bi-Laplacian GM operator does not exist' ) 
     254         ! 
     255 
    213256         IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes 
    214257         ! 
     
    222265         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background' 
    223266         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator' 
     267         CASE( np_bgm    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-Laplacian GM operator' 
    224268         END SELECT 
    225269         WRITE(numout,*) 
     
    233277      ! 
    234278      IF( ln_dynldf_OFF ) THEN 
    235          IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated' 
     279         !CW 
     280         !IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated' 
     281         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt, ahmf, bhmu, bhmv are not allocated' 
     282         ! 
    236283         RETURN 
    237284         ! 
     
    239286         ! 
    240287         !                                         ! allocate the ahm arrays 
    241          ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 
     288         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , bhmu(jpi,jpj,jpk) , bhmv(jpi,jpj,jpk) , STAT=ierr ) 
     289         ! 
    242290         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 
    243291         ! 
    244292         ahmt(:,:,:) = 0._wp                       ! init to 0 needed  
    245293         ahmf(:,:,:) = 0._wp 
    246          ! 
    247          !                                         ! value of lap/blp eddy mixing coef. 
     294         bhmu(:,:,:) = 0._wp                       ! init to 0 needed  
     295         bhmv(:,:,:) = 0._wp 
     296         !                                        ! value of lap/blp eddy mixing coef. 
    248297         IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian 
    249298         ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian 
     299         ELSEIF( ln_dynldf_bgm ) THEN   ;   bhm0 = rn_bhm_b       ;             ;   cl_Units = ' m4/s'   ! bi-Laplacian GM 
    250300         ENDIF 
    251301         zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient 
     
    265315            ahmf(:,:,1) = zah0 
    266316            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 
    267             ! 
     317            !        
    268318         CASE ( -20 )      !== fixed horizontal shape read in file  ==! 
    269319            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file' 
     
    323373            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 
    324374         END SELECT 
     375!CW      Add separate structure options for bi-Laplacian GM, to allow it to be used in combination with other types of diffusion, above. 
     376 
     377         IF( ln_dynldf_bgm ) THEN 
     378            SELECT CASE (nn_bhm_ijk_t) 
     379               !CW Define bi-Laplacian GM diffusivity and test the related computational stability criterion    
     380 
     381            CASE(  11  )      !==  fixed profile: constant except for zero at top and bottom, so that eddy-induced velocity, w*=0  ==! 
     382               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity = constant in interior, zero at top and bottom' 
     383               IF(lwp) WRITE(numout,*) '           interior viscous coef. = constant = ', bhm0, cl_Units 
     384 
     385               !     First set to constant in interior, all levels 
     386               DO jj = 2, jpjm1 
     387                  DO ji = 2, jpim1 
     388                     bhmu(ji,jj,:) = bhm0 
     389                     bhmv(ji,jj,:) = bhm0 
     390                  ENDDO 
     391               ENDDO 
     392 
     393               !     Surface BC : set to zero 
     394               bhmu(:,:,1)   = 0._wp   ;   bhmv(:,:,1)   = 0._wp 
     395               !     Flat bottom BC : set to zero 
     396               bhmu(:,:,jpk) = 0._wp   ;   bhmv(:,:,jpk) = 0._wp 
     397               !     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 
     398               !     and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp 
     399               DO jj = 2, jpjm1 
     400                  DO ji = 2, jpim1 
     401                     iku = mbku(ji,jj)+1 
     402                     ikv = mbkv(ji,jj)+1 
     403                     bhmu(:,:,iku) = 0._wp 
     404                     bhmv(:,:,ikv) = 0._wp 
     405                  ENDDO 
     406               ENDDO 
     407 
     408               !CW Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping 
     409               ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32  
     410               ! test at T-point - may need to revise this choice!!!! 
     411 
     412               ! Find domain maximum value of LHS, then test inequality. 
     413 
     414               ! initialise 
     415               lhscritmax=0._wp 
     416 
     417               DO jj = 2, jpjm1 
     418                  DO ji = 2, jpim1 
     419                     ikt = mbkt(ji,jj) !bottom last wet T-level 
     420                     DO jk = 2, ikt 
     421                        lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2) 
     422                        IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit 
     423                     ENDDO 
     424                  ENDDO 
     425               ENDDO 
     426 
     427               IF ( lhscrit .ge. 1._wp/32._wp) THEN 
     428                  IF(lwp) THEN 
     429                     WRITE(numout,*) 
     430                     WRITE(numout,*) ' WARNING : Bi-Laplacian GM eddy viscosity is not',              & 
     431                          &            ' consistent with the model time step and grid setup: ', & 
     432                          &            ' LHS=',lhscrit, &  
     433                          &            'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 
     434                     WRITE(numout,*) ' =========' 
     435                     WRITE(numout,*) 
     436                  ENDIF 
     437                  ! CW: warn, don't stop, as we may wish to violate this theoretical limit. 
     438                  ! CALL ctl_stop( 'ldf_dyn_init', 'Inconsistent Bi-Lap. GM diffusivity') 
     439               ELSE 
     440                  IF(lwp) THEN 
     441                     WRITE(numout,*) 
     442                     WRITE(numout,*) ' INFORMATION : Bi-Laplacian GM eddy viscosity is ',              & 
     443                          &            ' consistent with the model time step and grid setup: ', & 
     444                          &            ' LHS=',lhscrit, &  
     445                          &            'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 
     446                     WRITE(numout,*) ' =========' 
     447                     WRITE(numout,*) 
     448                  ENDIF 
     449               ENDIF 
     450               !----------------------------- 
     451            CASE(  12  )      !==  fixed profile: constant in interior, zero at bottom and linearly-tapering to zero at top, over depth shallower than rn_bgmzcrit  ==! 
     452               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity = constant in interior, zero at bottom and' 
     453               IF(lwp) WRITE(numout,*) '           linearly-tapering to zero at top, over depth shallower than ', rn_bgmzcrit, ' metres.' 
     454               IF(lwp) WRITE(numout,*) '           Interior viscous coef. = constant = ', bhm0, cl_Units 
     455 
     456               !     Impose linear taper from interior value to zero at zero depth, for depths < rn_bgmzcrit 
     457               !     N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 
     458 
     459               DO jk = 1, jpk 
     460 
     461                  IF (gdept_1d(jk) .lt. rn_bgmzcrit) THEN 
     462 
     463                     DO jj = 2, jpjm1 
     464                        DO ji = 2, jpim1 
     465                           bhmu(ji,jj,jk) = bhm0 * gdept_1d(jk)/rn_bgmzcrit 
     466                           bhmv(ji,jj,jk) = bhm0 * gdept_1d(jk)/rn_bgmzcrit 
     467                        ENDDO 
     468                     ENDDO 
     469 
     470                  ELSE 
     471 
     472                     ! constant at depth rn_bgmzcrit and larger 
     473                     DO jj = 2, jpjm1 
     474                        DO ji = 2, jpim1 
     475                           bhmu(ji,jj,jk) = bhm0 
     476                           bhmv(ji,jj,jk) = bhm0 
     477                        ENDDO 
     478                     ENDDO 
     479 
     480                  ENDIF 
     481 
     482               ENDDO 
     483 
     484 
     485               !     Surface BC : set to zero 
     486               bhmu(:,:,1)   = 0._wp   ;   bhmv(:,:,1)   = 0._wp 
     487               !     Flat bottom BC : set to zero 
     488               bhmu(:,:,jpk) = 0._wp   ;   bhmv(:,:,jpk) = 0._wp 
     489               !     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 
     490               !     and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp 
     491               DO jj = 2, jpjm1 
     492                  DO ji = 2, jpim1 
     493                     iku = mbku(ji,jj)+1 
     494                     ikv = mbkv(ji,jj)+1 
     495                     bhmu(:,:,iku) = 0._wp 
     496                     bhmv(:,:,ikv) = 0._wp 
     497                  ENDDO 
     498               ENDDO 
     499 
     500               !-------------------------------- 
     501            CASE(  13  )      !==  fixed profile: steady profile of the form bhm0 X dx^2 X dz^2 in interior, zero at bottom and top  ==! 
     502 
     503               cl_Units = ' 1/s' 
     504 
     505               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity : steady profile of the form' 
     506               IF(lwp) WRITE(numout,*) '           bhm0 X dx^2 X dz^2 in interior, zero boundary conds. at bottom and top.' 
     507               IF(lwp) WRITE(numout,*) '           In this case, bhm0 is the constant of proportionality,' 
     508               IF(lwp) WRITE(numout,*) '           dimensioned as [diffusivity]/L^4, i.e. bhm0=', bhm0, cl_Units 
     509 
     510               cl_Units = ' m4/s' 
     511 
     512               !     N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 
     513 
     514               DO jk = 1, jpk 
     515                  ! constant at depth rn_bgmzcrit and larger 
     516                  DO jj = 2, jpjm1 
     517                     DO ji = 2, jpim1 
     518                        bhmu(ji,jj,jk) = bhm0*e1u(ji,jj)**2*e3t_n(ji,jj,jk)**2 
     519                        bhmv(ji,jj,jk) = bhm0*e1v(ji,jj)**2*e3t_n(ji,jj,jk)**2 
     520                     ENDDO 
     521                  ENDDO 
     522               ENDDO 
     523 
     524 
     525               !     Surface BC : set to zero 
     526               bhmu(:,:,1)   = 0._wp   ;   bhmv(:,:,1)   = 0._wp 
     527               !     Flat bottom BC : set to zero 
     528               bhmu(:,:,jpk) = 0._wp   ;   bhmv(:,:,jpk) = 0._wp 
     529               !     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 
     530               !     and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp 
     531               DO jj = 2, jpjm1 
     532                  DO ji = 2, jpim1 
     533                     iku = mbku(ji,jj)+1 
     534                     ikv = mbkv(ji,jj)+1 
     535                     bhmu(:,:,iku) = 0._wp 
     536                     bhmv(:,:,ikv) = 0._wp 
     537                  ENDDO 
     538               ENDDO 
     539 
     540               !-------------------------------- 
     541 
     542            CASE DEFAULT 
     543               CALL ctl_stop('ldf_dyn_init: wrong choice for nn_bhm_ijk_t, the type of space-time variation of bhm') 
     544            END SELECT 
     545         ENDIF ! ln_dynldf_bgm 
    325546         ! 
    326547         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation  
     
    331552               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 
    332553               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 
     554               !CW 
     555            ELSEIF( ln_dynldf_bgm ) THEN                 !bi-Laplacian GM operator (mask only) 
     556               bhmu(:,:,1:jpkm1) =       bhmu(:,:,1:jpkm1)   * umask(:,:,1:jpkm1) 
     557               bhmv(:,:,1:jpkm1) =       bhmv(:,:,1:jpkm1)   * vmask(:,:,1:jpkm1) 
     558               ! 
    333559            ENDIF 
    334560         ENDIF 
     
    336562      ENDIF 
    337563      ! 
    338    END SUBROUTINE ldf_dyn_init 
    339  
    340  
    341    SUBROUTINE ldf_dyn( kt ) 
     564    END SUBROUTINE ldf_dyn_init 
     565 
     566 
     567    SUBROUTINE ldf_dyn( kt ) 
    342568      !!---------------------------------------------------------------------- 
    343569      !!                  ***  ROUTINE ldf_dyn  *** 
     
    366592      ! 
    367593      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==! 
    368       ! 
     594         ! 
    369595      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity ) 
    370596         ! 
     
    426652                  DO ji = 2, jpim1 
    427653                     zdb =   ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj)  & 
    428                         &  - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj) 
     654                          &  - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj) 
    429655                     dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk) 
    430656                  END DO 
     
    434660                  DO ji = 1, jpim1 
    435661                     zdb =   ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj)  & 
    436                         &  + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj) 
     662                          &  + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj) 
    437663                     dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk) 
    438664                  END DO 
     
    444670            ! 
    445671            DO jk = 1, jpkm1 
    446               ! 
     672               ! 
    447673               DO jj = 2, jpjm1                                ! T-point value 
    448674                  DO ji = fs_2, fs_jpim1 
     
    453679                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
    454680                     ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         & 
    455                         &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  & 
    456                         &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 
     681                          &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  & 
     682                          &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 
    457683                     ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
    458684                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt) 
     
    469695                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
    470696                     ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         & 
    471                         &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  & 
    472                         &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 
     697                          &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  & 
     698                          &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 
    473699                     ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
    474700                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt) 
Note: See TracChangeset for help on using the changeset viewer.