Changeset 791 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90
- Timestamp:
- 2008-01-12T21:33:34+01:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r786 r791 5 5 !!====================================================================== 6 6 !! History : 1.0 ! 06-09 (G. Reffray) Original code 7 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !!---------------------------------------------------------------------- 9 9 … … 30 30 PUBLIC tra_adv_qck ! routine called by traadv.F90 31 31 32 REAL(wp), DIMENSION(jpi,jpj) :: zbtr233 32 REAL(wp), DIMENSION(jpi,jpj,jpk) :: sl 34 REAL(wp) :: cst1, cst2, dt, coef1 ! temporary scalars35 INTEGER :: ji, jj, jk ! dummy loop indices36 33 37 34 !! * Substitutions … … 40 37 !!---------------------------------------------------------------------- 41 38 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 42 !! $Id :$39 !! $Id$ 43 40 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 41 !!---------------------------------------------------------------------- … … 94 91 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 95 92 !! 96 REAL(wp) :: z2 ! temporary scalar 93 INTEGER :: ji, jj, jk ! dummy loop indices 94 REAL(wp) :: zdt, z2, zcoef1 ! temporary scalars 95 REAL(wp) :: zbtr ! temporary scalars 97 96 !!---------------------------------------------------------------------- 98 97 … … 101 100 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme' 102 101 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 103 104 zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )105 cst1 = 1./12.106 cst2 = 2./3.107 102 ENDIF 108 103 109 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 110 ELSE ; z2 = 2. 104 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ! euler time-stepping 105 ELSE ; z2 = 2. ! leap-frog time-stepping 111 106 ENDIF 112 107 … … 116 111 !--------------------------------------------------------------- 117 112 113 !!gm : optimisation: sl compute twice (for t & s, and even more for trc) 114 !!gm note that sl is in permanent memory be used as workspace in the vertical part ! 118 115 sl(:,:,:) = 100. 119 116 120 ! =============== 121 DO jk = 1, jpkm1 ! Horizontal slab 122 ! ! =============== 123 dt = z2 * rdttra(jk) 124 DO jj = 2, jpjm1 125 DO ji = 2, fs_jpim1 ! vector opt. 126 coef1 = 1.e-12 127 IF (pun(ji-1,jj ,jk ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj ,jk ))*dt/e1t(ji,jj) 128 IF (pun(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pun(ji ,jj ,jk ))*dt/e1t(ji,jj) 129 IF (pvn(ji ,jj-1,jk ).LT.0.) coef1 = coef1 + ABS(pvn(ji ,jj-1,jk ))*dt/e2t(ji,jj) 130 IF (pvn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pvn(ji ,jj ,jk ))*dt/e2t(ji,jj) 131 IF (pwn(ji ,jj ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk+1))*dt/(fse3t(ji,jj,jk)) 132 IF (pwn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk ))*dt/(fse3t(ji,jj,jk)) 133 sl(ji,jj,jk) = 1./coef1 134 sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk)) 135 sl(ji,jj,jk) = MAX(1. ,sl(ji,jj,jk)) 136 ENDDO 137 ENDDO 138 ENDDO 139 !--- Lateral boundary conditions on the limiter slope 140 CALL lbc_lnk( sl(:,:,:), 'T', 1. ) 117 DO jk = 1, jpkm1 118 zdt = z2 * rdttra(jk) 119 DO jj = 2, jpjm1 120 DO ji = 2, fs_jpim1 ! vector opt. 121 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) *fse3t(ji,jj,jk) ) 122 zcoef1 = 1.e-12 123 IF( pun(ji-1,jj ,jk ) < 0.e0 ) zcoef1 = zcoef1 + ABS( pun(ji-1,jj ,jk ) ) * zdt * zbtr 124 IF( pun(ji ,jj ,jk ) > 0.e0 ) zcoef1 = zcoef1 + ABS( pun(ji ,jj ,jk ) ) * zdt * zbtr 125 IF( pvn(ji ,jj-1,jk ) < 0.e0 ) zcoef1 = zcoef1 + ABS( pvn(ji ,jj-1,jk ) ) * zdt * zbtr 126 IF( pvn(ji ,jj ,jk ) > 0.e0 ) zcoef1 = zcoef1 + ABS( pvn(ji ,jj ,jk ) ) * zdt * zbtr 127 IF( pwn(ji ,jj ,jk+1) < 0.e0 ) zcoef1 = zcoef1 + ABS( pwn(ji ,jj ,jk+1) ) * zdt * zbtr 128 IF( pwn(ji ,jj ,jk ) > 0.e0 ) zcoef1 = zcoef1 + ABS( pwn(ji ,jj ,jk ) ) * zdt * zbtr 129 sl(ji,jj,jk) = 1. / zcoef1 130 sl(ji,jj,jk) = MIN( 100., sl(ji,jj,jk) ) 131 sl(ji,jj,jk) = MAX( 1. , sl(ji,jj,jk) ) 132 END DO 133 END DO 134 END DO 135 CALL lbc_lnk( sl(:,:,:), 'T', 1. ) ! Lateral boundary conditions on the limiter slope 141 136 142 137 … … 146 141 CALL tra_adv_qck_hor( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) 147 142 143 ! ! control print 148 144 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - had: ', mask1=tmask, clinfo3=cdtype ) 149 145 … … 152 148 !------------------------------------------------------------------------- 153 149 154 CALL tra_adv_qck_ver( pwn, ptn , pta, z2 ) 155 150 CALL tra_adv_qck_ver( pwn, ptn , pta ) 151 152 ! ! control print 156 153 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - zad: ', mask1=tmask, clinfo3=cdtype ) 157 154 ! … … 159 156 160 157 161 SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) 162 !!---------------------------------------------------------------------- 158 SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, pz2 ) 159 !!---------------------------------------------------------------------- 160 !! *** ROUTINE tra_adv_qck_hor *** 161 !! 162 !! ** Purpose : 163 163 !! 164 164 !!---------------------------------------------------------------------- … … 166 166 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 167 167 INTEGER , INTENT(in ) :: ktra ! tracer index 168 REAL(wp) , INTENT(in ) :: z2 ! ???168 REAL(wp) , INTENT(in ) :: pz2 ! =1 or 2 (euler or leap-frog) 169 169 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! horizontal effective velocity 170 170 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptb ! before tracer field 171 171 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 172 173 REAL(wp) :: za, zbtr, e1, e2, c, dir, fu, fc, fd ! temporary scalars 174 REAL(wp) :: coef2, coef3, fho, mask, dx 175 REAL(wp), DIMENSION(jpi,jpj) :: zee 176 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmask, zlap, dwst, lim 172 !! 173 INTEGER :: ji, jj, jk ! dummy loop indices 174 REAL(wp) :: zdt ! temporary scalars 175 REAL(wp) :: zbtr, zc, zdir, zfu, zfc, zfd ! temporary scalars 176 REAL(wp) :: zcoef1, zcoef2, zcoef3, zfho, zmsk, zdx 177 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmask, zlap, zdwst, zlim 177 178 !!---------------------------------------------------------------------- 178 179 … … 182 183 zmask(:,:,:)= 0.e0 183 184 zlap (:,:,:)= 0.e0 184 dwst(:,:,:)= 0.e0185 lim(:,:,:)= 0.e0185 zdwst(:,:,:)= 0.e0 186 zlim (:,:,:)= 0.e0 186 187 187 188 !---------------------------------------------------------------------- … … 189 190 !---------------------------------------------------------------------- 190 191 191 ! =============== 192 DO jk = 1, jpkm1 ! Horizontal slab 193 ! ! =============== 194 ! Initialization of metric arrays (for z- or s-coordinates) 195 ! --------------------------------------------------------- 196 DO jj = 1, jpjm1 197 DO ji = 1, fs_jpim1 ! vector opt. 198 #if defined key_zco 199 ! z-coordinates, no vertical scale factors 200 zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) 201 #else 202 ! vertical scale factor are used 203 zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 204 #endif 205 END DO 206 END DO 192 DO jk = 1, jpkm1 207 193 208 194 ! Laplacian of tracers (at before time step) … … 211 197 DO jj = 1, jpjm1 212 198 DO ji = 1, fs_jpim1 ! vector opt. 213 zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji+1,jj ,jk) - ptb(ji,jj,jk) ) 214 END DO 215 END DO 216 DO jj = 2, jpjm1 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 #if defined key_zco 219 zee(ji,jj) = e1t(ji,jj) / e2t(ji,jj) 220 #else 221 zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk)) 222 #endif 223 zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji-1,jj,jk) ) 199 zmask(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) & 200 & * ( ptb(ji+1,jj ,jk) - ptb(ji,jj,jk) ) / e1u(ji,jj) 201 END DO 202 END DO 203 DO jj = 2, jpjm1 204 DO ji = fs_2, fs_jpim1 ! vector opt. 205 zlap(ji,jj,jk) = e1t(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji-1,jj,jk) ) & 206 & / ( e2t(ji,jj) * fse3t(ji,jj,jk) ) 224 207 END DO 225 208 END DO … … 231 214 DO ji = 2, fs_jpim1 ! vector opt. 232 215 ! Upstream in the x-direction for the tracer 233 zmask(ji,jj,jk) = ptb(ji-1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji-1,jj,jk))216 zmask(ji,jj,jk) = ptb(ji-1,jj,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji-1,jj,jk) ) 234 217 ! Downstream in the x-direction for the tracer 235 dwst (ji,jj,jk) =ptb(ji+1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji+1,jj,jk)) 236 ENDDO 237 ENDDO 238 END DO 239 !--- Lateral boundary conditions on the laplacian (unchanged sgn) 240 CALL lbc_lnk( zlap(:,:,:), 'T', 1. ) 241 !--- Lateral boundary conditions for the lim function 242 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. ) 243 ! =============== 244 DO jk = 1, jpkm1 ! Horizontal slab 245 ! ! =============== 246 ! --- lim at the U-points in function of the direction of the flow 247 ! ---------------------------------------------------------------- 218 zdwst(ji,jj,jk) = ptb(ji+1,jj,jk) + sl(ji,jj,jk) * ( ptb(ji,jj,jk) - ptb(ji+1,jj,jk) ) 219 END DO 220 END DO 221 END DO 222 223 !--- Lateral boundary conditions on the laplacian and lim functions (unchanged sgn) 224 CALL lbc_lnk( zlap (:,:,:), 'T', 1. ) 225 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zdwst(:,:,:), 'T', 1. ) 226 227 ! --- lim at the U-points in function of the direction of the flow 228 ! ---------------------------------------------------------------- 229 DO jk = 1, jpkm1 248 230 DO jj = 1, jpjm1 249 231 DO ji = 1, fs_jpim1 ! vector opt. 250 dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0251 lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk)232 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun>0 : diru = 1 otherwise diru = 0 233 zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji+1,jj,jk) 252 234 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 253 zmask(ji,jj,jk) =tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2254 END DO 255 END DO 256 END DO 257 !--- Lateral boundary conditions for the mask258 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) 235 zmask(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2 236 END DO 237 END DO 238 END DO 239 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) !--- Lateral boundary conditions for the mask 240 259 241 260 242 ! Horizontal advective fluxes 261 243 ! --------------------------- 262 ! =============== 263 DO jk = 1, jpkm1 ! Horizontal slab 264 ! =============== 265 dt = z2 * rdttra(jk) 244 DO jk = 1, jpkm1 245 zdt = pz2 * rdttra(jk) 266 246 !--- tracer flux at u and v-points 267 247 DO jj = 1, jpjm1 268 248 DO ji = 1, fs_jpim1 ! vector opt. 269 #if defined key_zco 270 e2 = e2u(ji,jj) 271 #else 272 e2 = e2u(ji,jj) * fse3u(ji,jj,jk) 273 #endif 274 dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0 275 276 dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj) 277 c = ABS(pun(ji,jj,jk))*dt/dx ! (0<cx<1 : Courant number on x-direction) 278 279 fu = lim(ji,jj,jk) ! FU + sl(FC-FU) in the x-direction for T 280 fc = dir*ptb(ji ,jj,jk)+(1-dir)*ptb(ji+1,jj,jk) ! FC in the x-direction for T 281 fd = dir*ptb(ji+1,jj,jk)+(1-dir)*ptb(ji ,jj,jk) ! FD in the x-direction for T 249 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun>0 : diru = 1 otherwise diru = 0 250 zdx = ( zdir * e1t(ji,jj) + (1-zdir)* e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk) 251 zc = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<cx<1 : Courant number on x-direction) 252 253 zfu = zlim(ji,jj,jk) ! FU + sl(FC-FU) in the x-direction for T 254 zfc = zdir * ptb(ji ,jj,jk) + (1-zdir) * ptb(ji+1,jj,jk) ! FC in the x-direction for T 255 zfd = zdir * ptb(ji+1,jj,jk) + (1-zdir) * ptb(ji ,jj,jk) ! FD in the x-direction for T 282 256 283 257 !--- QUICKEST scheme 284 258 ! Temperature on the x-direction 285 coef1 = 0.5*(fc+fd)286 coef2 = 0.5*c*(fd-fc)287 coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji+1,jj,jk) )288 fho = coef1-coef2-coef3289 fho = bound(fu,fd,fc,fho)259 zcoef1 = 0.5 * ( zfc + zfd ) 260 zcoef2 = 0.5 * zc * ( zfd - zfc ) 261 zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji+1,jj,jk) ) 262 zfho = zcoef1 - zcoef2 - zcoef3 263 zfho = bound( zfu, zfd, zfc, zfho ) 290 264 !--- If the second ustream point is a land point 291 265 !--- the flux is computed by the 1st order UPWIND scheme 292 mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji+1,jj,jk)293 fho = mask*fho + (1-mask)*fc294 dwst(ji,jj,jk)=e2*pun(ji,jj,jk)*fho266 zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji+1,jj,jk) 267 zfho = zmsk * zfho + (1-zmsk) * zfc 268 zdwst(ji,jj,jk) = pun(ji,jj,jk) * zfho 295 269 END DO 296 270 END DO … … 300 274 DO ji = fs_2, fs_jpim1 ! vector opt. 301 275 !--- horizontal advective trends 302 #if defined key_zco 303 zbtr = zbtr2(ji,jj) 304 #else 305 zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) 306 #endif 307 za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj ,jk) ) 276 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 308 277 !--- add it to the general tracer trends 309 pta(ji,jj,jk) = pta(ji,jj,jk) + za 310 END DO 311 END DO 312 ! ! =============== 313 END DO ! End of slab 314 ! ! =============== 315 316 ! Save the horizontal advective trends for diagnostic 317 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, dwst, pun, ptb ) 278 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zdwst(ji,jj,jk) - zdwst(ji-1,jj ,jk) ) 279 END DO 280 END DO 281 END DO 282 283 ! !-- trend diagnostic 284 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zdwst, pun, ptb ) 318 285 319 286 … … 321 288 ! I. Part 2 : y-direction 322 289 !---------------------------------------------------------------------- 323 ! ============== 324 DO jk = 1, jpkm1 ! Horizontal slab 325 ! ! =============== 326 ! Initialization of metric arrays (for z- or s-coordinates) 327 ! --------------------------------------------------------- 328 DO jj = 1, jpjm1 329 DO ji = 1, fs_jpim1 ! vector opt. 330 #if defined key_zco 331 ! z-coordinates, no vertical scale factors 332 zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) 333 #else 334 ! s-coordinates, vertical scale factor are used 335 zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 336 #endif 337 END DO 338 END DO 290 291 DO jk = 1, jpkm1 339 292 340 293 ! Laplacian of tracers (at before time step) … … 343 296 DO jj = 1, jpjm1 344 297 DO ji = 1, fs_jpim1 ! vector opt. 345 zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji ,jj+1,jk) - ptb(ji,jj,jk) ) 298 zmask(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) & 299 & * ( ptb(ji ,jj+1,jk) - ptb(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk) 346 300 END DO 347 301 END DO … … 349 303 DO jj = 2, jpjm1 350 304 DO ji = fs_2, fs_jpim1 ! vector opt. 351 #if defined key_zco 352 zee(ji,jj) = e2t(ji,jj) / e1t(ji,jj) 353 #else 354 zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk)) 355 #endif 356 zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji,jj-1,jk) ) 305 zlap(ji,jj,jk) = e2t(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji,jj-1,jk) ) & 306 & / ( e1t(ji,jj) * fse3t(ji,jj,jk) ) 357 307 END DO 358 308 END DO … … 362 312 DO ji = 2, fs_jpim1 ! vector opt. 363 313 ! Upstream in the y-direction for the tracer 364 zmask(ji,jj,jk) =ptb(ji,jj-1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj-1,jk))314 zmask(ji,jj,jk) = ptb(ji,jj-1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj-1,jk) ) 365 315 ! Downstream in the y-direction for the tracer 366 dwst (ji,jj,jk)=ptb(ji,jj+1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj+1,jk))316 zdwst(ji,jj,jk) = ptb(ji,jj+1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj+1,jk) ) 367 317 ENDDO 368 318 ENDDO 369 319 END DO 370 !--- Lateral boundary conditions on the laplacian (unchanged sgn) 371 CALL lbc_lnk( zlap(:,:,:), 'T', 1. ) 372 !--- Lateral boundary conditions for the lim function 373 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. ) 374 375 DO jk = 1, jpkm1 ! Horizontal slab 376 ! ! =============== 320 !--- Lateral boundary conditions on the laplacian and lim function (unchanged sgn) 321 CALL lbc_lnk( zlap (:,:,:), 'T', 1. ) 322 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zdwst(:,:,:), 'T', 1. ) 323 324 DO jk = 1, jpkm1 325 377 326 ! --- lim at the V-points in function of the direction of the flow 378 327 ! ---------------------------------------------------------------- 379 328 DO jj = 1, jpjm1 380 329 DO ji = 1, fs_jpim1 ! vector opt. 381 dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0382 lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk)330 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pvn>0 : dirv = 1 otherwise dirv = 0 331 zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji,jj+1,jk) 383 332 ! Mask at the T-points in the y-direction (mask=0 or mask=1) 384 zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 385 END DO 386 END DO 387 END DO 388 !--- Lateral boundary conditions for the mask 389 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) 333 zmask(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 334 END DO 335 END DO 336 END DO 337 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) !--- Lateral boundary conditions for the mask 390 338 391 339 ! Horizontal advective fluxes 392 ! ------------------------------- 393 ! =============== 394 DO jk = 1, jpkm1 ! Horizontal slab 395 ! =============== 396 dt = z2 * rdttra(jk) 340 ! --------------------------- 341 342 DO jk = 1, jpkm1 343 zdt = pz2 * rdttra(jk) 397 344 !--- tracer flux at u and v-points 398 345 DO jj = 1, jpjm1 399 346 DO ji = 1, fs_jpim1 ! vector opt. 400 #if defined key_zco 401 e1 = e1v(ji,jj) 402 #else 403 e1 = e1v(ji,jj) * fse3v(ji,jj,jk) 404 #endif 405 dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0 406 407 dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1) 408 c = ABS(pvn(ji,jj,jk))*dt/dx ! (0<cy<1 : Courant number on y-direction) 409 410 fu = lim(ji,jj,jk) ! FU + sl(FC-FU) in the y-direction for T 411 fc = dir*ptb(ji,jj ,jk)+(1-dir)*ptb(ji,jj+1,jk) ! FC in the y-direction for T 412 fd = dir*ptb(ji,jj+1,jk)+(1-dir)*ptb(ji,jj ,jk) ! FD in the y-direction for T 347 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pvn>0 : dirv = 1 otherwise dirv = 0 348 zdx = ( zdir * e2t(ji,jj) + (1-zdir)* e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk) 349 zc = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<cy<1 : Courant number on y-direction) 350 351 zfu = zlim(ji,jj,jk) ! FU + sl(FC-FU) in the y-direction for T 352 zfc = zdir * ptb(ji,jj ,jk) + (1-zdir) * ptb(ji,jj+1,jk) ! FC in the y-direction for T 353 zfd = zdir * ptb(ji,jj+1,jk) + (1-zdir) * ptb(ji,jj ,jk) ! FD in the y-direction for T 413 354 414 355 !--- QUICKEST scheme 415 356 ! Temperature on the y-direction 416 coef1 = 0.5*(fc+fd)417 coef2 = 0.5*c*(fd-fc)418 coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji,jj+1,jk) )419 fho = coef1-coef2-coef3420 fho = bound(fu,fd,fc,fho)357 zcoef1 = 0.5 * ( zfc + zfd ) 358 zcoef2 = 0.5 * zc * ( zfd - zfc ) 359 zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji,jj+1,jk) ) 360 zfho = zcoef1 - zcoef2 - zcoef3 361 zfho = bound( zfu, zfd, zfc, zfho ) 421 362 !--- If the second ustream point is a land point 422 363 !--- the flux is computed by the 1st order UPWIND scheme 423 mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji,jj+1,jk)424 fho = mask*fho + (1-mask)*fc425 dwst(ji,jj,jk)=e1*pvn(ji,jj,jk)*fho364 zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji,jj+1,jk) 365 zfho = zmsk * zfho + (1-zmsk) * zfc 366 zdwst(ji,jj,jk) = pvn(ji,jj,jk) * zfho 426 367 END DO 427 368 END DO … … 430 371 DO jj = 2, jpjm1 431 372 DO ji = fs_2, fs_jpim1 ! vector opt. 432 !--- horizontal advective trends 433 #if defined key_zco 434 zbtr = zbtr2(ji,jj) 435 #else 436 zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) 437 #endif 438 za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji ,jj-1,jk) ) 439 !--- add it to the general tracer trends 440 pta(ji,jj,jk) = pta(ji,jj,jk) + za 441 END DO 442 END DO 443 ! ! =============== 444 END DO ! End of slab 445 ! ! =============== 446 447 ! Save the horizontal advective trends for diagnostic 448 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, dwst, pvn, ptb ) 449 450 ! "zonal" mean advective heat and salt transport 373 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 374 ! horizontal advective trends added to the general tracer trends 375 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zdwst(ji,jj,jk) - zdwst(ji,jj-1,jk) ) 376 END DO 377 END DO 378 END DO 379 380 ! !-- trend diagnostic 381 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zdwst, pvn, ptb ) 382 ! ! "Poleward" heat or salt transport 451 383 IF( ln_diaptr .AND. cdtype == 'TRA' .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 452 #if defined key_zco 453 DO jk = 1, jpkm1 454 DO jj = 2, jpjm1 455 DO ji = fs_2, fs_jpim1 ! vector opt. 456 dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk) 457 END DO 458 END DO 459 END DO 460 # endif 461 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( dwst(:,:,:) ) 462 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( dwst(:,:,:) ) 384 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zdwst(:,:,:) ) 385 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zdwst(:,:,:) ) 463 386 ENDIF 464 387 ! … … 466 389 467 390 468 SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta, z2 ) 469 !!---------------------------------------------------------------------- 470 !! 471 !!---------------------------------------------------------------------- 472 REAL(wp), INTENT(in ) :: z2 473 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn 474 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn 475 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta 476 !! 477 REAL(wp) :: za, ze3tr, dt, dir, fc, fd ! temporary scalars 391 SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta ) 392 !!---------------------------------------------------------------------- 393 !! *** ROUTINE tra_adv_qck_ver *** 394 !! 395 !! ** Purpose : 396 !! 397 !!---------------------------------------------------------------------- 398 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical transport 399 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now tracer 400 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 401 !! 402 INTEGER :: ji, jj , jk ! dummy loop indices 403 REAL(wp) :: zbtr, zdir, zfc, zfd ! temporary scalars 404 !!---------------------------------------------------------------------- 478 405 479 406 ! Vertical advection … … 483 410 ! ---------------------------- 484 411 485 sl(:,:,jpk) = 0.e0 !Bottom value : flux set to zero 486 487 ! Surface value 488 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid : flux set to zero 489 sl(:,:, 1 ) = 0.e0 490 ELSE ! free surface-constant volume 491 sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) 412 sl(:,:,jpk) = 0.e0 ! bottom value 413 414 ! ! Surface value 415 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ; sl(:,:, 1 ) = 0.e0 ! rigid lid or non-linear fs 416 ELSE ; sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! linear free surface 492 417 ENDIF 418 419 !!gm bug: code au true 2nd order scheme 420 !!gm : sl used as work array: not good idea for optimisation (sl compute once for all tracers...) 493 421 494 422 ! Second order centered tracer flux at w-point 495 423 DO jk = 2, jpkm1 496 dt = z2 * rdttra(jk) 497 DO jj = 2, jpjm1 498 DO ji = fs_2, fs_jpim1 ! vector opt. 499 dir = 0.5 + sign(0.5,pwn(ji,jj,jk)) ! if pwn>0 : dirw = 1 otherwise dirw = 0 500 fc = dir*ptn(ji,jj,jk )*fse3t(ji,jj,jk-1)+(1-dir)*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk ) ! FC in the z-direction for T 501 fd = dir*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk )+(1-dir)*ptn(ji,jj,jk )*fse3t(ji,jj,jk-1) ! FD in the z-direction for T 424 DO jj = 2, jpjm1 425 DO ji = fs_2, fs_jpim1 ! vector opt. 426 zdir = 0.5 + SIGN( 0.5, pwn(ji,jj,jk) ) ! if pwn>0 : dirw = 1 otherwise dirw = 0 427 zfc = zdir * ptn(ji,jj,jk ) * fse3t(ji,jj,jk-1) + (1-zdir) * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk ) ! FC 428 zfd = zdir * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk ) + (1-zdir) * ptn(ji,jj,jk ) * fse3t(ji,jj,jk-1) ! FD 502 429 !--- Second order centered scheme 503 sl(ji,jj,jk) =pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk))430 sl(ji,jj,jk) = pwn(ji,jj,jk) * ( zfc + zfd ) / ( fse3t(ji,jj,jk-1) + fse3t(ji,jj,jk) ) 504 431 END DO 505 432 END DO … … 511 438 DO jj = 2, jpjm1 512 439 DO ji = fs_2, fs_jpim1 ! vector opt. 513 ze3tr = 1. / fse3t(ji,jj,jk) 514 ! vertical advective trends 515 za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) 516 ! add it to the general tracer trends 517 pta(ji,jj,jk) = pta(ji,jj,jk) + za 440 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 441 ! vertical advective trends added to the general tracer trends 442 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) 518 443 END DO 519 444 END DO … … 523 448 524 449 525 REAL FUNCTION bound(fu,fd,fc,fho) 526 REAL(wp) :: fu, fd, fc, fho, fref1, fref2 527 fref1 = fu 528 fref2 = MAX( MIN( fc , fd ), MIN( MAX( fc , fd ), fref1 ) ) 529 bound = MAX( MIN( fho, fc ), MIN( MAX( fho, fc ), fref2 ) ) 450 FUNCTION bound( pfu, pfd, pfc, pfho ) RESULT( pbound ) 451 !!---------------------------------------------------------------------- 452 !! *** FUNCTION bound *** 453 !! 454 !! ** Purpose : ??? 455 !!---------------------------------------------------------------------- 456 REAL(wp), INTENT(in) :: pfu, pfd, pfc, pfho 457 REAL(wp) :: pbound 458 !! 459 REAL(wp) :: zfref1, zfref2 460 !!---------------------------------------------------------------------- 461 zfref1 = pfu 462 zfref2 = MAX( MIN( pfc , pfd ), MIN( MAX( pfc , pfd ), zfref1 ) ) 463 pbound = MAX( MIN( pfho, pfc ), MIN( MAX( pfho, pfc ), zfref2 ) ) 464 ! 530 465 END FUNCTION 531 466
Note: See TracChangeset
for help on using the changeset viewer.
