- Timestamp:
- 2020-06-02T17:46:26+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
r12667 r13005 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) … … 97 117 CONTAINS 98 118 99 SUBROUTINE dyn_vor( kt, K mm, puu, pvv, Krhs)119 SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs) 100 120 !!---------------------------------------------------------------------- 101 121 !! … … 107 127 !! for futher diagnostics (l_trddyn=T) 108 128 !!---------------------------------------------------------------------- 109 INTEGER , INTENT( in ) :: kt ! ocean time-step index 110 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 111 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 129 INTEGER :: ji, jj, jk ! dummy loop indice 130 INTEGER , INTENT( in ) :: kt ! ocean time-step index 131 INTEGER , INTENT( in ) :: Kmm, Krhs, Kbb ! ocean time level indices 132 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 112 133 ! 113 134 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuu, zvv 114 136 !!---------------------------------------------------------------------- 115 137 ! … … 165 187 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 166 188 CASE( np_ENE ) !* energy conserving scheme 189 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 190 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 191 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 192 193 CASE( np_ENE_adVO ) !* energy conserving scheme with alternating direction 194 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 195 196 !== Alternative direction - VOR only ==! 197 198 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 199 200 ALLOCATE( zuu(jpi,jpj,jpk) ) 201 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 202 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 203 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 204 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 205 DEALLOCATE( zuu ) 206 ELSE 207 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 208 209 ALLOCATE( zvv(jpi,jpj,jpk) ) 210 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 211 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 212 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 213 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 214 DEALLOCATE( zvv ) 215 ENDIF 216 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 217 CASE( np_ENE_adKE ) !* energy conserving scheme with alternating direction 218 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 219 220 !== Alternative direction - KEG only ==! 221 167 222 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 223 224 ALLOCATE( zuu(jpi,jpj,jpk) ) 225 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 226 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 227 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 228 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 229 DEALLOCATE( zuu ) 230 ELSE 231 232 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 233 234 ALLOCATE( zvv(jpi,jpj,jpk) ) 235 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 236 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 237 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 238 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 239 DEALLOCATE( zvv ) 240 ENDIF 168 241 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 242 243 CASE( np_ENE_adKEVO ) !* energy conserving scheme with alternating direction 244 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 245 246 !== Alternative direction - KE + VOR ==! 247 248 ALLOCATE( zuu(jpi,jpj,jpk) ) 249 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 250 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 251 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 252 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 253 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 254 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 255 DEALLOCATE( zuu ) 256 ELSE 257 258 ALLOCATE( zvv(jpi,jpj,jpk) ) 259 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 260 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 261 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 262 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 263 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 264 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 265 DEALLOCATE( zvv ) 266 ENDIF 169 267 CASE( np_ENS ) !* enstrophy conserving scheme 268 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 269 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 270 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 271 CASE( np_ENS_adVO ) !* enstrophy conserving scheme with alternating direction 272 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 273 274 !== Alternative direction - VOR only ==! 275 276 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 277 278 ALLOCATE( zuu(jpi,jpj,jpk) ) 279 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 280 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 281 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 282 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 283 DEALLOCATE( zuu ) 284 ELSE 285 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 286 287 ALLOCATE( zvv(jpi,jpj,jpk) ) 288 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 289 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 290 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 291 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 292 DEALLOCATE( zvv ) 293 ENDIF 294 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 295 CASE( np_ENS_adKE ) !* enstrophy conserving scheme with alternating direction 296 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 297 298 !== Alternative direction - KEG only ==! 299 170 300 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 301 302 ALLOCATE( zuu(jpi,jpj,jpk) ) 303 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 304 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 305 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 306 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 307 DEALLOCATE( zuu ) 308 ELSE 309 310 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 311 312 ALLOCATE( zvv(jpi,jpj,jpk) ) 313 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 314 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 315 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 316 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 317 DEALLOCATE( zvv ) 318 ENDIF 319 CASE( np_ENS_adKEVO ) !* enstrophy conserving scheme with alternating direction 320 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 321 322 !== Alternative direction - KE + VOR ==! 323 324 ALLOCATE( zuu(jpi,jpj,jpk) ) 325 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 326 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 327 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 328 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 329 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 330 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 331 DEALLOCATE( zuu ) 332 ELSE 333 334 ALLOCATE( zvv(jpi,jpj,jpk) ) 335 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 336 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 337 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 338 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 339 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 340 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 341 DEALLOCATE( zvv ) 342 ENDIF 171 343 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 172 344 CASE( np_MIX ) !* mixed ene-ens scheme … … 313 485 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 314 486 !!---------------------------------------------------------------------- 315 INTEGER , INTENT(in ) :: kt ! ocean time-step index316 INTEGER , INTENT(in ) :: Kmm ! ocean time level index317 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric318 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities319 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend487 INTEGER , INTENT(in ) :: kt ! ocean time-step index 488 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 489 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 490 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 491 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 320 492 ! 321 493 INTEGER :: ji, jj, jk ! dummy loop indices … … 371 543 END SELECT 372 544 ! 373 ! !== horizontal fluxes and potential vorticity ==! 374 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 375 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 376 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 377 ! 378 ! !== compute and add the vorticity term trend =! 379 DO_2D_00_00 380 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 381 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 382 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 383 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 384 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 ) 385 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 ) 386 END_2D 545 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 546 ! 547 ! !== horizontal fluxes and potential vorticity ==! 548 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 549 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 550 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 551 ! 552 ! !== compute and add the vorticity term trend =! 553 DO_2D_00_00 554 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 555 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 556 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 557 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 558 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 ) 559 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 ) 560 END_2D 561 ! 562 ! 563 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 564 ! 565 ! 566 ! !== horizontal fluxes and potential vorticity ==! 567 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 568 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 569 ! 570 ! !== compute and add the vorticity term trend =! 571 DO_2D_00_00 572 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 573 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 574 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 ) 575 END_2D 576 ! 577 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 578 ! 579 ! 580 ! !== horizontal fluxes and potential vorticity ==! 581 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 582 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 583 ! 584 ! !== compute and add the vorticity term trend =! 585 DO_2D_00_00 586 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 587 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 588 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 ) 589 END_2D 590 ! 591 ENDIF 387 592 ! ! =============== 388 593 END DO ! End of slab … … 411 616 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 412 617 !!---------------------------------------------------------------------- 413 INTEGER , INTENT(in ) :: kt ! ocean time-step index414 INTEGER , INTENT(in ) :: Kmm ! ocean time level index415 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric416 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities417 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend618 INTEGER , INTENT(in ) :: kt ! ocean time-step index 619 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 620 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 621 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 622 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 418 623 ! 419 624 INTEGER :: ji, jj, jk ! dummy loop indices … … 469 674 ! 470 675 ! 471 ! !== horizontal fluxes and potential vorticity ==! 472 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 473 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 474 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 475 ! 476 ! !== compute and add the vorticity term trend =! 477 DO_2D_00_00 478 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 479 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 480 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 481 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 482 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 483 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 484 END_2D 676 !!an wut ? v et u 677 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 678 ! 679 ! !== horizontal fluxes and potential vorticity ==! 680 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 681 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 682 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 683 ! 684 ! !== compute and add the vorticity term trend =! 685 DO_2D_00_00 686 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 687 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 688 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 689 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 690 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 691 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 692 END_2D 693 ! 694 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 695 ! 696 ! 697 ! !== horizontal fluxes and potential vorticity ==! 698 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 699 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 700 ! 701 ! !== compute and add the vorticity term trend =! 702 DO_2D_00_00 703 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 704 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 705 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 706 END_2D 707 ! 708 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 709 ! 710 ! 711 ! !== horizontal fluxes and potential vorticity ==! 712 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 713 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 714 ! 715 ! !== compute and add the vorticity term trend =! 716 DO_2D_00_00 717 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 718 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 719 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 720 END_2D 721 ! 722 ENDIF 485 723 ! ! =============== 486 724 END DO ! End of slab … … 762 1000 !! 763 1001 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 764 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 1002 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk, & 1003 & ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & ! Alternative direction parameters 1004 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 765 1005 !!---------------------------------------------------------------------- 766 1006 ! … … 779 1019 IF(lwp) THEN ! Namelist print 780 1020 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 781 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens782 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene783 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT784 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT785 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een786 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f787 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix788 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk1021 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 1022 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 1023 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT 1024 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT 1025 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 1026 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 1027 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 1028 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 789 1029 ENDIF 790 1030 … … 799 1039 DO_3D_10_10( 1, jpk ) 800 1040 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 801 & + tmask(ji,jj ,jk) + tmask(ji+1,jj +1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp1041 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 802 1042 END_3D 803 1043 ! … … 809 1049 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 810 1050 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 1051 IF( ln_dynvor_ens_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adVO ; ENDIF 1052 IF( ln_dynvor_ens_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKE ; ENDIF 1053 IF( ln_dynvor_ens_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKEVO ; ENDIF 811 1054 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 1055 IF( ln_dynvor_ene_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adVO ; ENDIF 1056 IF( ln_dynvor_ene_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKE ; ENDIF 1057 IF( ln_dynvor_ene_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKEVO ; ENDIF 812 1058 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 813 1059 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF … … 826 1072 CASE( np_VEC_c2 ) 827 1073 IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity' 828 nrvm = np_RVO ! relative vorticity 1074 nrvm = np_RVO ! relative vorticity 829 1075 ntot = np_CRV ! relative + planetary vorticity 830 1076 CASE( np_FLX_c2 , np_FLX_ubs ) … … 856 1102 WRITE(numout,*) 857 1103 SELECT CASE( nvor_scheme ) 858 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 859 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 860 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 861 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 862 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 863 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 1104 1105 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 1106 CASE( np_ENS_adVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 1107 CASE( np_ENS_adKE ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 1108 CASE( np_ENS_adKEVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 1109 1110 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 1111 CASE( np_ENE_adVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 1112 CASE( np_ENE_adKE ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 1113 CASE( np_ENE_adKEVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 1114 1115 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 1116 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 1117 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 1118 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 864 1119 END SELECT 865 1120 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
