- Timestamp:
- 2020-06-24T15:31:32+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
r13151 r13154 22 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 23 23 !! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity 24 !! 4.x ! 2020-03 (G. Madec, A. Nasser) alternate direction computation of vorticity tendancy 25 !! ! for ENS, ENE 24 26 !!---------------------------------------------------------------------- 25 27 … … 45 47 USE lib_mpp ! MPP library 46 48 USE timing ! Timing 49 !!anAD only 50 USE dynkeg, ONLY : dyn_kegAD 51 USE dynadv, ONLY : nn_dynkeg 47 52 48 53 IMPLICIT NONE … … 53 58 54 59 ! !!* Namelist namdyn_vor: vorticity term 55 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 56 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 57 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 58 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 59 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 60 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 61 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 62 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 60 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 61 LOGICAL, PUBLIC :: ln_dynvor_ens_adVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 62 LOGICAL, PUBLIC :: ln_dynvor_ens_adKE = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 63 LOGICAL, PUBLIC :: ln_dynvor_ens_adKEVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 64 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 65 LOGICAL, PUBLIC :: ln_dynvor_ene_adVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 66 LOGICAL, PUBLIC :: ln_dynvor_ene_adKE = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 67 LOGICAL, PUBLIC :: ln_dynvor_ene_adKEVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 68 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 69 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 70 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 71 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 72 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 73 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 63 74 64 75 INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme … … 70 81 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 71 82 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 72 73 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 83 !!an 84 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKE = 11 ! ENS scheme - AD scheme (KE only) 85 INTEGER, PUBLIC, PARAMETER :: np_ENS_adVO = 12 ! ENS scheme - AD scheme (VOR only) 86 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKEVO = 13 ! ENS scheme - AD scheme (KE+VOR) 87 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKE = 21 ! ENE scheme - AD scheme (KE only) 88 INTEGER, PUBLIC, PARAMETER :: np_ENE_adVO = 22 ! ENE scheme - AD scheme (VOR only) 89 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKEVO = 23 ! ENE scheme - AD scheme (KE+VOR) 90 !!an 91 92 !!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO 93 INTEGER, PUBLIC :: ncor, nrvm, ntot ! choice of calculated vorticity 74 94 ! ! associated indices: 75 95 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) … … 99 119 CONTAINS 100 120 101 SUBROUTINE dyn_vor( kt, K mm, puu, pvv, Krhs)121 SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs) 102 122 !!---------------------------------------------------------------------- 103 123 !! … … 109 129 !! for futher diagnostics (l_trddyn=T) 110 130 !!---------------------------------------------------------------------- 111 INTEGER , INTENT( in ) :: kt ! ocean time-step index 112 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 113 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 131 INTEGER :: ji, jj, jk ! dummy loop indice 132 INTEGER , INTENT( in ) :: kt ! ocean time-step index 133 INTEGER , INTENT( in ) :: Kmm, Krhs, Kbb ! ocean time level indices 134 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 114 135 ! 115 136 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 137 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuu, zvv 116 138 !!---------------------------------------------------------------------- 117 139 ! … … 167 189 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 168 190 CASE( np_ENE ) !* energy conserving scheme 191 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 192 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 193 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 194 195 CASE( np_ENE_adVO ) !* energy conserving scheme with alternating direction 196 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 197 198 !== Alternative direction - VOR only ==! 199 200 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 201 202 ALLOCATE( zuu(jpi,jpj,jpk) ) 203 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 204 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 205 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 206 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 207 DEALLOCATE( zuu ) 208 ELSE 209 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 210 211 ALLOCATE( zvv(jpi,jpj,jpk) ) 212 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 213 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 214 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 215 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 216 DEALLOCATE( zvv ) 217 ENDIF 218 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 219 CASE( np_ENE_adKE ) !* energy conserving scheme with alternating direction 220 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 221 222 !== Alternative direction - KEG only ==! 223 169 224 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 225 226 ALLOCATE( zuu(jpi,jpj,jpk) ) 227 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 228 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 229 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 230 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 231 DEALLOCATE( zuu ) 232 ELSE 233 234 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 235 236 ALLOCATE( zvv(jpi,jpj,jpk) ) 237 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 238 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 239 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 240 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 241 DEALLOCATE( zvv ) 242 ENDIF 170 243 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 244 245 CASE( np_ENE_adKEVO ) !* energy conserving scheme with alternating direction 246 ! equivalent ADE pas sur coriolis 247 ! appel vor_ene(ncor) 248 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 249 250 !== Alternative direction - KE + VOR ==! 251 252 ALLOCATE( zuu(jpi,jpj,jpk) ) 253 !puis call que sur nrvm 254 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 255 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 256 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 257 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 258 !puis call que sur nrvm 259 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 260 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 261 DEALLOCATE( zuu ) 262 ELSE 263 264 ALLOCATE( zvv(jpi,jpj,jpk) ) 265 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 266 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 267 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 268 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 269 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 270 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 271 DEALLOCATE( zvv ) 272 ENDIF 171 273 CASE( np_ENS ) !* enstrophy conserving scheme 274 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 275 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 276 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 277 CASE( np_ENS_adVO ) !* enstrophy conserving scheme with alternating direction 278 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 279 280 !== Alternative direction - VOR only ==! 281 282 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 283 284 ALLOCATE( zuu(jpi,jpj,jpk) ) 285 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 286 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 287 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 288 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 289 DEALLOCATE( zuu ) 290 ELSE 291 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 292 293 ALLOCATE( zvv(jpi,jpj,jpk) ) 294 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 295 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 296 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 297 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 298 DEALLOCATE( zvv ) 299 ENDIF 300 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 301 CASE( np_ENS_adKE ) !* enstrophy conserving scheme with alternating direction 302 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 303 304 !== Alternative direction - KEG only ==! 305 172 306 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 307 308 ALLOCATE( zuu(jpi,jpj,jpk) ) 309 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 310 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 311 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 312 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 313 DEALLOCATE( zuu ) 314 ELSE 315 316 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 317 318 ALLOCATE( zvv(jpi,jpj,jpk) ) 319 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 320 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 321 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 322 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 323 DEALLOCATE( zvv ) 324 ENDIF 325 CASE( np_ENS_adKEVO ) !* enstrophy conserving scheme with alternating direction 326 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 327 328 !== Alternative direction - KE + VOR ==! 329 330 ALLOCATE( zuu(jpi,jpj,jpk) ) 331 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 332 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 333 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 334 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 335 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 336 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 337 DEALLOCATE( zuu ) 338 ELSE 339 340 ALLOCATE( zvv(jpi,jpj,jpk) ) 341 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 342 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 343 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 344 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 345 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 346 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 347 DEALLOCATE( zvv ) 348 ENDIF 173 349 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 174 350 CASE( np_MIX ) !* mixed ene-ens scheme … … 319 495 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 320 496 !!---------------------------------------------------------------------- 321 INTEGER , INTENT(in ) :: kt ! ocean time-step index322 INTEGER , INTENT(in ) :: Kmm ! ocean time level index323 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric324 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities325 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend497 INTEGER , INTENT(in ) :: kt ! ocean time-step index 498 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 499 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 500 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 501 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 326 502 ! 327 503 INTEGER :: ji, jj, jk ! dummy loop indices 328 504 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 329 505 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace 506 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! unpenalised scale factors 330 507 !!---------------------------------------------------------------------- 331 508 ! … … 377 554 END SELECT 378 555 ! 379 ! !== horizontal fluxes and potential vorticity ==! 380 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 381 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 382 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 383 ! 384 ! !== compute and add the vorticity term trend =! 385 DO_2D_00_00 386 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 387 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 388 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 389 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 390 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 391 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 392 END_2D 556 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 557 ! 558 ! !== horizontal fluxes and potential vorticity ==! 559 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 560 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 561 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 562 ! 563 ! !== compute and add the vorticity term trend =! 564 DO_2D_00_00 565 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 566 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 567 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 568 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 569 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 570 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 571 END_2D 572 ! 573 ! 574 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 575 ! 576 ! 577 ! !== horizontal fluxes and potential vorticity ==! 578 !!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée 579 ! on transportera (en AD) de la bonne manière 580 !zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 581 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm) 582 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 583 ! 584 ! !== compute and add the vorticity term trend =! 585 DO_2D_00_00 586 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 587 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 588 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 589 END_2D 590 ! 591 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 592 ! 593 ! 594 ! !== horizontal fluxes and potential vorticity ==! 595 !!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée 596 !zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 597 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm) 598 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 599 ! 600 ! !== compute and add the vorticity term trend =! 601 DO_2D_00_00 602 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 603 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 604 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 605 END_2D 606 ! 607 ENDIF 393 608 ! ! =============== 394 609 END DO ! End of slab … … 417 632 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 418 633 !!---------------------------------------------------------------------- 419 INTEGER , INTENT(in ) :: kt ! ocean time-step index420 INTEGER , INTENT(in ) :: Kmm ! ocean time level index421 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric422 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities423 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend634 INTEGER , INTENT(in ) :: kt ! ocean time-step index 635 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 636 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 637 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 638 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 424 639 ! 425 640 INTEGER :: ji, jj, jk ! dummy loop indices … … 475 690 ! 476 691 ! 477 ! !== horizontal fluxes and potential vorticity ==! 478 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 479 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 480 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 481 ! 482 ! !== compute and add the vorticity term trend =! 483 DO_2D_00_00 484 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 485 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 486 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 487 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 488 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 489 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 490 END_2D 692 !!an wut ? v et u 693 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 694 ! 695 ! !== horizontal fluxes and potential vorticity ==! 696 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 697 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 698 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 699 ! 700 ! !== compute and add the vorticity term trend =! 701 DO_2D_00_00 702 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 704 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 706 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 708 END_2D 709 ! 710 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 711 ! 712 ! 713 ! !== horizontal fluxes and potential vorticity ==! 714 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 715 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 716 ! 717 ! !== compute and add the vorticity term trend =! 718 DO_2D_00_00 719 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 720 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 721 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 722 END_2D 723 ! 724 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 725 ! 726 ! 727 ! !== horizontal fluxes and potential vorticity ==! 728 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 729 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 730 ! 731 ! !== compute and add the vorticity term trend =! 732 DO_2D_00_00 733 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 734 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 735 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 736 END_2D 737 ! 738 ENDIF 491 739 ! ! =============== 492 740 END DO ! End of slab … … 772 1020 !! 773 1021 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 774 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 1022 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk, & 1023 & ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & ! Alternative direction parameters 1024 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 775 1025 !!---------------------------------------------------------------------- 776 1026 ! … … 809 1059 DO_3D_10_10( 1, jpk ) 810 1060 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 811 & + tmask(ji,jj ,jk) + tmask(ji+1,jj +1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp1061 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 812 1062 END_3D 813 1063 ! … … 819 1069 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 820 1070 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 1071 IF( ln_dynvor_ens_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adVO ; ENDIF 1072 IF( ln_dynvor_ens_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKE ; ENDIF 1073 IF( ln_dynvor_ens_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKEVO ; ENDIF 821 1074 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 1075 IF( ln_dynvor_ene_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adVO ; ENDIF 1076 IF( ln_dynvor_ene_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKE ; ENDIF 1077 IF( ln_dynvor_ene_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKEVO ; ENDIF 822 1078 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 823 1079 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF … … 866 1122 WRITE(numout,*) 867 1123 SELECT CASE( nvor_scheme ) 868 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 869 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 870 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 871 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 872 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 873 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 1124 1125 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 1126 CASE( np_ENS_adVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 1127 CASE( np_ENS_adKE ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 1128 CASE( np_ENS_adKEVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 1129 1130 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 1131 CASE( np_ENE_adVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 1132 CASE( np_ENE_adKE ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 1133 CASE( np_ENE_adKEVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 1134 1135 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 1136 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 1137 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 1138 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 874 1139 END SELECT 875 1140 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
