| 1 | MODULE traadv_qck |
|---|
| 2 | !!====================================================================== |
|---|
| 3 | !! *** MODULE traadv_qck *** |
|---|
| 4 | !! Ocean tracers: horizontal & vertical advective trend using QUICKEST scheme |
|---|
| 5 | !!====================================================================== |
|---|
| 6 | !! History : 1.0 ! 06-09 (G. Reffray) Original code |
|---|
| 7 | !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport |
|---|
| 8 | !!---------------------------------------------------------------------- |
|---|
| 9 | |
|---|
| 10 | !!---------------------------------------------------------------------- |
|---|
| 11 | !! tra_adv_qck : update the tracer trend with the horizontal advection |
|---|
| 12 | !! trends using a 3rd order finite difference scheme ????? |
|---|
| 13 | !! The vertical advection scheme is the 2nd centered scheme ????? |
|---|
| 14 | !!---------------------------------------------------------------------- |
|---|
| 15 | USE dom_oce ! ocean space and time domain |
|---|
| 16 | USE trdmod ! ocean active tracers trends |
|---|
| 17 | USE trdmod_oce ! ocean variables trends |
|---|
| 18 | USE flxrnf ! |
|---|
| 19 | USE ocfzpt ! |
|---|
| 20 | USE lib_mpp |
|---|
| 21 | USE lbclnk ! ocean lateral boundary condition (or mpp link) |
|---|
| 22 | USE in_out_manager ! I/O manager |
|---|
| 23 | USE diaptr ! poleward transport diagnostics |
|---|
| 24 | USE dynspg_oce ! |
|---|
| 25 | USE prtctl ! Print control |
|---|
| 26 | |
|---|
| 27 | IMPLICIT NONE |
|---|
| 28 | PRIVATE |
|---|
| 29 | |
|---|
| 30 | PUBLIC tra_adv_qck ! routine called by traadv.F90 |
|---|
| 31 | |
|---|
| 32 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: sl |
|---|
| 33 | |
|---|
| 34 | !! * Substitutions |
|---|
| 35 | # include "domzgr_substitute.h90" |
|---|
| 36 | # include "vectopt_loop_substitute.h90" |
|---|
| 37 | !!---------------------------------------------------------------------- |
|---|
| 38 | !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) |
|---|
| 39 | !! $Id$ |
|---|
| 40 | !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) |
|---|
| 41 | !!---------------------------------------------------------------------- |
|---|
| 42 | |
|---|
| 43 | CONTAINS |
|---|
| 44 | |
|---|
| 45 | SUBROUTINE tra_adv_qck( kt, cdtype, ktra, pun, pvn, pwn, & |
|---|
| 46 | & ptb, ptn, pta ) |
|---|
| 47 | !!---------------------------------------------------------------------- |
|---|
| 48 | !! *** ROUTINE tra_adv_qck *** |
|---|
| 49 | !! |
|---|
| 50 | !! ** Purpose : Compute the now trend due to the advection of tracers |
|---|
| 51 | !! and add it to the general trend of passive tracer equations. |
|---|
| 52 | !! |
|---|
| 53 | !! ** Method : The advection is evaluated by a third order scheme |
|---|
| 54 | !! For a positive velocity u : |
|---|
| 55 | !! |
|---|
| 56 | !! i-1 i i+1 i+2 |
|---|
| 57 | !! |
|---|
| 58 | !! |--FU--|--FC--|--FD--|------| |
|---|
| 59 | !! u(i)>0 |
|---|
| 60 | !! |
|---|
| 61 | !! For a negative velocity u : |
|---|
| 62 | !! |
|---|
| 63 | !! |------|--FD--|--FC--|--FU--| |
|---|
| 64 | !! u(i)<0 |
|---|
| 65 | !! |
|---|
| 66 | !! FU is the second upwind point |
|---|
| 67 | !! FD is the first douwning point |
|---|
| 68 | !! FC is the central point (or the first upwind point) |
|---|
| 69 | !! |
|---|
| 70 | !! Flux(i) = u(i) * {0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i)Â?)/6)(FU+FD-2FC)} |
|---|
| 71 | !! with C(i)=|u(i)|dx(i)/dt (Courant number) |
|---|
| 72 | !! |
|---|
| 73 | !! dt = 2*rdtra and the scalar values are tb and sb |
|---|
| 74 | !! |
|---|
| 75 | !! On the vertical, the simple centered scheme used tn and sn |
|---|
| 76 | !! |
|---|
| 77 | !! The fluxes are bounded by the ULTIMATE limiter |
|---|
| 78 | !! to guarantee the monotonicity of the solution and to |
|---|
| 79 | !! prevent the appearance of spurious numerical oscillations |
|---|
| 80 | !! |
|---|
| 81 | !! ** Action : - update (ta,sa) with the now advective tracer trends |
|---|
| 82 | !! - save the trends in (ttrdh,strdh) ('key_trdtra') |
|---|
| 83 | !! |
|---|
| 84 | !! References : Leonard (1979, 1991) |
|---|
| 85 | !!---------------------------------------------------------------------- |
|---|
| 86 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
|---|
| 87 | CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) |
|---|
| 88 | INTEGER , INTENT(in ) :: ktra ! tracer index |
|---|
| 89 | REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn, pwn ! 3 ocean velocity components |
|---|
| 90 | REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptb, ptn ! before and now tracer fields |
|---|
| 91 | REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend |
|---|
| 92 | !! |
|---|
| 93 | INTEGER :: ji, jj, jk ! dummy loop indices |
|---|
| 94 | REAL(wp) :: zdt, z2, zcoef1 ! temporary scalars |
|---|
| 95 | REAL(wp) :: zbtr ! temporary scalars |
|---|
| 96 | !!---------------------------------------------------------------------- |
|---|
| 97 | |
|---|
| 98 | IF( kt == nit000 ) THEN |
|---|
| 99 | IF(lwp) WRITE(numout,*) |
|---|
| 100 | IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme' |
|---|
| 101 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' |
|---|
| 102 | ENDIF |
|---|
| 103 | |
|---|
| 104 | IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ! euler time-stepping |
|---|
| 105 | ELSE ; z2 = 2. ! leap-frog time-stepping |
|---|
| 106 | ENDIF |
|---|
| 107 | |
|---|
| 108 | |
|---|
| 109 | ! I. Slope estimation at the T-point for the limiter ULTIMATE |
|---|
| 110 | ! SL = Sum(1/C_out) with C_out : Courant number for the outflows |
|---|
| 111 | !--------------------------------------------------------------- |
|---|
| 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 ! |
|---|
| 115 | sl(:,:,:) = 100. |
|---|
| 116 | |
|---|
| 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 |
|---|
| 136 | |
|---|
| 137 | |
|---|
| 138 | ! II. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme |
|---|
| 139 | !--------------------------------------------------------------------------- |
|---|
| 140 | |
|---|
| 141 | CALL tra_adv_qck_hor( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) |
|---|
| 142 | |
|---|
| 143 | ! ! control print |
|---|
| 144 | IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - had: ', mask1=tmask, clinfo3=cdtype ) |
|---|
| 145 | |
|---|
| 146 | |
|---|
| 147 | ! III. The vertical fluxes are computed with the 2nd order centered scheme |
|---|
| 148 | !------------------------------------------------------------------------- |
|---|
| 149 | |
|---|
| 150 | CALL tra_adv_qck_ver( pwn, ptn , pta ) |
|---|
| 151 | |
|---|
| 152 | ! ! control print |
|---|
| 153 | IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - zad: ', mask1=tmask, clinfo3=cdtype ) |
|---|
| 154 | ! |
|---|
| 155 | END SUBROUTINE tra_adv_qck |
|---|
| 156 | |
|---|
| 157 | |
|---|
| 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 | !! |
|---|
| 164 | !!---------------------------------------------------------------------- |
|---|
| 165 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
|---|
| 166 | CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) |
|---|
| 167 | INTEGER , INTENT(in ) :: ktra ! tracer index |
|---|
| 168 | REAL(wp) , INTENT(in ) :: pz2 ! =1 or 2 (euler or leap-frog) |
|---|
| 169 | REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! horizontal effective velocity |
|---|
| 170 | REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptb ! before tracer field |
|---|
| 171 | REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend |
|---|
| 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 |
|---|
| 178 | !!---------------------------------------------------------------------- |
|---|
| 179 | |
|---|
| 180 | !---------------------------------------------------------------------- |
|---|
| 181 | ! 0. Initialization (should ot be needed on the whole array ???) |
|---|
| 182 | !---------------------------------------------------------------------- |
|---|
| 183 | zmask(:,:,:)= 0.e0 |
|---|
| 184 | zlap (:,:,:)= 0.e0 |
|---|
| 185 | zdwst(:,:,:)= 0.e0 |
|---|
| 186 | zlim (:,:,:)= 0.e0 |
|---|
| 187 | |
|---|
| 188 | !---------------------------------------------------------------------- |
|---|
| 189 | ! I. Part 1 : x-direction |
|---|
| 190 | !---------------------------------------------------------------------- |
|---|
| 191 | |
|---|
| 192 | DO jk = 1, jpkm1 |
|---|
| 193 | |
|---|
| 194 | ! Laplacian of tracers (at before time step) |
|---|
| 195 | ! ------------------------------------------ |
|---|
| 196 | !--- First derivative (gradient) |
|---|
| 197 | DO jj = 1, jpjm1 |
|---|
| 198 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 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) ) |
|---|
| 207 | END DO |
|---|
| 208 | END DO |
|---|
| 209 | !--- Function lim=FU+SL*(FC-FU) used by the limiter |
|---|
| 210 | !--- Computation of the ustream and downstream lim at the T-points |
|---|
| 211 | !!gm bug : fs_2 instead of 2 ... |
|---|
| 212 | !!gm a lot of optimisation to be done in this routine.... |
|---|
| 213 | DO jj = 2, jpjm1 |
|---|
| 214 | DO ji = 2, fs_jpim1 ! vector opt. |
|---|
| 215 | ! Upstream in the x-direction for the tracer |
|---|
| 216 | zmask(ji,jj,jk) = ptb(ji-1,jj,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji-1,jj,jk) ) |
|---|
| 217 | ! Downstream in the x-direction for the tracer |
|---|
| 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 |
|---|
| 230 | DO jj = 1, jpjm1 |
|---|
| 231 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 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) |
|---|
| 234 | ! Mask at the T-points in the x-direction (mask=0 or mask=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 | |
|---|
| 241 | |
|---|
| 242 | ! Horizontal advective fluxes |
|---|
| 243 | ! --------------------------- |
|---|
| 244 | DO jk = 1, jpkm1 |
|---|
| 245 | zdt = pz2 * rdttra(jk) |
|---|
| 246 | !--- tracer flux at u and v-points |
|---|
| 247 | DO jj = 1, jpjm1 |
|---|
| 248 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 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 |
|---|
| 256 | |
|---|
| 257 | !--- QUICKEST scheme |
|---|
| 258 | ! Temperature on the x-direction |
|---|
| 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 ) |
|---|
| 264 | !--- If the second ustream point is a land point |
|---|
| 265 | !--- the flux is computed by the 1st order UPWIND scheme |
|---|
| 266 | 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 |
|---|
| 269 | END DO |
|---|
| 270 | END DO |
|---|
| 271 | |
|---|
| 272 | !--- Tracer flux divergence at t-point added to the general trend |
|---|
| 273 | DO jj = 2, jpjm1 |
|---|
| 274 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 275 | !--- horizontal advective trends |
|---|
| 276 | zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) |
|---|
| 277 | !--- add it to the general tracer trends |
|---|
| 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 ) |
|---|
| 285 | |
|---|
| 286 | |
|---|
| 287 | !---------------------------------------------------------------------- |
|---|
| 288 | ! I. Part 2 : y-direction |
|---|
| 289 | !---------------------------------------------------------------------- |
|---|
| 290 | |
|---|
| 291 | DO jk = 1, jpkm1 |
|---|
| 292 | |
|---|
| 293 | ! Laplacian of tracers (at before time step) |
|---|
| 294 | ! ------------------------------------------ |
|---|
| 295 | !--- First derivative (gradient) |
|---|
| 296 | DO jj = 1, jpjm1 |
|---|
| 297 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 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) |
|---|
| 300 | END DO |
|---|
| 301 | END DO |
|---|
| 302 | !--- Second derivative (divergence) |
|---|
| 303 | DO jj = 2, jpjm1 |
|---|
| 304 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 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) ) |
|---|
| 307 | END DO |
|---|
| 308 | END DO |
|---|
| 309 | !--- Function lim=FU+SL*(FC-FU) used by the limiter |
|---|
| 310 | !--- Computation of the ustream and downstream lim at the T-points |
|---|
| 311 | DO jj = 2, jpjm1 |
|---|
| 312 | DO ji = 2, fs_jpim1 ! vector opt. |
|---|
| 313 | ! Upstream in the y-direction for the tracer |
|---|
| 314 | zmask(ji,jj,jk) = ptb(ji,jj-1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj-1,jk) ) |
|---|
| 315 | ! Downstream in the y-direction for the tracer |
|---|
| 316 | zdwst(ji,jj,jk) = ptb(ji,jj+1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj+1,jk) ) |
|---|
| 317 | ENDDO |
|---|
| 318 | ENDDO |
|---|
| 319 | END DO |
|---|
| 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 | |
|---|
| 326 | ! --- lim at the V-points in function of the direction of the flow |
|---|
| 327 | ! ---------------------------------------------------------------- |
|---|
| 328 | DO jj = 1, jpjm1 |
|---|
| 329 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 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) |
|---|
| 332 | ! Mask at the T-points in the y-direction (mask=0 or mask=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 |
|---|
| 338 | |
|---|
| 339 | ! Horizontal advective fluxes |
|---|
| 340 | ! --------------------------- |
|---|
| 341 | |
|---|
| 342 | DO jk = 1, jpkm1 |
|---|
| 343 | zdt = pz2 * rdttra(jk) |
|---|
| 344 | !--- tracer flux at u and v-points |
|---|
| 345 | DO jj = 1, jpjm1 |
|---|
| 346 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 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 |
|---|
| 354 | |
|---|
| 355 | !--- QUICKEST scheme |
|---|
| 356 | ! Temperature on the y-direction |
|---|
| 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 ) |
|---|
| 362 | !--- If the second ustream point is a land point |
|---|
| 363 | !--- the flux is computed by the 1st order UPWIND scheme |
|---|
| 364 | 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 |
|---|
| 367 | END DO |
|---|
| 368 | END DO |
|---|
| 369 | |
|---|
| 370 | !--- Tracer flux divergence at t-point added to the general trend |
|---|
| 371 | DO jj = 2, jpjm1 |
|---|
| 372 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 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 |
|---|
| 383 | IF( ln_diaptr .AND. cdtype == 'TRA' .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN |
|---|
| 384 | IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zdwst(:,:,:) ) |
|---|
| 385 | IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zdwst(:,:,:) ) |
|---|
| 386 | ENDIF |
|---|
| 387 | ! |
|---|
| 388 | END SUBROUTINE tra_adv_qck_hor |
|---|
| 389 | |
|---|
| 390 | |
|---|
| 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 | !!---------------------------------------------------------------------- |
|---|
| 405 | |
|---|
| 406 | ! Vertical advection |
|---|
| 407 | ! ------------------ |
|---|
| 408 | |
|---|
| 409 | ! 1. Vertical advective fluxes |
|---|
| 410 | ! ---------------------------- |
|---|
| 411 | |
|---|
| 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 |
|---|
| 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...) |
|---|
| 421 | |
|---|
| 422 | ! Second order centered tracer flux at w-point |
|---|
| 423 | DO jk = 2, jpkm1 |
|---|
| 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 |
|---|
| 429 | !--- Second order centered scheme |
|---|
| 430 | sl(ji,jj,jk) = pwn(ji,jj,jk) * ( zfc + zfd ) / ( fse3t(ji,jj,jk-1) + fse3t(ji,jj,jk) ) |
|---|
| 431 | END DO |
|---|
| 432 | END DO |
|---|
| 433 | END DO |
|---|
| 434 | |
|---|
| 435 | ! 2. Tracer flux divergence at t-point added to the general trend |
|---|
| 436 | ! --------------------------------------------------------------- |
|---|
| 437 | DO jk = 1, jpkm1 |
|---|
| 438 | DO jj = 2, jpjm1 |
|---|
| 439 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 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) ) |
|---|
| 443 | END DO |
|---|
| 444 | END DO |
|---|
| 445 | END DO |
|---|
| 446 | ! |
|---|
| 447 | END SUBROUTINE tra_adv_qck_ver |
|---|
| 448 | |
|---|
| 449 | |
|---|
| 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 | ! |
|---|
| 465 | END FUNCTION |
|---|
| 466 | |
|---|
| 467 | !!====================================================================== |
|---|
| 468 | END MODULE traadv_qck |
|---|
