Changeset 5974 for branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
- Timestamp:
- 2015-12-02T11:52:05+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5029 r5974 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 USE c1d ! 1D vertical configuration 35 ! 34 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 37 USE prtctl ! Print control … … 44 46 45 47 PUBLIC dyn_vor ! routine called by step.F90 46 PUBLIC dyn_vor_init ! routine called by opa.F9048 PUBLIC dyn_vor_init ! routine called by nemogcm.F90 47 49 48 50 ! !!* 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 51 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE) 52 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 53 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 54 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN) 55 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 56 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 57 58 INTEGER :: nvor_scheme ! choice of the type of advection scheme 59 ! ! associated indices: 60 INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme 61 INTEGER, PUBLIC, PARAMETER :: np_ENS = 2 ! ENS scheme 62 INTEGER, PUBLIC, PARAMETER :: np_MIX = 3 ! MIX scheme 63 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 64 65 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 66 ! ! associated indices: 67 INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 68 INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity 69 INTEGER, PARAMETER :: np_MET = 3 ! metric term 70 INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 71 INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term 72 73 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 74 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 75 REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12 76 60 77 !! * Substitutions 61 78 # include "domzgr_substitute.h90" 62 79 # include "vectopt_loop_substitute.h90" 63 80 !!---------------------------------------------------------------------- 64 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)81 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 65 82 !! $Id$ 66 83 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 87 104 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 88 105 ! 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 106 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 107 ! 108 CASE ( np_ENE ) !* energy conserving scheme 109 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 110 ztrdu(:,:,:) = ua(:,:,:) 111 ztrdv(:,:,:) = va(:,:,:) 112 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 111 113 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 112 114 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 114 116 ztrdu(:,:,:) = ua(:,:,:) 115 117 ztrdv(:,:,:) = va(:,:,:) 116 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend118 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend 117 119 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 118 120 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 119 121 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 120 122 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 trend123 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity trend 124 ENDIF 125 ! 126 CASE ( np_ENS ) !* enstrophy conserving scheme 127 IF( l_trddyn ) THEN ! trend diagnostics: splitthe trend in two 128 ztrdu(:,:,:) = ua(:,:,:) 129 ztrdv(:,:,:) = va(:,:,:) 130 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend 129 131 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 130 132 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 132 134 ztrdu(:,:,:) = ua(:,:,:) 133 135 ztrdv(:,:,:) = va(:,:,:) 134 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend136 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend 135 137 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 136 138 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 137 139 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 138 140 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)141 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity trend 142 ENDIF 143 ! 144 CASE ( np_MIX ) !* mixed ene-ens scheme 145 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 146 ztrdu(:,:,:) = ua(:,:,:) 147 ztrdv(:,:,:) = va(:,:,:) 148 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 147 149 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 148 150 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 150 152 ztrdu(:,:,:) = ua(:,:,:) 151 153 ztrdv(:,:,:) = va(:,:,:) 152 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)154 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 153 155 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 154 156 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 155 157 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 158 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 159 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 160 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 161 ENDIF 162 ! 163 CASE ( np_EEN ) !* energy and enstrophy conserving scheme 164 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 165 ztrdu(:,:,:) = ua(:,:,:) 166 ztrdv(:,:,:) = va(:,:,:) 167 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 165 168 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 169 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 168 171 ztrdu(:,:,:) = ua(:,:,:) 169 172 ztrdv(:,:,:) = va(:,:,:) 170 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend173 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend 171 174 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 175 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 176 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 174 177 ELSE 175 CALL vor_een( kt, ntot, ua, va ) ! total vorticity 178 CALL vor_een( kt, ntot, ua, va ) ! total vorticity trend 176 179 ENDIF 177 180 ! … … 197 200 !! 198 201 !! ** 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 ) 202 !! and the Sadourny (1975) flux form formulation : conserves the 203 !! horizontal kinetic energy. 204 !! The general trend of momentum is increased due to the vorticity 205 !! term which is given by: 206 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v vn) ] 207 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u un) ] 208 !! where rvor is the relative vorticity 210 209 !! 211 210 !! ** Action : - Update (ua,va) with the now vorticity term trend … … 213 212 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 214 213 !!---------------------------------------------------------------------- 215 !216 214 INTEGER , INTENT(in ) :: kt ! ocean time-step index 217 215 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 220 218 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 221 219 ! 222 INTEGER :: ji, jj, jk ! dummy loop indices223 REAL(wp) :: zx1, zy1, z fact2, zx2, zy2 ! local scalars224 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz220 INTEGER :: ji, jj, jk ! dummy loop indices 221 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 222 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! 2D workspace 225 223 !!---------------------------------------------------------------------- 226 224 ! … … 234 232 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 235 233 ENDIF 236 237 zfact2 = 0.5 * 0.5 ! Local constant initialization 238 239 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 234 ! 240 235 ! ! =============== 241 236 DO jk = 1, jpkm1 ! Horizontal slab 242 237 ! ! =============== 243 238 ! 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 239 SELECT CASE( kvor ) !== vorticity considered ==! 240 CASE ( np_COR ) !* Coriolis (planetary vorticity) 241 zwz(:,:) = ff(:,:) 242 CASE ( np_RVO ) !* relative vorticity 243 DO jj = 1, jpjm1 244 DO ji = 1, fs_jpim1 ! vector opt. 245 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 246 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 247 END DO 248 END DO 249 CASE ( np_MET ) !* metric term 250 250 DO jj = 1, jpjm1 251 251 DO ji = 1, fs_jpim1 ! vector opt. 252 252 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 253 253 & - ( 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 254 & * 0.5 * r1_e1e2f(ji,jj) 255 END DO 256 END DO 257 CASE ( np_CRV ) !* Coriolis + relative vorticity 258 DO jj = 1, jpjm1 259 DO ji = 1, fs_jpim1 ! vector opt. 260 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 261 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 262 & * r1_e1e2f(ji,jj) 263 END DO 264 END DO 265 CASE ( np_CME ) !* Coriolis + metric 266 DO jj = 1, jpjm1 267 DO ji = 1, fs_jpim1 ! vector opt. 268 zwz(ji,jj) = ff(ji,jj) & 269 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 270 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 271 & * 0.5 * r1_e1e2f(ji,jj) 272 END DO 273 END DO 274 CASE DEFAULT ! error 275 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 268 276 END SELECT 277 ! 278 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 279 DO jj = 1, jpjm1 280 DO ji = 1, fs_jpim1 ! vector opt. 281 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 282 END DO 283 END DO 284 ENDIF 269 285 270 286 IF( ln_sco ) THEN … … 276 292 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 277 293 ENDIF 278 279 ! Compute and add the vorticity term trend 280 ! ---------------------------------------- 294 ! !== compute and add the vorticity term trend =! 281 295 DO jj = 2, jpjm1 282 296 DO ji = fs_2, fs_jpim1 ! vector opt. … … 285 299 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 286 300 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 )301 pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 302 pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 289 303 END DO 290 304 END DO … … 299 313 300 314 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 315 SUBROUTINE vor_ens( kt, kvor, pua, pva ) 421 316 !!---------------------------------------------------------------------- … … 429 324 !! potential enstrophy of a horizontally non-divergent flow. the 430 325 !! 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) ] 326 !! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 327 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 437 328 !! Add this trend to the general momentum trend (ua,va): 438 329 !! (ua,va) = (ua,va) + ( voru , vorv ) … … 442 333 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 443 334 !!---------------------------------------------------------------------- 444 !445 335 INTEGER , INTENT(in ) :: kt ! ocean time-step index 446 336 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 449 339 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 450 340 ! 451 INTEGER :: ji, jj, jk 452 REAL(wp) :: z fact1, zuav, zvau ! temporaryscalars453 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww341 INTEGER :: ji, jj, jk ! dummy loop indices 342 REAL(wp) :: zuav, zvau ! local scalars 343 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww ! 2D workspace 454 344 !!---------------------------------------------------------------------- 455 345 ! … … 463 353 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 464 354 ENDIF 465 466 zfact1 = 0.5 * 0.25 ! Local constant initialization467 468 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )469 355 ! ! =============== 470 356 DO jk = 1, jpkm1 ! Horizontal slab 471 357 ! ! =============== 472 358 ! 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 359 SELECT CASE( kvor ) !== vorticity considered ==! 360 CASE ( np_COR ) !* Coriolis (planetary vorticity) 361 zwz(:,:) = ff(:,:) 362 CASE ( np_RVO ) !* relative vorticity 363 DO jj = 1, jpjm1 364 DO ji = 1, fs_jpim1 ! vector opt. 365 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 366 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 367 END DO 368 END DO 369 CASE ( np_MET ) !* metric term 479 370 DO jj = 1, jpjm1 480 371 DO ji = 1, fs_jpim1 ! vector opt. 481 372 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 482 373 & - ( 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 374 & * 0.5 * r1_e1e2f(ji,jj) 375 END DO 376 END DO 377 CASE ( np_CRV ) !* Coriolis + relative vorticity 378 DO jj = 1, jpjm1 379 DO ji = 1, fs_jpim1 ! vector opt. 380 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 381 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 382 & * r1_e1e2f(ji,jj) 383 END DO 384 END DO 385 CASE ( np_CME ) !* Coriolis + metric 386 DO jj = 1, jpjm1 387 DO ji = 1, fs_jpim1 ! vector opt. 388 zwz(ji,jj) = ff(ji,jj) & 389 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 390 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 391 & * 0.5 * r1_e1e2f(ji,jj) 392 END DO 393 END DO 394 CASE DEFAULT ! error 395 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 497 396 END SELECT 498 397 ! 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 398 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 399 DO jj = 1, jpjm1 400 DO ji = 1, fs_jpim1 ! vector opt. 401 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 402 END DO 403 END DO 404 ENDIF 405 ! 406 IF( ln_sco ) THEN !== horizontal fluxes ==! 407 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 408 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 409 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 507 410 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 ! ---------------------------------------- 411 zwx(:,:) = e2u(:,:) * un(:,:,jk) 412 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 413 ENDIF 414 ! !== compute and add the vorticity term trend =! 518 415 DO jj = 2, jpjm1 519 416 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 & 417 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 418 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 419 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 420 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 524 421 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 525 422 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 553 450 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 554 451 !!---------------------------------------------------------------------- 555 !556 452 INTEGER , INTENT(in ) :: kt ! ocean time-step index 557 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 558 ! ! =nrvm (relative vorticity or metric) 453 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 559 454 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 560 455 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 456 ! 457 INTEGER :: ji, jj, jk ! dummy loop indices 458 INTEGER :: ierr ! local integer 459 REAL(wp) :: zua, zva ! local scalars 460 REAL(wp) :: zmsk, ze3 ! local scalars 461 ! 462 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, z1_e3f 463 REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse 574 464 !!---------------------------------------------------------------------- 575 465 ! 576 466 IF( nn_timing == 1 ) CALL timing_start('vor_een') 577 467 ! 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 468 CALL wrk_alloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 469 CALL wrk_alloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 583 470 ! 584 471 IF( kt == nit000 ) THEN … … 586 473 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 587 474 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 475 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 ) 476 ! 631 477 ! ! =============== 632 478 DO jk = 1, jpkm1 ! Horizontal slab 633 479 ! ! =============== 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 480 ! 481 SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point 482 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 483 DO jj = 1, jpjm1 484 DO ji = 1, fs_jpim1 ! vector opt. 485 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 486 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 487 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4.0_wp / ze3 488 ELSE ; z1_e3f(ji,jj) = 0.0_wp 489 ENDIF 490 END DO 491 END DO 492 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 493 DO jj = 1, jpjm1 494 DO ji = 1, fs_jpim1 ! vector opt. 495 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 496 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 497 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 498 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 499 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3 500 ELSE ; z1_e3f(ji,jj) = 0.0_wp 501 ENDIF 502 END DO 503 END DO 504 END SELECT 505 ! 506 SELECT CASE( kvor ) !== vorticity considered ==! 507 CASE ( np_COR ) !* Coriolis (planetary vorticity) 508 DO jj = 1, jpjm1 509 DO ji = 1, fs_jpim1 ! vector opt. 510 zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj) 511 END DO 512 END DO 513 CASE ( np_RVO ) !* relative vorticity 514 DO jj = 1, jpjm1 515 DO ji = 1, fs_jpim1 ! vector opt. 516 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 517 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 518 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 519 END DO 520 END DO 521 CASE ( np_MET ) !* metric term 643 522 DO jj = 1, jpjm1 644 523 DO ji = 1, fs_jpim1 ! vector opt. 645 524 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 646 525 & - ( 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. ) 526 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 527 END DO 528 END DO 529 CASE ( np_CRV ) !* Coriolis + relative vorticity 530 DO jj = 1, jpjm1 531 DO ji = 1, fs_jpim1 ! vector opt. 532 zwz(ji,jj) = ( ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 533 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 534 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 535 END DO 536 END DO 537 CASE ( np_CME ) !* Coriolis + metric 538 DO jj = 1, jpjm1 539 DO ji = 1, fs_jpim1 ! vector opt. 540 zwz(ji,jj) = ( ff(ji,jj) & 541 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 542 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 543 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 544 END DO 545 END DO 546 CASE DEFAULT ! error 547 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 664 548 END SELECT 665 549 ! 550 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 551 DO jj = 1, jpjm1 552 DO ji = 1, fs_jpim1 ! vector opt. 553 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 554 END DO 555 END DO 556 ENDIF 557 ! 558 CALL lbc_lnk( zwz, 'F', 1. ) 559 ! 560 ! !== horizontal fluxes ==! 666 561 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 667 562 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 668 563 669 ! Compute and add the vorticity term trend 670 ! ---------------------------------------- 564 ! !== compute and add the vorticity term trend =! 671 565 jj = 2 672 566 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 673 DO ji = 2, jpi 567 DO ji = 2, jpi ! split in 2 parts due to vector opt. 674 568 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 675 569 ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 687 581 DO jj = 2, jpjm1 688 582 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 ) )583 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 584 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 585 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 586 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 693 587 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 694 588 pva(ji,jj,jk) = pva(ji,jj,jk) + zva … … 698 592 END DO ! End of slab 699 593 ! ! =============== 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 594 ! 595 CALL wrk_dealloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 596 CALL wrk_dealloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 705 597 ! 706 598 IF( nn_timing == 1 ) CALL timing_stop('vor_een') … … 720 612 INTEGER :: ios ! Local integer output status for namelist read 721 613 !! 722 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old614 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk 723 615 !!---------------------------------------------------------------------- 724 616 … … 737 629 WRITE(numout,*) '~~~~~~~~~~~~' 738 630 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 631 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 632 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 633 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 634 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 635 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 636 WRITE(numout,*) ' masked (=1) or unmasked(=0) vorticity ln_dynvor_msk = ', ln_dynvor_msk 744 637 ENDIF 745 638 639 !!gm this should be removed when choosing a unique strategy for fmask at the coast 746 640 ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 747 641 ! at angles with three ocean points and one land point 642 IF(lwp) WRITE(numout,*) 643 IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 748 644 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 749 645 DO jk = 1, jpk … … 759 655 ! 760 656 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 770 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 657 !!gm end 658 659 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 660 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 661 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 662 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 663 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 664 ! 665 IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 666 ! 667 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) 668 ncor = np_COR 782 669 IF( ln_dynadv_vec ) THEN 783 670 IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity' 784 nrvm = 2785 ntot = 4671 nrvm = np_RVO ! relative vorticity 672 ntot = np_CRV ! relative + planetary vorticity 786 673 ELSE 787 674 IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term' 788 nrvm = 3789 ntot = 5675 nrvm = np_MET ! metric term 676 ntot = np_CME ! Coriolis + metric term 790 677 ENDIF 791 678 792 679 IF(lwp) THEN ! Print the choice 793 680 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' 681 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' vorticity scheme ==>> energy conserving scheme' 682 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' vorticity scheme ==>> enstrophy conserving scheme' 683 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' vorticity scheme ==>> mixed enstrophy/energy conserving scheme' 684 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' vorticity scheme ==>> energy and enstrophy conserving scheme' 799 685 ENDIF 800 686 !
Note: See TracChangeset
for help on using the changeset viewer.
