- Timestamp:
- 2015-11-20T09:39:06+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5038 r5901 16 16 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 17 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 18 !! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory 18 19 !!---------------------------------------------------------------------- 19 20 … … 22 23 !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) 23 24 !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) 24 !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)25 25 !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) 26 26 !! dyn_vor_init : set and control of the different vorticity option … … 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 ! 34 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 36 USE prtctl ! Print control … … 44 45 45 46 PUBLIC dyn_vor ! routine called by step.F90 46 PUBLIC dyn_vor_init ! routine called by opa.F9047 PUBLIC dyn_vor_init ! routine called by nemogcm.F90 47 48 48 49 ! !!* Namelist namdyn_vor: vorticity term 49 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme 50 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme 51 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme 52 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme 53 LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) 54 55 INTEGER :: nvor = 0 ! type of vorticity trend used 56 INTEGER :: ncor = 1 ! coriolis 57 INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term 58 INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term 59 50 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE) 51 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 52 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 53 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN) 54 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 55 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 56 57 INTEGER :: nvor_scheme ! choice of the type of advection scheme 58 ! ! associated indices: 59 INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme 60 INTEGER, PUBLIC, PARAMETER :: np_ENS = 2 ! ENS scheme 61 INTEGER, PUBLIC, PARAMETER :: np_MIX = 3 ! MIX scheme 62 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 63 64 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 65 ! ! associated indices: 66 INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 67 INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity 68 INTEGER, PARAMETER :: np_MET = 3 ! metric term 69 INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 70 INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term 71 72 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 73 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 74 REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12 75 60 76 !! * Substitutions 61 77 # include "domzgr_substitute.h90" 62 78 # include "vectopt_loop_substitute.h90" 63 79 !!---------------------------------------------------------------------- 64 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)80 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 65 81 !! $Id$ 66 82 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 87 103 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 88 104 ! 89 ! ! vorticity term 90 SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend 91 ! 92 CASE ( -1 ) ! esopa: test all possibility with control print 93 CALL vor_ene( kt, ntot, ua, va ) 94 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 95 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 96 CALL vor_ens( kt, ntot, ua, va ) 97 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 98 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 99 CALL vor_mix( kt ) 100 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 101 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 102 CALL vor_een( kt, ntot, ua, va ) 103 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 104 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 105 ! 106 CASE ( 0 ) ! energy conserving scheme 107 IF( l_trddyn ) THEN 108 ztrdu(:,:,:) = ua(:,:,:) 109 ztrdv(:,:,:) = va(:,:,:) 110 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 105 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 106 ! 107 CASE ( np_ENE ) !* energy conserving scheme 108 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 109 ztrdu(:,:,:) = ua(:,:,:) 110 ztrdv(:,:,:) = va(:,:,:) 111 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 111 112 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 112 113 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 114 115 ztrdu(:,:,:) = ua(:,:,:) 115 116 ztrdv(:,:,:) = va(:,:,:) 116 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend117 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend 117 118 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 118 119 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 119 120 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 120 121 ELSE 121 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity 122 ENDIF 123 ! 124 CASE ( 1 ) !enstrophy conserving scheme125 IF( l_trddyn ) THEN126 ztrdu(:,:,:) = ua(:,:,:) 127 ztrdv(:,:,:) = va(:,:,:) 128 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend122 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity trend 123 ENDIF 124 ! 125 CASE ( np_ENS ) !* enstrophy conserving scheme 126 IF( l_trddyn ) THEN ! trend diagnostics: splitthe trend in two 127 ztrdu(:,:,:) = ua(:,:,:) 128 ztrdv(:,:,:) = va(:,:,:) 129 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend 129 130 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 130 131 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 132 133 ztrdu(:,:,:) = ua(:,:,:) 133 134 ztrdv(:,:,:) = va(:,:,:) 134 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend135 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend 135 136 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 136 137 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 137 138 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 138 139 ELSE 139 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity 140 ENDIF 141 ! 142 CASE ( 2 ) !mixed ene-ens scheme143 IF( l_trddyn ) THEN144 ztrdu(:,:,:) = ua(:,:,:) 145 ztrdv(:,:,:) = va(:,:,:) 146 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)140 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity trend 141 ENDIF 142 ! 143 CASE ( np_MIX ) !* mixed ene-ens scheme 144 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 145 ztrdu(:,:,:) = ua(:,:,:) 146 ztrdv(:,:,:) = va(:,:,:) 147 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 147 148 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 148 149 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 150 151 ztrdu(:,:,:) = ua(:,:,:) 151 152 ztrdv(:,:,:) = va(:,:,:) 152 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)153 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 153 154 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 154 155 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 155 156 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 157 ELSE 157 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) 158 ENDIF 159 ! 160 CASE ( 3 ) ! energy and enstrophy conserving scheme 161 IF( l_trddyn ) THEN 162 ztrdu(:,:,:) = ua(:,:,:) 163 ztrdv(:,:,:) = va(:,:,:) 164 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 158 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 159 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 160 ENDIF 161 ! 162 CASE ( np_EEN ) !* energy and enstrophy conserving scheme 163 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 164 ztrdu(:,:,:) = ua(:,:,:) 165 ztrdv(:,:,:) = va(:,:,:) 166 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 165 167 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 168 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 168 170 ztrdu(:,:,:) = ua(:,:,:) 169 171 ztrdv(:,:,:) = va(:,:,:) 170 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend172 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend 171 173 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 174 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 175 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 174 176 ELSE 175 CALL vor_een( kt, ntot, ua, va ) ! total vorticity 177 CALL vor_een( kt, ntot, ua, va ) ! total vorticity trend 176 178 ENDIF 177 179 ! … … 197 199 !! 198 200 !! ** Method : Trend evaluated using now fields (centered in time) 199 !! and the Sadourny (1975) flux form formulation : conserves the 200 !! horizontal kinetic energy. 201 !! The trend of the vorticity term is given by: 202 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 203 !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] 204 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] 205 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 206 !! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ] 207 !! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ] 208 !! Add this trend to the general momentum trend (ua,va): 209 !! (ua,va) = (ua,va) + ( voru , vorv ) 201 !! and the Sadourny (1975) flux form formulation : conserves the 202 !! horizontal kinetic energy. 203 !! The general trend of momentum is increased due to the vorticity 204 !! term which is given by: 205 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v vn) ] 206 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u un) ] 207 !! where rvor is the relative vorticity 210 208 !! 211 209 !! ** Action : - Update (ua,va) with the now vorticity term trend … … 213 211 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 214 212 !!---------------------------------------------------------------------- 215 !216 213 INTEGER , INTENT(in ) :: kt ! ocean time-step index 217 214 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 220 217 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 221 218 ! 222 INTEGER :: ji, jj, jk ! dummy loop indices223 REAL(wp) :: zx1, zy1, z fact2, zx2, zy2 ! local scalars224 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz219 INTEGER :: ji, jj, jk ! dummy loop indices 220 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 221 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! 2D workspace 225 222 !!---------------------------------------------------------------------- 226 223 ! … … 234 231 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 235 232 ENDIF 236 237 zfact2 = 0.5 * 0.5 ! Local constant initialization 238 239 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 233 ! 240 234 ! ! =============== 241 235 DO jk = 1, jpkm1 ! Horizontal slab 242 236 ! ! =============== 243 237 ! 244 ! Potential vorticity and horizontal fluxes 245 ! ----------------------------------------- 246 SELECT CASE( kvor ) ! vorticity considered 247 CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) 248 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity 249 CASE ( 3 ) ! metric term 238 SELECT CASE( kvor ) !== vorticity considered ==! 239 CASE ( np_COR ) !* Coriolis (planetary vorticity) 240 zwz(:,:) = ff(:,:) 241 CASE ( np_RVO ) !* relative vorticity 242 DO jj = 1, jpjm1 243 DO ji = 1, fs_jpim1 ! vector opt. 244 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 245 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 246 END DO 247 END DO 248 CASE ( np_MET ) !* metric term 250 249 DO jj = 1, jpjm1 251 250 DO ji = 1, fs_jpim1 ! vector opt. 252 251 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 253 252 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 254 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 255 END DO 256 END DO 257 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) 258 CASE ( 5 ) ! total (coriolis + metric) 259 DO jj = 1, jpjm1 260 DO ji = 1, fs_jpim1 ! vector opt. 261 zwz(ji,jj) = ( ff (ji,jj) & 262 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 263 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 264 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 265 & ) 266 END DO 267 END DO 253 & * 0.5 * r1_e1e2f(ji,jj) 254 END DO 255 END DO 256 CASE ( np_CRV ) !* Coriolis + relative vorticity 257 DO jj = 1, jpjm1 258 DO ji = 1, fs_jpim1 ! vector opt. 259 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 260 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 261 & * r1_e1e2f(ji,jj) 262 END DO 263 END DO 264 CASE ( np_CME ) !* Coriolis + metric 265 DO jj = 1, jpjm1 266 DO ji = 1, fs_jpim1 ! vector opt. 267 zwz(ji,jj) = ff(ji,jj) & 268 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 269 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 270 & * 0.5 * r1_e1e2f(ji,jj) 271 END DO 272 END DO 273 CASE DEFAULT ! error 274 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 268 275 END SELECT 276 ! 277 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 278 DO jj = 1, jpjm1 279 DO ji = 1, fs_jpim1 ! vector opt. 280 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 281 END DO 282 END DO 283 ENDIF 269 284 270 285 IF( ln_sco ) THEN … … 276 291 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 277 292 ENDIF 278 279 ! Compute and add the vorticity term trend 280 ! ---------------------------------------- 293 ! !== compute and add the vorticity term trend =! 281 294 DO jj = 2, jpjm1 282 295 DO ji = fs_2, fs_jpim1 ! vector opt. … … 285 298 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 286 299 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 287 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 /e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )288 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 /e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )300 pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 301 pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 289 302 END DO 290 303 END DO … … 299 312 300 313 301 SUBROUTINE vor_mix( kt )302 !!----------------------------------------------------------------------303 !! *** ROUTINE vor_mix ***304 !!305 !! ** Purpose : Compute the now total vorticity trend and add it to306 !! the general trend of the momentum equation.307 !!308 !! ** Method : Trend evaluated using now fields (centered in time)309 !! Mixte formulation : conserves the potential enstrophy of a hori-310 !! zontally non-divergent flow for (rotzu x uh), the relative vor-311 !! ticity term and the horizontal kinetic energy for (f x uh), the312 !! coriolis term. the now trend of the vorticity term is given by:313 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:314 !! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]315 !! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]316 !! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]317 !! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]318 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:319 !! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]320 !! +1/e1u mj-1[ f mi(e1v vn) ]321 !! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]322 !! +1/e2v mi-1[ f mj(e2u un) ]323 !! Add this now trend to the general momentum trend (ua,va):324 !! (ua,va) = (ua,va) + ( voru , vorv )325 !!326 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend327 !!328 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.329 !!----------------------------------------------------------------------330 !331 INTEGER, INTENT(in) :: kt ! ocean timestep index332 !333 INTEGER :: ji, jj, jk ! dummy loop indices334 REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars335 REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -336 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww337 !!----------------------------------------------------------------------338 !339 IF( nn_timing == 1 ) CALL timing_start('vor_mix')340 !341 CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )342 !343 IF( kt == nit000 ) THEN344 IF(lwp) WRITE(numout,*)345 IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'346 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'347 ENDIF348 349 zfact1 = 0.5 * 0.25 ! Local constant initialization350 zfact2 = 0.5 * 0.5351 352 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )353 ! ! ===============354 DO jk = 1, jpkm1 ! Horizontal slab355 ! ! ===============356 !357 ! Relative and planetary potential vorticity and horizontal fluxes358 ! ----------------------------------------------------------------359 IF( ln_sco ) THEN360 IF( ln_dynadv_vec ) THEN361 zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)362 ELSE363 DO jj = 1, jpjm1364 DO ji = 1, fs_jpim1 ! vector opt.365 zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &366 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &367 & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )368 END DO369 END DO370 ENDIF371 zwz(:,:) = ff (:,:) / fse3f(:,:,jk)372 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)373 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)374 ELSE375 IF( ln_dynadv_vec ) THEN376 zww(:,:) = rotn(:,:,jk)377 ELSE378 DO jj = 1, jpjm1379 DO ji = 1, fs_jpim1 ! vector opt.380 zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &381 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &382 & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )383 END DO384 END DO385 ENDIF386 zwz(:,:) = ff (:,:)387 zwx(:,:) = e2u(:,:) * un(:,:,jk)388 zwy(:,:) = e1v(:,:) * vn(:,:,jk)389 ENDIF390 391 ! Compute and add the vorticity term trend392 ! ----------------------------------------393 DO jj = 2, jpjm1394 DO ji = fs_2, fs_jpim1 ! vector opt.395 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)396 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)397 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)398 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)399 ! enstrophy conserving formulation for relative vorticity term400 zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )401 zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )402 ! energy conserving formulation for planetary vorticity term403 zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )404 zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )405 ! mixed vorticity trend added to the momentum trends406 ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua407 va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva408 END DO409 END DO410 ! ! ===============411 END DO ! End of slab412 ! ! ===============413 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )414 !415 IF( nn_timing == 1 ) CALL timing_stop('vor_mix')416 !417 END SUBROUTINE vor_mix418 419 420 314 SUBROUTINE vor_ens( kt, kvor, pua, pva ) 421 315 !!---------------------------------------------------------------------- … … 429 323 !! potential enstrophy of a horizontally non-divergent flow. the 430 324 !! trend of the vorticity term is given by: 431 !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: 432 !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 433 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 434 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 435 !! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ] 436 !! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ] 325 !! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 326 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 437 327 !! Add this trend to the general momentum trend (ua,va): 438 328 !! (ua,va) = (ua,va) + ( voru , vorv ) … … 442 332 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 443 333 !!---------------------------------------------------------------------- 444 !445 334 INTEGER , INTENT(in ) :: kt ! ocean time-step index 446 335 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 449 338 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 450 339 ! 451 INTEGER :: ji, jj, jk 452 REAL(wp) :: z fact1, zuav, zvau ! temporaryscalars453 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww340 INTEGER :: ji, jj, jk ! dummy loop indices 341 REAL(wp) :: zuav, zvau ! local scalars 342 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww ! 2D workspace 454 343 !!---------------------------------------------------------------------- 455 344 ! … … 463 352 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 464 353 ENDIF 465 466 zfact1 = 0.5 * 0.25 ! Local constant initialization467 468 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )469 354 ! ! =============== 470 355 DO jk = 1, jpkm1 ! Horizontal slab 471 356 ! ! =============== 472 357 ! 473 ! Potential vorticity and horizontal fluxes 474 ! ----------------------------------------- 475 SELECT CASE( kvor ) ! vorticity considered 476 CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) 477 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity 478 CASE ( 3 ) ! metric term 358 SELECT CASE( kvor ) !== vorticity considered ==! 359 CASE ( np_COR ) !* Coriolis (planetary vorticity) 360 zwz(:,:) = ff(:,:) 361 CASE ( np_RVO ) !* relative vorticity 362 DO jj = 1, jpjm1 363 DO ji = 1, fs_jpim1 ! vector opt. 364 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 365 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 366 END DO 367 END DO 368 CASE ( np_MET ) !* metric term 479 369 DO jj = 1, jpjm1 480 370 DO ji = 1, fs_jpim1 ! vector opt. 481 371 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 482 372 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 483 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 484 END DO 485 END DO 486 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) 487 CASE ( 5 ) ! total (coriolis + metric) 488 DO jj = 1, jpjm1 489 DO ji = 1, fs_jpim1 ! vector opt. 490 zwz(ji,jj) = ( ff (ji,jj) & 491 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 492 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 493 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 494 & ) 495 END DO 496 END DO 373 & * 0.5 * r1_e1e2f(ji,jj) 374 END DO 375 END DO 376 CASE ( np_CRV ) !* Coriolis + relative vorticity 377 DO jj = 1, jpjm1 378 DO ji = 1, fs_jpim1 ! vector opt. 379 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 380 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 381 & * r1_e1e2f(ji,jj) 382 END DO 383 END DO 384 CASE ( np_CME ) !* Coriolis + metric 385 DO jj = 1, jpjm1 386 DO ji = 1, fs_jpim1 ! vector opt. 387 zwz(ji,jj) = ff(ji,jj) & 388 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 389 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 390 & * 0.5 * r1_e1e2f(ji,jj) 391 END DO 392 END DO 393 CASE DEFAULT ! error 394 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 497 395 END SELECT 498 396 ! 499 IF( ln_sco ) THEN 500 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 501 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 502 zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 503 zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 504 zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 505 END DO 506 END DO 397 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 398 DO jj = 1, jpjm1 399 DO ji = 1, fs_jpim1 ! vector opt. 400 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 401 END DO 402 END DO 403 ENDIF 404 ! 405 IF( ln_sco ) THEN !== horizontal fluxes ==! 406 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 407 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 408 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 507 409 ELSE 508 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 509 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 510 zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 511 zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 512 END DO 513 END DO 514 ENDIF 515 ! 516 ! Compute and add the vorticity term trend 517 ! ---------------------------------------- 410 zwx(:,:) = e2u(:,:) * un(:,:,jk) 411 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 412 ENDIF 413 ! !== compute and add the vorticity term trend =! 518 414 DO jj = 2, jpjm1 519 415 DO ji = fs_2, fs_jpim1 ! vector opt. 520 zuav = zfact1 /e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &521 & 522 zvau =- zfact1 /e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &523 & 416 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 417 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 418 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 419 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 524 420 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 525 421 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 553 449 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 554 450 !!---------------------------------------------------------------------- 555 !556 451 INTEGER , INTENT(in ) :: kt ! ocean time-step index 557 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 558 ! ! =nrvm (relative vorticity or metric) 452 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 559 453 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 560 454 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 561 !! 562 INTEGER :: ji, jj, jk ! dummy loop indices 563 INTEGER :: ierr ! local integer 564 REAL(wp) :: zfac12, zua, zva ! local scalars 565 REAL(wp) :: zmsk, ze3 ! local scalars 566 ! ! 3D workspace 567 REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz 568 REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse 569 #if defined key_vvl 570 REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T) 571 #else 572 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all 573 #endif 455 ! 456 INTEGER :: ji, jj, jk ! dummy loop indices 457 INTEGER :: ierr ! local integer 458 REAL(wp) :: zua, zva ! local scalars 459 REAL(wp) :: zmsk, ze3 ! local scalars 460 ! 461 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, z1_e3f 462 REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse 574 463 !!---------------------------------------------------------------------- 575 464 ! 576 465 IF( nn_timing == 1 ) CALL timing_start('vor_een') 577 466 ! 578 CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz ) 579 CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) 580 #if defined key_vvl 581 CALL wrk_alloc( jpi, jpj, jpk, ze3f ) 582 #endif 467 CALL wrk_alloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 468 CALL wrk_alloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 583 469 ! 584 470 IF( kt == nit000 ) THEN … … 586 472 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 587 473 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 588 #if ! defined key_vvl589 IF( .NOT.ALLOCATED(ze3f) ) THEN590 ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )591 IF( lk_mpp ) CALL mpp_sum ( ierr )592 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )593 ENDIF594 ze3f(:,:,:) = 0._wp595 #endif596 474 ENDIF 597 598 IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 599 600 IF( ln_dynvor_een_old ) THEN ! original formulation 601 DO jk = 1, jpk 602 DO jj = 1, jpjm1 603 DO ji = 1, jpim1 604 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 605 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 606 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = 4.0_wp / ze3 607 END DO 608 END DO 609 END DO 610 ELSE ! new formulation from NEMO 3.6 611 DO jk = 1, jpk 612 DO jj = 1, jpjm1 613 DO ji = 1, jpim1 614 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 615 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 616 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 617 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 618 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3 619 END DO 620 END DO 621 END DO 622 ENDIF 623 624 CALL lbc_lnk( ze3f, 'F', 1. ) 625 ENDIF 626 627 zfac12 = 1._wp / 12._wp ! Local constant initialization 628 629 630 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 475 ! 631 476 ! ! =============== 632 477 DO jk = 1, jpkm1 ! Horizontal slab 633 478 ! ! =============== 634 635 ! Potential vorticity and horizontal fluxes 636 ! ----------------------------------------- 637 SELECT CASE( kvor ) ! vorticity considered 638 CASE ( 1 ) ! planetary vorticity (Coriolis) 639 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 640 CASE ( 2 ) ! relative vorticity 641 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 642 CASE ( 3 ) ! metric term 479 ! 480 SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point 481 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 482 DO jj = 1, jpjm1 483 DO ji = 1, fs_jpim1 ! vector opt. 484 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 485 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 486 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4.0_wp / ze3 487 ELSE ; z1_e3f(ji,jj) = 0.0_wp 488 ENDIF 489 END DO 490 END DO 491 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 492 DO jj = 1, jpjm1 493 DO ji = 1, fs_jpim1 ! vector opt. 494 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 495 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 496 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 497 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 498 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3 499 ELSE ; z1_e3f(ji,jj) = 0.0_wp 500 ENDIF 501 END DO 502 END DO 503 END SELECT 504 ! 505 SELECT CASE( kvor ) !== vorticity considered ==! 506 CASE ( np_COR ) !* Coriolis (planetary vorticity) 507 DO jj = 1, jpjm1 508 DO ji = 1, fs_jpim1 ! vector opt. 509 zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj) 510 END DO 511 END DO 512 CASE ( np_RVO ) !* relative vorticity 513 DO jj = 1, jpjm1 514 DO ji = 1, fs_jpim1 ! vector opt. 515 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 516 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 517 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 518 END DO 519 END DO 520 CASE ( np_MET ) !* metric term 643 521 DO jj = 1, jpjm1 644 522 DO ji = 1, fs_jpim1 ! vector opt. 645 523 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 646 524 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 647 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) 648 END DO 649 END DO 650 CALL lbc_lnk( zwz, 'F', 1. ) 651 CASE ( 4 ) ! total (relative + planetary vorticity) 652 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 653 CASE ( 5 ) ! total (coriolis + metric) 654 DO jj = 1, jpjm1 655 DO ji = 1, fs_jpim1 ! vector opt. 656 zwz(ji,jj) = ( ff (ji,jj) & 657 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 658 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 659 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 660 & ) * ze3f(ji,jj,jk) 661 END DO 662 END DO 663 CALL lbc_lnk( zwz, 'F', 1. ) 525 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 526 END DO 527 END DO 528 CASE ( np_CRV ) !* Coriolis + relative vorticity 529 DO jj = 1, jpjm1 530 DO ji = 1, fs_jpim1 ! vector opt. 531 zwz(ji,jj) = ( ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 532 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 533 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 534 END DO 535 END DO 536 CASE ( np_CME ) !* Coriolis + metric 537 DO jj = 1, jpjm1 538 DO ji = 1, fs_jpim1 ! vector opt. 539 zwz(ji,jj) = ( ff(ji,jj) & 540 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 541 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 542 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 543 END DO 544 END DO 545 CASE DEFAULT ! error 546 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 664 547 END SELECT 665 548 ! 549 CALL lbc_lnk( zwz, 'F', 1. ) 550 ! 551 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 552 DO jj = 1, jpjm1 553 DO ji = 1, fs_jpim1 ! vector opt. 554 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 555 END DO 556 END DO 557 ENDIF 558 ! 559 ! !== horizontal fluxes ==! 666 560 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 667 561 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 668 562 669 ! Compute and add the vorticity term trend 670 ! ---------------------------------------- 563 ! !== compute and add the vorticity term trend =! 671 564 jj = 2 672 565 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 673 DO ji = 2, jpi 566 DO ji = 2, jpi ! split in 2 parts due to vector opt. 674 567 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 675 568 ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 687 580 DO jj = 2, jpjm1 688 581 DO ji = fs_2, fs_jpim1 ! vector opt. 689 zua = + zfac12 /e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &690 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )691 zva = - zfac12 /e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &692 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )582 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 583 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 584 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 585 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 693 586 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 694 587 pva(ji,jj,jk) = pva(ji,jj,jk) + zva … … 698 591 END DO ! End of slab 699 592 ! ! =============== 700 CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz ) 701 CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) 702 #if defined key_vvl 703 CALL wrk_dealloc( jpi, jpj, jpk, ze3f ) 704 #endif 593 ! 594 CALL wrk_dealloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 595 CALL wrk_dealloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 705 596 ! 706 597 IF( nn_timing == 1 ) CALL timing_stop('vor_een') … … 720 611 INTEGER :: ios ! Local integer output status for namelist read 721 612 !! 722 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old613 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk 723 614 !!---------------------------------------------------------------------- 724 615 … … 737 628 WRITE(numout,*) '~~~~~~~~~~~~' 738 629 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 739 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 740 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 741 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 742 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 743 WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old 630 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 631 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 632 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 633 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 634 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 635 WRITE(numout,*) ' masked (=1) or unmasked(=0) vorticity ln_dynvor_msk = ', ln_dynvor_msk 744 636 ENDIF 745 637 638 !!gm this should be removed when choosing a unique strategy for fmask at the coast 746 639 ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 747 640 ! at angles with three ocean points and one land point 641 IF(lwp) WRITE(numout,*) 642 IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 748 643 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 749 644 DO jk = 1, jpk … … 759 654 ! 760 655 ENDIF 761 762 ioptio = 0 ! Control of vorticity scheme options 763 IF( ln_dynvor_ene ) ioptio = ioptio + 1 764 IF( ln_dynvor_ens ) ioptio = ioptio + 1 765 IF( ln_dynvor_mix ) ioptio = ioptio + 1 766 IF( ln_dynvor_een ) ioptio = ioptio + 1 767 IF( ln_dynvor_een_old ) ioptio = ioptio + 1 768 IF( lk_esopa ) ioptio = 1 769 656 !!gm end 657 658 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 659 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 660 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 661 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 662 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 663 ! 770 664 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 771 772 ! ! Set nvor (type of scheme for vorticity) 773 IF( ln_dynvor_ene ) nvor = 0 774 IF( ln_dynvor_ens ) nvor = 1 775 IF( ln_dynvor_mix ) nvor = 2 776 IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3 777 IF( lk_esopa ) nvor = -1 778 779 ! ! Set ncor, nrvm, ntot (type of vorticity) 780 IF(lwp) WRITE(numout,*) 781 ncor = 1 665 ! 666 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) 667 ncor = np_COR 782 668 IF( ln_dynadv_vec ) THEN 783 669 IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity' 784 nrvm = 2785 ntot = 4670 nrvm = np_RVO ! relative vorticity 671 ntot = np_CRV ! relative + planetary vorticity 786 672 ELSE 787 673 IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term' 788 nrvm = 3789 ntot = 5674 nrvm = np_MET ! metric term 675 ntot = np_CME ! Coriolis + metric term 790 676 ENDIF 791 677 792 678 IF(lwp) THEN ! Print the choice 793 679 WRITE(numout,*) 794 IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme' 795 IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme' 796 IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme' 797 IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme' 798 IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options' 680 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' vorticity scheme ==>> energy conserving scheme' 681 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' vorticity scheme ==>> enstrophy conserving scheme' 682 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' vorticity scheme ==>> mixed enstrophy/energy conserving scheme' 683 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' vorticity scheme ==>> energy and enstrophy conserving scheme' 799 684 ENDIF 800 685 !
Note: See TracChangeset
for help on using the changeset viewer.
