- Timestamp:
- 2020-03-11T10:26:18+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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 )
Note: See TracChangeset
for help on using the changeset viewer.
