Changeset 12532
- Timestamp:
- 2020-03-11T10:26:18+01:00 (5 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM
- Files:
-
- 5 edited
-
cfgs/SHARED/namelist_ref (modified) (2 diffs)
-
src/OCE/DIA/diaar5.F90 (modified) (2 diffs)
-
src/OCE/DYN/dynldf.F90 (modified) (4 diffs)
-
src/OCE/DYN/dynldf_lap_blp.F90 (modified) (4 diffs)
-
src/OCE/LDF/ldfdyn.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/cfgs/SHARED/namelist_ref
r12083 r12532 923 923 ln_dynldf_lap = .false. ! laplacian operator 924 924 ln_dynldf_blp = .false. ! bilaplacian operator 925 ln_dynldf_bgm = .false. ! bilaplician Gent-McWilliams 925 926 ! ! Direction of action : 926 927 ln_dynldf_lev = .false. ! iso-level … … 947 948 ! ! iso-neutral laplacian operator (ln_dynldf_iso=T) : 948 949 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 949 956 / 950 957 !----------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90
r11715 r12532 42 42 !!---------------------------------------------------------------------- 43 43 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 44 !! $Id $44 !! $Id: diaar5.F90 10425 2018-12-19 21:54:16Z smasson $ 45 45 !! Software governed by the CeCILL license (see ./LICENSE) 46 46 !!---------------------------------------------------------------------- … … 136 136 ! ! steric sea surface height 137 137 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 139 140 CALL iom_put( 'rhop', zrhop ) 140 141 ! -
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90
r11715 r12532 7 7 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 8 8 !! ! add velocity dependent coefficient and optional read in file 9 !! 4.0 ! 2020-02 (C. Wilson, ...) add bi-Laplacian GM implementation via momentum 9 10 !!---------------------------------------------------------------------- 10 11 … … 38 39 !!---------------------------------------------------------------------- 39 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 40 !! $Id $41 !! $Id: dynldf.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ 41 42 !! Software governed by the CeCILL license (see ./LICENSE) 42 43 !!---------------------------------------------------------------------- … … 66 67 CASE ( np_lap ) ; CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 67 68 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 68 70 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 69 71 ! … … 104 106 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> iso-level laplacian operator' 105 107 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' 106 109 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 107 110 END SELECT -
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90
r11715 r12532 5 5 !!====================================================================== 6 6 !! 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 7 8 !!---------------------------------------------------------------------- 8 9 … … 25 26 PUBLIC dyn_ldf_lap ! called by dynldf.F90 26 27 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 27 30 28 31 !! * Substitutions … … 30 33 !!---------------------------------------------------------------------- 31 34 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 32 !! $Id $35 !! $Id: dynldf_lap_blp.F90 10425 2018-12-19 21:54:16Z smasson $ 33 36 !! Software governed by the CeCILL license (see ./LICENSE) 34 37 !!---------------------------------------------------------------------- … … 104 107 END SUBROUTINE dyn_ldf_lap 105 108 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 106 273 107 274 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 8 8 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 9 9 !! ! 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 10 11 !!---------------------------------------------------------------------- 11 12 … … 15 16 !!---------------------------------------------------------------------- 16 17 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE dom_oce ! ocean space and time domain 18 19 USE phycst ! physical constants 19 20 USE ldfslp ! lateral diffusion: slopes of mixing orientation … … 35 36 LOGICAL , PUBLIC :: ln_dynldf_OFF !: No operator (i.e. no explicit diffusion) 36 37 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator 38 LOGICAL , PUBLIC :: ln_dynldf_bgm !: bi-Laplacian GM implementation via momentum 37 39 LOGICAL , PUBLIC :: ln_dynldf_blp !: bilaplacian operator 38 40 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction … … 42 44 ! ! time invariant coefficients: aht = 1/2 Ud*Ld (lap case) 43 45 ! ! 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. 44 47 REAL(wp), PUBLIC :: rn_Uv !: lateral viscous velocity [m/s] 45 48 REAL(wp), PUBLIC :: rn_Lv !: lateral viscous length [m] … … 50 53 ! ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 51 54 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] 53 57 ! !!* Parameter to control the type of lateral viscous operator 54 58 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator … … 57 61 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 58 62 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 60 64 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 61 65 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 62 66 63 67 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] 64 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 65 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) … … 76 81 !!---------------------------------------------------------------------- 77 82 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 78 !! $Id $83 !! $Id: ldfdyn.F90 11653 2019-10-04 12:34:02Z davestorkey $ 79 84 !! Software governed by the CeCILL license (see ./LICENSE) 80 85 !!---------------------------------------------------------------------- … … 91 96 !! ln_dynldf_lap = T laplacian operator 92 97 !! ln_dynldf_blp = T bilaplacian operator 98 !! ln_dynldf_bgm = T bi-Laplacian GM operator 93 99 !! - the parameter nn_ahm_ijk_t: 94 100 !! nn_ahm_ijk_t = 0 => = constant … … 103 109 !! ( L^2|D| laplacian operator 104 110 !! 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 105 116 !!---------------------------------------------------------------------- 106 117 INTEGER :: ji, jj, jk ! dummy loop indices 107 118 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 119 INTEGER :: iku, ikv, ikt ! local integer 108 120 REAL(wp) :: zah0, zah_max, zUfac ! local scalar 121 REAL(wp) :: bhm0, lhscrit, lhscritmax ! local scalar 109 122 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 110 123 !! 111 124 NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp, & ! type of operator 125 & ln_dynldf_bgm, & ! type of operator 112 126 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator 113 127 & 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 114 129 & rn_csmc , rn_minfac , rn_maxfac ! Smagorinsky settings 115 130 !!---------------------------------------------------------------------- … … 134 149 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 135 150 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 136 153 ! 137 154 WRITE(numout,*) ' direction of action :' … … 145 162 WRITE(numout,*) ' lateral viscous length (if cst) rn_Lv = ', rn_Lv, ' m' 146 163 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' 148 165 WRITE(numout,*) ' Smagorinsky settings (nn_ahm_ijk_t = 32) :' 149 166 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc … … 159 176 nldf_dyn = np_ERROR 160 177 ioptio = 0 178 !CW exclude combination of lap+blp (as in existing code), but allow other combinations 161 179 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 166 187 IF(.NOT.ln_dynldf_OFF ) THEN !== direction ==>> type of operator ==! 167 188 ioptio = 0 … … 209 230 ENDIF 210 231 ! 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 211 251 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 212 252 ! 253 IF( ierr == 3 ) CALL ctl_stop( 'rotated bi-Laplacian GM operator does not exist' ) 254 ! 255 213 256 IF( nldf_dyn == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 214 257 ! … … 222 265 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 223 266 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 267 CASE( np_bgm ) ; WRITE(numout,*) ' ==>>> iso-level bi-Laplacian GM operator' 224 268 END SELECT 225 269 WRITE(numout,*) … … 233 277 ! 234 278 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 ! 236 283 RETURN 237 284 ! … … 239 286 ! 240 287 ! ! 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 ! 242 290 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 243 291 ! 244 292 ahmt(:,:,:) = 0._wp ! init to 0 needed 245 293 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. 248 297 IF( ln_dynldf_lap ) THEN ; zUfac = r1_2 *rn_Uv ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 249 298 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 250 300 ENDIF 251 301 zah0 = zUfac * rn_Lv**inn ! mixing coefficient … … 265 315 ahmf(:,:,1) = zah0 266 316 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 267 ! 317 ! 268 318 CASE ( -20 ) !== fixed horizontal shape read in file ==! 269 319 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j) read in eddy_viscosity.nc file' … … 323 373 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 324 374 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 325 546 ! 326 547 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation … … 331 552 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 332 553 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 ! 333 559 ENDIF 334 560 ENDIF … … 336 562 ENDIF 337 563 ! 338 END SUBROUTINE ldf_dyn_init339 340 341 SUBROUTINE ldf_dyn( kt )564 END SUBROUTINE ldf_dyn_init 565 566 567 SUBROUTINE ldf_dyn( kt ) 342 568 !!---------------------------------------------------------------------- 343 569 !! *** ROUTINE ldf_dyn *** … … 366 592 ! 367 593 SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==! 368 !594 ! 369 595 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 370 596 ! … … 426 652 DO ji = 2, jpim1 427 653 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) 429 655 dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk) 430 656 END DO … … 434 660 DO ji = 1, jpim1 435 661 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) 437 663 dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk) 438 664 END DO … … 444 670 ! 445 671 DO jk = 1, jpkm1 446 !672 ! 447 673 DO jj = 2, jpjm1 ! T-point value 448 674 DO ji = fs_2, fs_jpim1 … … 453 679 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 454 680 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) ) ) 457 683 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 458 684 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) … … 469 695 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 470 696 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) ) ) 473 699 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 474 700 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.
