Changeset 455 for trunk/NEMO/OPA_SRC/DYN/dynvor.F90
- Timestamp:
- 2006-05-10T18:53:54+02:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynvor.F90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynvor.F90
r418 r455 7 7 8 8 !!---------------------------------------------------------------------- 9 !! dyn_vor_enstrophy: enstrophy conserving scheme (ln_dynvor_ens=T) 10 !! dyn_vor_energy : energy conserving scheme (ln_dynvor_ene=T) 11 !! dyn_vor_mixed : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 12 !! dyn_vor_ene_ens : energy and enstrophy conserving (ln_dynvor_een=T) 13 !! dyn_vor_ctl : control of the different vorticity option 9 !! dyn_vor : Update the momentum trend with the vorticity trend 10 !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) 11 !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) 12 !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 13 !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) 14 !! vor_ctl : control of the different vorticity option 14 15 !!---------------------------------------------------------------------- 15 16 !! * Modules used … … 26 27 27 28 !! * Routine accessibility 28 PUBLIC dyn_vor_enstrophy ! routine called by step.F90 29 PUBLIC dyn_vor_energy ! routine called by step.F90 30 PUBLIC dyn_vor_mixed ! routine called by step.F90 31 PUBLIC dyn_vor_ene_ens ! routine called by step.F90 32 PUBLIC dyn_vor_ctl ! routine called by step.F90 29 PUBLIC dyn_vor ! routine called by step.F90 33 30 34 31 !! * Shared module variables … … 38 35 LOGICAL, PUBLIC :: ln_dynvor_een = .FALSE. !: energy and enstrophy conserving scheme 39 36 37 !! * module variables 38 INTEGER :: & 39 nvor = 0 ! type of vorticity trend used 40 40 41 !! * Substitutions 41 42 # include "domzgr_substitute.h90" … … 43 44 !!---------------------------------------------------------------------- 44 45 !! OPA 9.0 , LOCEAN-IPSL (2005) 45 !! $Header$46 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt47 46 !!---------------------------------------------------------------------- 48 47 49 48 CONTAINS 50 49 51 SUBROUTINE dyn_vor_energy( kt ) 52 !!---------------------------------------------------------------------- 53 !! *** ROUTINE dyn_vor_energy *** 50 SUBROUTINE dyn_vor( kt ) 51 !!---------------------------------------------------------------------- 52 !! 53 !! ** Purpose : compute the lateral ocean tracer physics. 54 !! 55 !! ** Action : - Update (ua,va) with the now vorticity term trend 56 !! - save the trends in (utrd,vtrd) in 2 parts (relative 57 !! and planetary vorticity trends) ('key_trddyn') 58 !! 59 !! History : 60 !! 9.0 ! 05-11 (G. Madec) Original code 61 !!---------------------------------------------------------------------- 62 USE oce, ONLY : ztrdu => ta, & ! use ta as 3D workspace 63 ztrdv => sa ! use sa as 3D workspace 64 65 !! * Arguments 66 INTEGER, INTENT( in ) :: kt ! ocean time-step index 67 !!---------------------------------------------------------------------- 68 69 IF( kt == nit000 ) CALL vor_ctl ! initialisation & control of options 70 71 ! ! vorticity term including Coriolis 72 SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend 73 74 CASE ( -1 ) ! esopa: test all possibility with control print 75 CALL vor_ene( kt, 'TOT', ua, va ) 76 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 77 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 78 CALL vor_ens( kt, 'TOT', ua, va ) 79 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 80 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 81 CALL vor_mix( kt ) 82 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 83 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 84 CALL vor_een( kt, 'TOT', ua, va ) 85 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 86 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 87 88 CASE ( 0 ) ! energy conserving scheme 89 IF( l_trddyn ) THEN 90 ztrdu(:,:,:) = ua(:,:,:) 91 ztrdv(:,:,:) = va(:,:,:) 92 CALL vor_ene( kt, 'VOR', ua, va ) ! relative vorticity trend 93 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 94 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 95 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 96 ztrdu(:,:,:) = ua(:,:,:) 97 ztrdv(:,:,:) = va(:,:,:) 98 CALL vor_ene( kt, 'COR', ua, va ) ! planetary vorticity trend 99 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 102 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 103 ELSE 104 CALL vor_ene( kt, 'TOT', ua, va ) ! total vorticity 105 ENDIF 106 107 CASE ( 1 ) ! enstrophy conserving scheme 108 IF( l_trddyn ) THEN 109 ztrdu(:,:,:) = ua(:,:,:) 110 ztrdv(:,:,:) = va(:,:,:) 111 CALL vor_ens( kt, 'VOR', ua, va ) ! relative vorticity trend 112 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 113 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 114 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 115 ztrdu(:,:,:) = ua(:,:,:) 116 ztrdv(:,:,:) = va(:,:,:) 117 CALL vor_ens( kt, 'COR', ua, va ) ! planetary vorticity trend 118 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 119 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 120 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 121 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 122 ELSE 123 CALL vor_ens( kt, 'TOT', ua, va ) ! total vorticity 124 ENDIF 125 126 CASE ( 2 ) ! mixed ene-ens scheme 127 IF( l_trddyn ) THEN 128 ztrdu(:,:,:) = ua(:,:,:) 129 ztrdv(:,:,:) = va(:,:,:) 130 CALL vor_ens( kt, 'VOR', ua, va ) ! relative vorticity trend (ens) 131 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 132 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 133 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 134 ztrdu(:,:,:) = ua(:,:,:) 135 ztrdv(:,:,:) = va(:,:,:) 136 CALL vor_ene( kt, 'COR', ua, va ) ! planetary vorticity trend (ene) 137 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 138 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 139 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 140 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 141 ELSE 142 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) 143 ENDIF 144 145 CASE ( 3 ) ! energy and enstrophy conserving scheme 146 IF( l_trddyn ) THEN 147 ztrdu(:,:,:) = ua(:,:,:) 148 ztrdv(:,:,:) = va(:,:,:) 149 CALL vor_een( kt, 'VOR', ua, va ) ! relative vorticity trend 150 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 151 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 152 CALL trd_mod( ztrdu, ztrdv, jpdtdrvo, 'DYN', kt ) 153 ztrdu(:,:,:) = ua(:,:,:) 154 ztrdv(:,:,:) = va(:,:,:) 155 CALL vor_een( kt, 'COR', ua, va ) ! planetary vorticity trend 156 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 157 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 158 CALL trd_mod( ztrdu, ztrdu, jpdtdpvo, 'DYN', kt ) 159 CALL trd_mod( ztrdu, ztrdv, jpdtddat, 'DYN', kt ) 160 ELSE 161 CALL vor_een( kt, 'TOT', ua, va ) ! total vorticity 162 ENDIF 163 164 END SELECT 165 166 ! ! print sum trends (used for debugging) 167 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 168 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 169 170 END SUBROUTINE dyn_vor 171 172 173 SUBROUTINE vor_ene( kt, cd_vor, pua, pva ) 174 !!---------------------------------------------------------------------- 175 !! *** ROUTINE vor_ene *** 54 176 !! 55 177 !! ** Purpose : Compute the now total vorticity trend and add it to … … 60 182 !! horizontal kinetic energy. 61 183 !! The trend of the vorticity term is given by: 62 !! * s-coordinate (l k_sco=T), the e3. are inside the derivatives:184 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 63 185 !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] 64 186 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] … … 81 203 !! 9.0 ! 04-08 (C. Talandier) New trends organization 82 204 !!---------------------------------------------------------------------- 83 !! * Modules used84 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace85 ztdva => sa ! use sa as 3D workspace86 87 205 !! * Arguments 88 206 INTEGER, INTENT( in ) :: kt ! ocean time-step index 207 CHARACTER(len=3) , INTENT( in ) :: & 208 cd_vor ! define the vorticity considered 209 ! ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 210 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 211 pua, pva ! ???? 89 212 90 213 !! * Local declarations 91 214 INTEGER :: ji, jj, jk ! dummy loop indices 92 215 REAL(wp) :: & 93 zfact2, zua, zva,& ! temporary scalars216 zfact2, & ! temporary scalars 94 217 zx1, zx2, zy1, zy2 ! " " 95 218 REAL(wp), DIMENSION(jpi,jpj) :: & 96 219 zwx, zwy, zwz ! temporary workspace 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &98 zcu, zcv ! " "99 220 !!---------------------------------------------------------------------- 100 221 101 222 IF( kt == nit000 ) THEN 102 223 IF(lwp) WRITE(numout,*) 103 IF(lwp) WRITE(numout,*) 'dyn _vor_energy: vorticity term: energy conserving scheme'104 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~'224 IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 225 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 105 226 ENDIF 106 227 107 228 ! Local constant initialization 108 229 zfact2 = 0.5 * 0.5 109 zcu (:,:,:) = 0.e0 110 zcv (:,:,:) = 0.e0 111 112 ! Save ua and va trends 113 IF( l_trddyn ) THEN 114 ztdua(:,:,:) = ua(:,:,:) 115 ztdva(:,:,:) = va(:,:,:) 116 zcu(:,:,:) = 0.e0 117 zcv(:,:,:) = 0.e0 118 ENDIF 119 230 231 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 232 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz ) 120 233 ! ! =============== 121 234 DO jk = 1, jpkm1 ! Horizontal slab … … 124 237 ! Potential vorticity and horizontal fluxes 125 238 ! ----------------------------------------- 126 IF( lk_sco ) THEN 127 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 239 SELECT CASE( cd_vor ) ! vorticity considered 240 CASE ( 'COR' ) ! planetary vorticcity (Coriolis) 241 zwz(:,:) = ff(:,:) 242 CASE ( 'VOR' ) ! relative vorticity 243 zwz(:,:) = rotn(:,:,jk) 244 CASE ( 'TOT' ) ! total vorticity (planetary + relative) 245 zwz(:,:) = rotn(:,:,jk) + ff(:,:) 246 END SELECT 247 248 IF( ln_sco ) THEN 249 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 128 250 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 129 251 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 130 252 ELSE 131 zwz(:,:) = rotn(:,:,jk) + ff(:,:)132 253 zwx(:,:) = e2u(:,:) * un(:,:,jk) 133 254 zwy(:,:) = e1v(:,:) * vn(:,:,jk) … … 142 263 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 143 264 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 144 zua = zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 145 zva =-zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 146 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 147 va(ji,jj,jk) = va(ji,jj,jk) + zva 265 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 266 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 148 267 END DO 149 268 END DO … … 152 271 ! ! =============== 153 272 154 ! save the relative & planetary vorticity trends for diagnostic 155 ! momentum trends 156 IF( l_trddyn ) THEN 157 ! Compute the planetary vorticity term trend 158 ! ! =============== 159 DO jk = 1, jpkm1 ! Horizontal slab 160 ! ! =============== 161 DO jj = 2, jpjm1 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 164 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 165 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 166 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 167 # if defined key_s_coord 168 zcu(ji,jj,jk) = zfact2 / e1u(ji,jj) * ( ff(ji ,jj-1) / fse3f(ji,jj-1,jk) * zy1 & 169 & + ff(ji ,jj ) / fse3f(ji,jj ,jk) * zy2 ) 170 zcv(ji,jj,jk) =-zfact2 / e2v(ji,jj) * ( ff(ji-1,jj ) / fse3f(ji-1,jj,jk) * zx1 & 171 & + ff(ji ,jj ) / fse3f(ji ,jj,jk) * zx2 ) 172 # else 173 zcu(ji,jj,jk) = zfact2 / e1u(ji,jj) * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 174 zcv(ji,jj,jk) =-zfact2 / e2v(ji,jj) * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 175 # endif 176 END DO 177 END DO 178 ! ! =============== 179 END DO ! End of slab 180 ! ! =============== 181 182 ! Compute the relative vorticity term trend 183 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 184 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 185 186 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 187 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 188 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 189 ENDIF 190 191 IF(ln_ctl) THEN ! print sum trends (used for debugging) 192 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 193 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 194 ENDIF 195 196 END SUBROUTINE dyn_vor_energy 197 198 199 SUBROUTINE dyn_vor_mixed( kt ) 200 !!---------------------------------------------------------------------- 201 !! *** ROUTINE dyn_vor_mixed *** 273 END SUBROUTINE vor_ene 274 275 276 SUBROUTINE vor_mix( kt ) 277 !!---------------------------------------------------------------------- 278 !! *** ROUTINE vor_mix *** 202 279 !! 203 280 !! ** Purpose : Compute the now total vorticity trend and add it to … … 209 286 !! ticity term and the horizontal kinetic energy for (f x uh), the 210 287 !! coriolis term. the now trend of the vorticity term is given by: 211 !! * s-coordinate (l k_sco=T), the e3. are inside the derivatives:288 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 212 289 !! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ] 213 290 !! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ] … … 234 311 !! 9.0 ! 04-08 (C. Talandier) New trends organization 235 312 !!---------------------------------------------------------------------- 236 !! * Modules used237 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace238 ztdva => sa ! use sa as 3D workspace239 313 !! * Arguments 240 314 INTEGER, INTENT( in ) :: kt ! ocean timestep index … … 247 321 REAL(wp), DIMENSION(jpi,jpj) :: & 248 322 zwx, zwy, zwz, zww ! temporary workspace 249 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &250 zcu, zcv ! " "251 323 !!---------------------------------------------------------------------- 252 324 253 325 IF( kt == nit000 ) THEN 254 326 IF(lwp) WRITE(numout,*) 255 IF(lwp) WRITE(numout,*) 'dyn _vor_mixed: vorticity term: mixed energy/enstrophy conserving scheme'256 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~'327 IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme' 328 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 257 329 ENDIF 258 330 … … 261 333 zfact2 = 0.5 * 0.5 262 334 263 ! Save ua and va trends 264 IF( l_trddyn ) THEN 265 ztdua(:,:,:) = ua(:,:,:) 266 ztdva(:,:,:) = va(:,:,:) 267 zcu(:,:,:) = 0.e0 268 zcv(:,:,:) = 0.e0 269 ENDIF 270 335 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 336 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 271 337 ! ! =============== 272 338 DO jk = 1, jpkm1 ! Horizontal slab … … 275 341 ! Relative and planetary potential vorticity and horizontal fluxes 276 342 ! ---------------------------------------------------------------- 277 IF( l k_sco ) THEN343 IF( ln_sco ) THEN 278 344 zwz(:,:) = ff (:,:) / fse3f(:,:,jk) 279 345 zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk) … … 310 376 ! ! =============== 311 377 312 ! save the relative & planetary vorticity trends for diagnostic 313 ! momentum trends 314 IF( l_trddyn ) THEN 315 ! Compute the planetary vorticity term trend 316 ! ! =============== 317 DO jk = 1, jpkm1 ! Horizontal slab 318 ! ! =============== 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 322 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 323 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 324 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 325 ! energy conserving formulation for planetary vorticity term 326 zcu(ji,jj,jk) = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 327 zcv(ji,jj,jk) =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 328 END DO 329 END DO 330 ! ! =============== 331 END DO ! End of slab 332 ! ! =============== 333 334 ! Compute the relative vorticity term trend 335 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 336 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 337 338 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 339 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 340 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 341 ENDIF 342 343 IF(ln_ctl) THEN ! print sum trends (used for debugging) 344 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 345 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 346 ENDIF 347 348 END SUBROUTINE dyn_vor_mixed 349 350 351 SUBROUTINE dyn_vor_enstrophy( kt ) 352 !!---------------------------------------------------------------------- 353 !! *** ROUTINE dyn_vor_enstrophy *** 378 END SUBROUTINE vor_mix 379 380 381 SUBROUTINE vor_ens( kt, cd_vor, pua, pva ) 382 !!---------------------------------------------------------------------- 383 !! *** ROUTINE vor_ens *** 354 384 !! 355 385 !! ** Purpose : Compute the now total vorticity trend and add it to … … 360 390 !! potential enstrophy of a horizontally non-divergent flow. the 361 391 !! trend of the vorticity term is given by: 362 !! * s-coordinate (l k_sco=T), the e3. are inside the derivative:392 !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: 363 393 !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 364 394 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] … … 381 411 !! 9.0 ! 04-08 (C. Talandier) New trends organization 382 412 !!---------------------------------------------------------------------- 383 !! * modules used384 USE oce, ONLY: zwx => ta, & ! use ta as 3D workspace385 zwy => sa ! use sa as 3D workspace386 413 !! * Arguments 387 414 INTEGER, INTENT( in ) :: kt ! ocean timestep 415 CHARACTER(len=3) , INTENT( in ) :: & 416 cd_vor ! define the vorticity considered 417 ! ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 418 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 419 pua, pva ! ???? 388 420 389 421 !! * Local declarations 390 422 INTEGER :: ji, jj, jk ! dummy loop indices 391 423 REAL(wp) :: & 392 zfact1, zua, zva, zuav, zvau ! temporary scalars 393 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 394 zcu, zcv, zwz, & ! temporary workspace 395 ztdua, ztdva ! temporary workspace 424 zfact1, zuav, zvau ! temporary scalars 425 REAL(wp), DIMENSION(jpi,jpj) :: & 426 zwx, zwy, zwz ! temporary workspace 396 427 !!---------------------------------------------------------------------- 397 428 398 429 IF( kt == nit000 ) THEN 399 430 IF(lwp) WRITE(numout,*) 400 IF(lwp) WRITE(numout,*) 'dyn _vor_enstrophy: vorticity term: enstrophy conserving scheme'401 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~~~~'431 IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 432 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 402 433 ENDIF 403 434 … … 405 436 zfact1 = 0.5 * 0.25 406 437 407 ! Save ua and va trends 408 IF( l_trddyn ) THEN 409 ztdua(:,:,:) = ua(:,:,:) 410 ztdva(:,:,:) = va(:,:,:) 411 zcu(:,:,:) = 0.e0 412 zcv(:,:,:) = 0.e0 413 ENDIF 414 438 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 439 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz ) 415 440 ! ! =============== 416 441 DO jk = 1, jpkm1 ! Horizontal slab … … 419 444 ! Potential vorticity and horizontal fluxes 420 445 ! ----------------------------------------- 421 IF( lk_sco ) THEN 446 SELECT CASE( cd_vor ) ! vorticity considered 447 CASE ( 'COR' ) ! planetary vorticcity (Coriolis) 448 zwz(:,:) = ff(:,:) 449 CASE ( 'VOR' ) ! relative vorticity 450 zwz(:,:) = rotn(:,:,jk) 451 CASE ( 'TOT' ) ! total vorticity (planetary + relative) 452 zwz(:,:) = rotn(:,:,jk) + ff(:,:) 453 END SELECT 454 455 IF( ln_sco ) THEN 422 456 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 423 457 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 424 zwz(ji,jj ,jk) = ( rotn(ji,jj,jk) + ff(ji,jj)) / fse3f(ji,jj,jk)425 zwx(ji,jj ,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)426 zwy(ji,jj ,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)458 zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 459 zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 460 zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 427 461 END DO 428 462 END DO … … 430 464 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 431 465 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 432 zwz(ji,jj,jk) = rotn(ji,jj,jk) + ff(ji,jj) 433 zwx(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 434 zwy(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 466 zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 467 zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 435 468 END DO 436 469 END DO … … 442 475 DO jj = 2, jpjm1 443 476 DO ji = fs_2, fs_jpim1 ! vector opt. 444 zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1,jk) + zwy(ji+1,jj-1,jk) & 445 + zwy(ji ,jj ,jk) + zwy(ji+1,jj ,jk) ) 446 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ,jk) + zwx(ji-1,jj+1,jk) & 447 + zwx(ji ,jj ,jk) + zwx(ji ,jj+1,jk) ) 448 449 zua = zuav * ( zwz(ji ,jj-1,jk) + zwz(ji,jj,jk) ) 450 zva = zvau * ( zwz(ji-1,jj ,jk) + zwz(ji,jj,jk) ) 451 452 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 453 va(ji,jj,jk) = va(ji,jj,jk) + zva 477 zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 478 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 479 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 480 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 481 482 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 483 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 454 484 END DO 455 485 END DO … … 458 488 ! ! =============== 459 489 460 461 ! save the relative & planetary vorticity trends for diagnostic 462 ! momentum trends 463 IF( l_trddyn ) THEN 464 ! Compute the planetary vorticity term trend 465 ! ! =============== 466 DO jk = 1, jpkm1 ! Horizontal slab 467 ! ! =============== 468 DO jj = 2, jpjm1 469 DO ji = fs_2, fs_jpim1 ! vector opt. 470 zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1,jk) + zwy(ji+1,jj-1,jk) & 471 & + zwy(ji ,jj ,jk) + zwy(ji+1,jj ,jk) ) 472 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ,jk) + zwx(ji-1,jj+1,jk) & 473 & + zwx(ji ,jj ,jk) + zwx(ji ,jj+1,jk) ) 474 # if defined key_s_coord 475 zcu(ji,jj,jk) = zuav * ( ff(ji ,jj-1) / fse3f(ji ,jj-1,jk) & 476 & + ff(ji ,jj ) / fse3f(ji ,jj ,jk) ) 477 zcv(ji,jj,jk) = zvau * ( ff(ji-1,jj ) / fse3f(ji-1,jj ,jk) & 478 & + ff(ji ,jj ) / fse3f(ji ,jj ,jk) ) 479 # else 480 zcu(ji,jj,jk) = zuav * ( ff(ji ,jj-1) + ff(ji,jj) ) 481 zcv(ji,jj,jk) = zvau * ( ff(ji-1,jj ) + ff(ji,jj) ) 482 # endif 483 END DO 484 END DO 485 ! ! =============== 486 END DO ! End of slab 487 ! ! =============== 488 489 ! Compute the relative vorticity term trend 490 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 491 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 492 493 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 494 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 495 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 496 ENDIF 497 498 IF(ln_ctl) THEN ! print sum trends (used for debugging) 499 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 500 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 501 ENDIF 502 503 END SUBROUTINE dyn_vor_enstrophy 504 505 506 SUBROUTINE dyn_vor_ene_ens( kt ) 507 !!---------------------------------------------------------------------- 508 !! *** ROUTINE dyn_vor_ene_ens *** 490 END SUBROUTINE vor_ens 491 492 493 SUBROUTINE vor_een( kt, cd_vor, pua, pva ) 494 !!---------------------------------------------------------------------- 495 !! *** ROUTINE vor_een *** 509 496 !! 510 497 !! ** Purpose : Compute the now total vorticity trend and add it to … … 516 503 !! when horizontal divergence is zero. 517 504 !! The trend of the vorticity term is given by: 518 !! * s-coordinate (l k_sco=T), the e3. are inside the derivatives:505 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 519 506 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 520 507 !! Add this trend to the general momentum trend (ua,va): … … 534 521 !! 9.0 ! 04-08 (C. Talandier) New trends organization 535 522 !!---------------------------------------------------------------------- 536 !! * Modules used537 USE oce, ONLY : ztdua => ta, & ! use ta as 3D workspace538 ztdva => sa ! use sa as 3D workspace539 540 523 !! * Arguments 541 524 INTEGER, INTENT( in ) :: kt ! ocean time-step index 525 CHARACTER(len=3) , INTENT( in ) :: & 526 cd_vor ! define the vorticity considered 527 ! ! ='COR' (planetary) ; ='VOR' (relative) ; ='TOT' (total) 528 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 529 pua, pva ! ???? 542 530 543 531 !! * Local declarations … … 547 535 REAL(wp), DIMENSION(jpi,jpj) :: & 548 536 zwx, zwy, zwz, & ! temporary workspace 549 ztnw, ztne, ztsw, ztse, & ! " " 550 zcor ! potential planetary vorticity (f/e3) 551 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 552 zcu, zcv ! temporary workspace 537 ztnw, ztne, ztsw, ztse ! " " 553 538 REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: & 554 539 ze3f … … 557 542 IF( kt == nit000 ) THEN 558 543 IF(lwp) WRITE(numout,*) 559 IF(lwp) WRITE(numout,*) 'dyn _vor_ene_ens: vorticity term: energy and enstrophy conserving scheme'560 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ~~~~'544 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 545 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 561 546 562 547 DO jk = 1, jpk … … 565 550 ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 566 551 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25 567 !!! ze3f(ji,jj,jk) = MAX( ze3f(ji,jj,jk) , 1.e-20)568 552 IF( ze3f(ji,jj,jk) /= 0.e0 ) ze3f(ji,jj,jk) = 1.e0 / ze3f(ji,jj,jk) 569 553 END DO … … 576 560 zfac12 = 1.e0 / 12.e0 577 561 578 ! Save ua and va trends579 IF( l_trddyn ) THEN580 ztdua(:,:,:) = ua(:,:,:)581 ztdva(:,:,:) = va(:,:,:)582 zcu(:,:,:) = 0.e0583 zcv(:,:,:) = 0.e0584 ENDIF585 562 563 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 564 !$OMP PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 586 565 ! ! =============== 587 566 DO jk = 1, jpkm1 ! Horizontal slab … … 590 569 ! Potential vorticity and horizontal fluxes 591 570 ! ----------------------------------------- 592 !!!bug zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 593 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 571 SELECT CASE( cd_vor ) ! vorticity considered 572 CASE ( 'COR' ) ! planetary vorticcity (Coriolis) 573 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 574 CASE ( 'VOR' ) ! relative vorticity 575 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 576 CASE ( 'TOT' ) ! total vorticity (planetary + relative) 577 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 578 END SELECT 579 594 580 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 595 581 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 596 zcor(:,:) = ff(:,:) * ze3f(:,:,jk)597 582 598 583 ! Compute and add the vorticity term trend … … 621 606 zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 622 607 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 623 ua(ji,jj,jk) =ua(ji,jj,jk) + zua624 va(ji,jj,jk) =va(ji,jj,jk) + zva608 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 609 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 625 610 END DO 626 611 END DO … … 629 614 ! ! =============== 630 615 631 ! save the relative & planetary vorticity trends for diagnostic 632 ! momentum trends 633 IF( l_trddyn ) THEN 634 ! Compute the planetary vorticity term trend 635 ! ! =============== 636 DO jk = 1, jpkm1 ! Horizontal slab 637 ! ! =============== 638 DO jj = 2, jpjm1 639 DO ji = fs_2, fs_jpim1 ! vector opt. 640 zcu(ji,jj,jk) = + zfac12 / e1u(ji,jj) * ( zcor(ji,jj ) * zwy(ji ,jj ) + zcor(ji+1,jj) * zwy(ji+1,jj ) & 641 & + zcor(ji,jj ) * zwy(ji ,jj-1) + zcor(ji+1,jj) * zwy(ji+1,jj-1) ) 642 zcv(ji,jj,jk) = - zfac12 / e2v(ji,jj) * ( zcor(ji,jj+1) * zwx(ji-1,jj+1) + zcor(ji,jj+1) * zwx(ji ,jj+1) & 643 & + zcor(ji,jj ) * zwx(ji-1,jj ) + zcor(ji,jj ) * zwx(ji ,jj ) ) 644 END DO 645 END DO 646 ! ! =============== 647 END DO ! End of slab 648 ! ! =============== 649 650 ! Compute the relative vorticity term trend 651 ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:) - zcu(:,:,:) 652 ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:) - zcv(:,:,:) 653 654 CALL trd_mod(zcu , zcv , jpdtdpvo, 'DYN', kt) 655 CALL trd_mod(zcu , zcv , jpdtddat, 'DYN', kt) 656 CALL trd_mod(ztdua, ztdva, jpdtdrvo, 'DYN', kt) 657 ENDIF 658 659 IF(ln_ctl) THEN ! print sum trends (used for debugging) 660 CALL prt_ctl(tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 661 & tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn') 662 ENDIF 663 664 END SUBROUTINE dyn_vor_ene_ens 665 666 667 SUBROUTINE dyn_vor_ctl 616 END SUBROUTINE vor_een 617 618 619 SUBROUTINE vor_ctl 668 620 !!--------------------------------------------------------------------- 669 !! *** ROUTINE dyn_vor_ctl ***621 !! *** ROUTINE vor_ctl *** 670 622 !! 671 623 !! ** Purpose : Control the consistency between cpp options for … … 691 643 IF(lwp) THEN 692 644 WRITE(numout,*) 693 WRITE(numout,*) 'dyn _vor_ctl : vorticity term : read namelist and control the consistency'645 WRITE(numout,*) 'dyn:vor_ctl : vorticity term : read namelist and control the consistency' 694 646 WRITE(numout,*) '~~~~~~~~~~~' 695 647 WRITE(numout,*) ' Namelist nam_dynvor : oice of the vorticity term scheme' 648 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 696 649 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 697 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene698 650 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 699 651 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een … … 701 653 702 654 ioptio = 0 703 655 IF( ln_dynvor_ene ) THEN 656 nvor = 0 657 IF(lwp) WRITE(numout,*) 658 IF(lwp) WRITE(numout,*) ' vorticity term : energy conserving scheme' 659 ioptio = ioptio + 1 660 ENDIF 704 661 IF( ln_dynvor_ens ) THEN 662 nvor = 1 705 663 IF(lwp) WRITE(numout,*) 706 664 IF(lwp) WRITE(numout,*) ' vorticity term : enstrophy conserving scheme' 707 665 ioptio = ioptio + 1 708 666 ENDIF 709 IF( ln_dynvor_ene ) THEN710 IF(lwp) WRITE(numout,*)711 IF(lwp) WRITE(numout,*) ' vorticity term : energy conserving scheme'712 ioptio = ioptio + 1713 ENDIF714 667 IF( ln_dynvor_mix ) THEN 668 nvor = 2 715 669 IF(lwp) WRITE(numout,*) 716 670 IF(lwp) WRITE(numout,*) ' vorticity term : mixed enstrophy/energy conserving scheme' … … 718 672 ENDIF 719 673 IF( ln_dynvor_een ) THEN 674 nvor = 3 720 675 IF(lwp) WRITE(numout,*) 721 676 IF(lwp) WRITE(numout,*) ' vorticity term : energy and enstrophy conserving scheme' 722 677 ioptio = ioptio + 1 723 678 ENDIF 724 IF ( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN 725 WRITE(numout,cform_err) 726 IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme' 727 nstop = nstop + 1 728 ENDIF 729 730 END SUBROUTINE dyn_vor_ctl 679 IF( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN 680 WRITE(numout,cform_err) 681 IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme' 682 nstop = nstop + 1 683 ENDIF 684 IF( lk_esopa ) THEN 685 nvor = -1 686 IF(lwp ) WRITE(numout,*) ' esopa test: use all lateral physics options' 687 ENDIF 688 IF(lwp) WRITE(numout,*) ' choice of vor_... nvor = ', nvor 689 690 END SUBROUTINE vor_ctl 731 691 732 692 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
