- Timestamp:
- 2018-04-30T15:39:18+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r9190 r9528 6 6 !!====================================================================== 7 7 !! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code 8 !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code8 !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code 9 9 !! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays 10 10 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module … … 19 19 !! - ! 2016-12 (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 20 20 !! 4.0 ! 2017-07 (G. Madec) linear dynamics + trends diag. with Stokes-Coriolis 21 !! - ! 2018-03 (G. Madec) add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 21 23 !!---------------------------------------------------------------------- 22 24 … … 50 52 51 53 ! !!* Namelist namdyn_vor: vorticity term 52 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE) 53 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 54 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 55 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN) 54 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 55 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 56 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 57 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 58 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 56 59 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 60 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 57 61 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 58 62 59 INTEGER :: nvor_scheme ! choice of the type of advection scheme 60 ! ! associated indices: 63 INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme 64 ! ! associated indices: 65 INTEGER, PUBLIC, PARAMETER :: np_ENS = 0 ! ENS scheme 61 66 INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme 62 INTEGER, PUBLIC, PARAMETER :: np_EN S = 2 ! ENS scheme63 INTEGER, PUBLIC, PARAMETER :: np_ MIX = 3 ! MIX scheme67 INTEGER, PUBLIC, PARAMETER :: np_ENT = 2 ! ENT scheme (t-point vorticity) 68 INTEGER, PUBLIC, PARAMETER :: np_EET = 3 ! EET scheme (EEN using e3t) 64 69 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 70 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 65 71 66 72 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 67 73 ! ! associated indices: 68 INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 69 INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity 70 INTEGER, PARAMETER :: np_MET = 3 ! metric term 71 INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 72 INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term 74 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 75 INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity 76 INTEGER, PUBLIC, PARAMETER :: np_MET = 3 ! metric term 77 INTEGER, PUBLIC, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 78 INTEGER, PUBLIC, PARAMETER :: np_CME = 5 ! Coriolis + metric term 79 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2u_2 ! = di(e2u)/2 used in T-point metric term calculation 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - - 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2u)/(2*e1e2f) used in F-point metric term calculation 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1v)/(2*e1e2f) - - - - 73 84 74 85 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 … … 109 120 ztrdv(:,:,:) = va(:,:,:) 110 121 SELECT CASE( nvor_scheme ) 122 CASE( np_ENS ) ; CALL vor_ens( kt, ncor, un , vn , ua, va ) ! enstrophy conserving scheme 123 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 111 124 CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, ncor, un , vn , ua, va ) ! energy conserving scheme 112 125 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 113 CASE( np_ENS ) ; CALL vor_ens( kt, ncor, un , vn , ua, va ) ! enstrophy conserving scheme 114 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 126 CASE( np_ENT ) ; CALL vor_enT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (T-pts) 127 IF( ln_stcor ) CALL vor_enT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 128 CASE( np_EET ) ; CALL vor_eeT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 129 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 115 130 CASE( np_EEN ) ; CALL vor_een( kt, ncor, un , vn , ua, va ) ! energy & enstrophy scheme 116 131 IF( ln_stcor ) CALL vor_een( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend … … 124 139 ztrdv(:,:,:) = va(:,:,:) 125 140 SELECT CASE( nvor_scheme ) 141 CASE( np_ENT ) ; CALL vor_enT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (T-pts) 142 CASE( np_EET ) ; CALL vor_eeT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 126 143 CASE( np_ENE ) ; CALL vor_ene( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme 127 144 CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! enstrophy conserving scheme … … 138 155 ! 139 156 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 157 CASE( np_ENT ) !* energy conserving scheme (T-pts) 158 CALL vor_enT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 159 IF( ln_stcor ) CALL vor_enT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 160 CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) 161 CALL vor_eeT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 162 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 140 163 CASE( np_ENE ) !* energy conserving scheme 141 164 CALL vor_ene( kt, ntot, un , vn , ua, va ) ! total vorticity trend … … 164 187 165 188 189 SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs ) 190 !!---------------------------------------------------------------------- 191 !! *** ROUTINE vor_enT *** 192 !! 193 !! ** Purpose : Compute the now total vorticity trend and add it to 194 !! the general trend of the momentum equation. 195 !! 196 !! ** Method : Trend evaluated using now fields (centered in time) 197 !! and t-point evaluation of vorticity (planetary and relative). 198 !! conserves the horizontal kinetic energy. 199 !! The general trend of momentum is increased due to the vorticity 200 !! term which is given by: 201 !! voru = 1/bu mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t mj[vn] ] 202 !! vorv = 1/bv mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f mj[un] ] 203 !! where rvor is the relative vorticity at f-point 204 !! 205 !! ** Action : - Update (ua,va) with the now vorticity term trend 206 !!---------------------------------------------------------------------- 207 INTEGER , INTENT(in ) :: kt ! ocean time-step index 208 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 209 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities 210 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 211 ! 212 INTEGER :: ji, jj, jk ! dummy loop indices 213 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 214 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zwt ! 2D workspace 215 !!---------------------------------------------------------------------- 216 ! 217 IF( kt == nit000 ) THEN 218 IF(lwp) WRITE(numout,*) 219 IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 220 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 221 ENDIF 222 ! 223 ! ! =============== 224 DO jk = 1, jpkm1 ! Horizontal slab 225 ! ! =============== 226 ! 227 SELECT CASE( kvor ) !== volume weighted vorticity considered ==! 228 CASE ( np_COR ) !* Coriolis (planetary vorticity) 229 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 230 CASE ( np_RVO ) !* relative vorticity 231 DO jj = 1, jpjm1 232 DO ji = 1, jpim1 233 zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 234 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 235 END DO 236 END DO 237 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 238 DO jj = 1, jpjm1 239 DO ji = 1, jpim1 240 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 241 END DO 242 END DO 243 ENDIF 244 CALL lbc_lnk( zwz, 'F', 1. ) 245 DO jj = 2, jpj 246 DO ji = 2, jpi ! vector opt. 247 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) & 248 & + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 249 END DO 250 END DO 251 CASE ( np_MET ) !* metric term 252 DO jj = 2, jpj 253 DO ji = 2, jpi 254 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 255 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t_n(ji,jj,jk) 256 END DO 257 END DO 258 CASE ( np_CRV ) !* Coriolis + relative vorticity 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 ! relative vorticity 261 zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 262 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 263 END DO 264 END DO 265 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 266 DO jj = 1, jpjm1 267 DO ji = 1, jpim1 268 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 269 END DO 270 END DO 271 ENDIF 272 CALL lbc_lnk( zwz, 'F', 1. ) 273 DO jj = 2, jpj 274 DO ji = 2, jpi ! vector opt. 275 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) & 276 & + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 277 END DO 278 END DO 279 CASE ( np_CME ) !* Coriolis + metric 280 DO jj = 2, jpj 281 DO ji = 2, jpi ! vector opt. 282 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 283 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 284 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t_n(ji,jj,jk) 285 END DO 286 END DO 287 CASE DEFAULT ! error 288 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 289 END SELECT 290 ! 291 ! !== compute and add the vorticity term trend =! 292 DO jj = 2, jpjm1 293 DO ji = 2, jpim1 ! vector opt. 294 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) & 295 & * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) & 296 & + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) ) 297 ! 298 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) & 299 & * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) & 300 & + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) ) 301 END DO 302 END DO 303 ! ! =============== 304 END DO ! End of slab 305 ! ! =============== 306 END SUBROUTINE vor_enT 307 308 166 309 SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 167 310 !!---------------------------------------------------------------------- … … 217 360 DO jj = 1, jpjm1 218 361 DO ji = 1, fs_jpim1 ! vector opt. 219 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 220 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 221 & * 0.5 * r1_e1e2f(ji,jj) 362 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 363 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 222 364 END DO 223 365 END DO … … 225 367 DO jj = 1, jpjm1 226 368 DO ji = 1, fs_jpim1 ! vector opt. 227 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 228 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 229 & * r1_e1e2f(ji,jj) 369 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 370 & - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 230 371 END DO 231 372 END DO … … 233 374 DO jj = 1, jpjm1 234 375 DO ji = 1, fs_jpim1 ! vector opt. 235 zwz(ji,jj) = ff_f(ji,jj) & 236 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 237 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 238 & * 0.5 * r1_e1e2f(ji,jj) 376 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 377 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 239 378 END DO 240 379 END DO … … 328 467 DO jj = 1, jpjm1 329 468 DO ji = 1, fs_jpim1 ! vector opt. 330 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 331 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 332 & * 0.5 * r1_e1e2f(ji,jj) 469 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 470 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 333 471 END DO 334 472 END DO … … 336 474 DO jj = 1, jpjm1 337 475 DO ji = 1, fs_jpim1 ! vector opt. 338 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 339 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 340 & * r1_e1e2f(ji,jj) 476 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 477 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 341 478 END DO 342 479 END DO … … 344 481 DO jj = 1, jpjm1 345 482 DO ji = 1, fs_jpim1 ! vector opt. 346 zwz(ji,jj) = ff_f(ji,jj) & 347 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 348 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 349 & * 0.5 * r1_e1e2f(ji,jj) 483 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 484 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 350 485 END DO 351 486 END DO … … 412 547 INTEGER :: ierr ! local integer 413 548 REAL(wp) :: zua, zva ! local scalars 414 REAL(wp) :: zmsk, ze3 415 REAL(wp), DIMENSION(jpi,jpj) ::zwx , zwy , zwz , z1_e3f416 REAL(wp), DIMENSION(jpi,jpj) ::ztnw, ztne, ztsw, ztse549 REAL(wp) :: zmsk, ze3f ! local scalars 550 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , zwz , z1_e3f 551 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 417 552 !!---------------------------------------------------------------------- 418 553 ! … … 431 566 DO jj = 1, jpjm1 432 567 DO ji = 1, fs_jpim1 ! vector opt. 433 ze3 568 ze3f = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 434 569 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 435 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3436 ELSE ; z1_e3f(ji,jj) = 0._wp570 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 571 ELSE ; z1_e3f(ji,jj) = 0._wp 437 572 ENDIF 438 573 END DO … … 441 576 DO jj = 1, jpjm1 442 577 DO ji = 1, fs_jpim1 ! vector opt. 443 ze3 578 ze3f = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 444 579 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 445 580 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 446 581 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 447 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3448 ELSE ; z1_e3f(ji,jj) = 0._wp582 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f 583 ELSE ; z1_e3f(ji,jj) = 0._wp 449 584 ENDIF 450 585 END DO … … 462 597 DO jj = 1, jpjm1 463 598 DO ji = 1, fs_jpim1 ! vector opt. 464 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 465 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 466 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 599 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 600 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 467 601 END DO 468 602 END DO … … 470 604 DO jj = 1, jpjm1 471 605 DO ji = 1, fs_jpim1 ! vector opt. 472 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 473 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 474 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 606 zwz(ji,jj) = ( ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 607 & - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 475 608 END DO 476 609 END DO … … 478 611 DO jj = 1, jpjm1 479 612 DO ji = 1, fs_jpim1 ! vector opt. 480 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)&481 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) )&482 & * r1_e1e2f(ji,jj)) * z1_e3f(ji,jj)613 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 614 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 615 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 483 616 END DO 484 617 END DO … … 486 619 DO jj = 1, jpjm1 487 620 DO ji = 1, fs_jpim1 ! vector opt. 488 zwz(ji,jj) = ( ff_f(ji,jj) & 489 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 490 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 491 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 621 zwz(ji,jj) = ( ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 622 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 492 623 END DO 493 624 END DO … … 543 674 544 675 676 677 SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva ) 678 !!---------------------------------------------------------------------- 679 !! *** ROUTINE vor_eeT *** 680 !! 681 !! ** Purpose : Compute the now total vorticity trend and add it to 682 !! the general trend of the momentum equation. 683 !! 684 !! ** Method : Trend evaluated using now fields (centered in time) 685 !! and the Arakawa and Lamb (1980) vector form formulation using 686 !! a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 687 !! The change consists in 688 !! Add this trend to the general momentum trend (ua,va). 689 !! 690 !! ** Action : - Update (ua,va) with the now vorticity term trend 691 !! 692 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 693 !!---------------------------------------------------------------------- 694 INTEGER , INTENT(in ) :: kt ! ocean time-step index 695 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 696 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 697 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 698 ! 699 INTEGER :: ji, jj, jk ! dummy loop indices 700 INTEGER :: ierr ! local integer 701 REAL(wp) :: zua, zva ! local scalars 702 REAL(wp) :: zmsk, z1_e3t ! local scalars 703 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , zwz 704 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 705 !!---------------------------------------------------------------------- 706 ! 707 IF( kt == nit000 ) THEN 708 IF(lwp) WRITE(numout,*) 709 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 710 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 711 ENDIF 712 ! 713 ! ! =============== 714 DO jk = 1, jpkm1 ! Horizontal slab 715 ! ! =============== 716 ! 717 ! 718 SELECT CASE( kvor ) !== vorticity considered ==! 719 CASE ( np_COR ) !* Coriolis (planetary vorticity) 720 DO jj = 1, jpjm1 721 DO ji = 1, fs_jpim1 ! vector opt. 722 zwz(ji,jj) = ff_f(ji,jj) 723 END DO 724 END DO 725 CASE ( np_RVO ) !* relative vorticity 726 DO jj = 1, jpjm1 727 DO ji = 1, fs_jpim1 ! vector opt. 728 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 729 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 730 & * r1_e1e2f(ji,jj) 731 END DO 732 END DO 733 CASE ( np_MET ) !* metric term 734 DO jj = 1, jpjm1 735 DO ji = 1, fs_jpim1 ! vector opt. 736 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 737 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 738 END DO 739 END DO 740 CASE ( np_CRV ) !* Coriolis + relative vorticity 741 DO jj = 1, jpjm1 742 DO ji = 1, fs_jpim1 ! vector opt. 743 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 744 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 745 & * r1_e1e2f(ji,jj) ) 746 END DO 747 END DO 748 CASE ( np_CME ) !* Coriolis + metric 749 DO jj = 1, jpjm1 750 DO ji = 1, fs_jpim1 ! vector opt. 751 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 752 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 753 END DO 754 END DO 755 CASE DEFAULT ! error 756 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 757 END SELECT 758 ! 759 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 760 DO jj = 1, jpjm1 761 DO ji = 1, fs_jpim1 ! vector opt. 762 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 763 END DO 764 END DO 765 ENDIF 766 ! 767 CALL lbc_lnk( zwz, 'F', 1. ) 768 ! 769 ! !== horizontal fluxes ==! 770 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 771 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 772 773 ! !== compute and add the vorticity term trend =! 774 jj = 2 775 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 776 DO ji = 2, jpi ! split in 2 parts due to vector opt. 777 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 778 ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t 779 ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t 780 ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t 781 ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t 782 END DO 783 DO jj = 3, jpj 784 DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3 785 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 786 ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t 787 ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t 788 ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t 789 ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t 790 END DO 791 END DO 792 DO jj = 2, jpjm1 793 DO ji = fs_2, fs_jpim1 ! vector opt. 794 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 795 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 796 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 797 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 798 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 799 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 800 END DO 801 END DO 802 ! ! =============== 803 END DO ! End of slab 804 ! ! =============== 805 END SUBROUTINE vor_eeT 806 807 545 808 SUBROUTINE dyn_vor_init 546 809 !!--------------------------------------------------------------------- … … 550 813 !! tracer advection schemes 551 814 !!---------------------------------------------------------------------- 552 INTEGER :: ioptio ! local integer 553 INTEGER :: ji, jj, jk ! dummy loop indices 554 INTEGER :: ios ! Local integer output status for namelist read 555 !! 556 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, & 557 & ln_dynvor_een, nn_een_e3f , ln_dynvor_msk 558 !!---------------------------------------------------------------------- 559 815 INTEGER :: ji, jj, jk ! dummy loop indices 816 INTEGER :: ioptio, ios ! local integer 817 !! 818 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 819 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 820 !!---------------------------------------------------------------------- 821 ! 822 IF(lwp) THEN 823 WRITE(numout,*) 824 WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' 825 WRITE(numout,*) '~~~~~~~~~~~~' 826 ENDIF 827 ! 560 828 REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options 561 829 READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901) … … 565 833 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp ) 566 834 IF(lwm) WRITE ( numond, namdyn_vor ) 567 835 ! 568 836 IF(lwp) THEN ! Namelist print 569 WRITE(numout,*)570 WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'571 WRITE(numout,*) '~~~~~~~~~~~~'572 837 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 573 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene574 838 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 575 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 839 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 840 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT 841 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT 576 842 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 577 843 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 844 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 578 845 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 579 846 ENDIF 847 848 IF( ln_dynvor_msk ) CALL ctl_stop( 'dyn_vor_init: masked vorticity is not currently not available') 580 849 581 850 !!gm this should be removed when choosing a unique strategy for fmask at the coast … … 586 855 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 587 856 DO jk = 1, jpk 588 DO jj = 2, jpjm1589 DO ji = 2, jpim1590 IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp )&591 fmask(ji,jj,jk) = 1._wp857 DO jj = 1, jpjm1 858 DO ji = 1, jpim1 859 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 860 & + tmask(ji,jj ,jk) + tmask(ji+1,jj+1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 592 861 END DO 593 862 END DO 594 863 END DO 595 596 597 864 ! 865 CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask 866 ! 598 867 ENDIF 599 868 !!gm end 600 869 601 870 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 602 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 603 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 604 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 605 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 871 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 872 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 873 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 874 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF 875 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 876 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 606 877 ! 607 878 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) … … 622 893 nrvm = np_MET ! metric term 623 894 ntot = np_CME ! Coriolis + metric term 895 ! 896 SELECT CASE( nvor_scheme ) ! pre-computed gradients for the metric term: 897 CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2 898 ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) ) 899 DO jj = 2, jpjm1 900 DO ji = 2, jpim1 901 di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp 902 dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp 903 END DO 904 END DO 905 CALL lbc_lnk_multi( di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. ) ! Lateral boundary conditions 906 ! 907 CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f) 908 ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) ) 909 DO jj = 1, jpjm1 910 DO ji = 1, jpim1 911 di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 912 dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 913 END DO 914 END DO 915 CALL lbc_lnk_multi( di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. ) ! Lateral boundary conditions 916 END SELECT 917 ! 624 918 END SELECT 625 919 … … 627 921 WRITE(numout,*) 628 922 SELECT CASE( nvor_scheme ) 629 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme' 630 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme' 631 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme' 632 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme' 923 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 924 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 925 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 926 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 927 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 928 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 633 929 END SELECT 634 930 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
