Changeset 5948 for branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
- Timestamp:
- 2015-11-30T11:47:24+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r5947 r5948 4 4 !! Ocean tracers: horizontal component of the lateral tracer mixing trend 5 5 !!====================================================================== 6 !! History : OPA ! 1994-08 (G. Madec, M. Imbard) 7 !! 8.0 ! 1997-05 (G. Madec) split into traldf and trazdf 8 !! NEMO ! 2002-08 (G. Madec) Free form, F90 9 !! 1.0 ! 2005-11 (G. Madec) merge traldf and trazdf :-) 10 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 6 !! History : OPA ! 1994-08 (G. Madec, M. Imbard) 7 !! 8.0 ! 1997-05 (G. Madec) split into traldf and trazdf 8 !! NEMO ! 2002-08 (G. Madec) Free form, F90 9 !! 1.0 ! 2005-11 (G. Madec) merge traldf and trazdf :-) 10 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 11 !! 3.7 ! 2014-01 (G. Madec, S. Masson) restructuration/simplification of aht/aeiv specification 12 !! - ! 2014-02 (F. Lemarie, G. Madec) triad operator (Griffies) + Method of Stabilizing Correction 11 13 !!---------------------------------------------------------------------- 12 #if defined key_ldfslp || defined key_esopa 14 13 15 !!---------------------------------------------------------------------- 14 !! 'key_ldfslp' slope of the lateral diffusive direction 15 !!---------------------------------------------------------------------- 16 !! tra_ldf_iso : update the tracer trend with the horizontal 17 !! component of a iso-neutral laplacian operator 18 !! and with the vertical part of 19 !! the isopycnal or geopotential s-coord. operator 16 !! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator 20 18 !!---------------------------------------------------------------------- 21 19 USE oce ! ocean dynamics and active tracers … … 23 21 USE trc_oce ! share passive tracers/Ocean variables 24 22 USE zdf_oce ! ocean vertical physics 25 USE ldftra _oce ! ocean active tracers: lateral physics23 USE ldftra ! lateral diffusion: tracer eddy coefficients 26 24 USE ldfslp ! iso-neutral slopes 27 25 USE diaptr ! poleward transport diagnostics 26 ! 28 27 USE in_out_manager ! I/O manager 29 28 USE iom ! I/O library … … 40 39 !! * Substitutions 41 40 # include "domzgr_substitute.h90" 42 # include "ldftra_substitute.h90"43 41 # include "vectopt_loop_substitute.h90" 44 42 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)43 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 46 44 !! $Id$ 47 45 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 49 47 CONTAINS 50 48 51 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv,&52 & pgui, pgvi,&53 & ptb, pta, kjpt, pahtb0)49 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 50 & pgui, pgvi, & 51 & ptb , ptbb, pta , kjpt, kpass ) 54 52 !!---------------------------------------------------------------------- 55 53 !! *** ROUTINE tra_ldf_iso *** … … 66 64 !! 67 65 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 68 !! ======== with partial cell update if ln_zps=T. 66 !! ======== with partial cell update if ln_zps=T 67 !! with top cell update if ln_isfcav 69 68 !! 70 69 !! 2nd part : horizontal fluxes of the lateral mixing operator 71 70 !! ======== 72 !! zftu = (aht+ahtb0)e2u*e3u/e1u di[ tb ]73 !! - ahte2u*uslp dk[ mi(mk(tb)) ]74 !! zftv = (aht+ahtb0)e1v*e3v/e2v dj[ tb ]75 !! - ahte2u*vslp dk[ mj(mk(tb)) ]71 !! zftu = pahu e2u*e3u/e1u di[ tb ] 72 !! - pahu e2u*uslp dk[ mi(mk(tb)) ] 73 !! zftv = pahv e1v*e3v/e2v dj[ tb ] 74 !! - pahv e2u*vslp dk[ mj(mk(tb)) ] 76 75 !! take the horizontal divergence of the fluxes: 77 !! difft = 1/(e1 t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] }76 !! difft = 1/(e1e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 78 77 !! Add this trend to the general trend (ta,sa): 79 78 !! ta = ta + difft … … 82 81 !! ======== (excluding the vertical flux proportional to dk[t] ) 83 82 !! vertical fluxes associated with the rotated lateral mixing: 84 !! zftw = -aht {e2t*wslpi di[ mi(mk(tb)) ]85 !! +e1t*wslpj dj[ mj(mk(tb)) ] }83 !! zftw = - { mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ] 84 !! + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ] } 86 85 !! take the horizontal divergence of the fluxes: 87 !! difft = 1/(e1 t*e2t*e3t) dk[ zftw ]86 !! difft = 1/(e1e2t*e3t) dk[ zftw ] 88 87 !! Add this trend to the general trend (ta,sa): 89 88 !! pta = pta + difft … … 91 90 !! ** Action : Update pta arrays with the before rotated diffusion 92 91 !!---------------------------------------------------------------------- 93 USE oce , ONLY: zftu => ua , zftv => va ! (ua,va) used as workspace94 !95 92 INTEGER , INTENT(in ) :: kt ! ocean time-step index 96 INTEGER , INTENT(in ) :: kit000 93 INTEGER , INTENT(in ) :: kit000 ! first time step index 97 94 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 98 95 INTEGER , INTENT(in ) :: kjpt ! number of tracers 99 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at pstep levels 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 103 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 96 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 97 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 98 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 99 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 100 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! tracer (kpass=1) or laplacian of tracer (kpass=2) 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptbb ! tracer (only used in kpass=2) 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 104 103 ! 105 104 INTEGER :: ji, jj, jk, jn ! dummy loop indices 106 INTEGER :: ikt 107 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 108 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 109 REAL(wp) :: zcoef0, zbtr, ztra ! - - 110 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdkt, zdk1t, zdit, zdjt, ztfw 105 INTEGER :: ierr ! local integer 106 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 107 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 108 REAL(wp) :: zcoef0, ze3w_2, zsign, z2dt, z1_2dt ! - - 109 #if defined key_diaar5 110 REAL(wp) :: zztmp ! local scalar 111 #endif 112 REAL(wp), POINTER, DIMENSION(:,:) :: zdkt, zdk1t, z2d 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw 112 114 !!---------------------------------------------------------------------- 113 115 ! 114 116 IF( nn_timing == 1 ) CALL timing_start('tra_ldf_iso') 115 117 ! 116 CALL wrk_alloc( jpi, jpj, z2d ) 117 CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 118 ! 119 118 CALL wrk_alloc( jpi,jpj, zdkt, zdk1t, z2d ) 119 CALL wrk_alloc( jpi,jpj,jpk, zdit, zdjt , zftu, zftv, ztfw ) 120 ! 120 121 IF( kt == kit000 ) THEN 121 122 IF(lwp) WRITE(numout,*) 122 123 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 123 124 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 125 ! 126 akz (:,:,:) = 0._wp 127 ah_wslp2(:,:,:) = 0._wp 128 ENDIF 129 ! 130 ! ! set time step size (Euler/Leapfrog) 131 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdttra(1) ! at nit000 (Euler) 132 ELSE ; z2dt = 2.* rdttra(1) ! (Leapfrog) 133 ENDIF 134 z1_2dt = 1._wp / z2dt 135 ! 136 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) 137 ELSE ; zsign = -1._wp 138 ENDIF 139 140 141 !!---------------------------------------------------------------------- 142 !! 0 - calculate ah_wslp2 and akz 143 !!---------------------------------------------------------------------- 144 ! 145 IF( kpass == 1 ) THEN !== first pass only ==! 146 ! 147 DO jk = 2, jpkm1 148 DO jj = 2, jpjm1 149 DO ji = fs_2, fs_jpim1 ! vector opt. 150 ! 151 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 152 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 153 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 154 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 155 ! 156 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 157 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 158 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 159 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 160 ! 161 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 162 & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 163 END DO 164 END DO 165 END DO 166 ! 167 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient 168 DO jk = 2, jpkm1 169 DO jj = 2, jpjm1 170 DO ji = fs_2, fs_jpim1 171 akz(ji,jj,jk) = 0.25_wp * ( & 172 & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 173 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 174 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 175 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 176 END DO 177 END DO 178 END DO 179 ! 180 IF( ln_traldf_blp ) THEN ! bilaplacian operator 181 DO jk = 2, jpkm1 182 DO jj = 1, jpjm1 183 DO ji = 1, fs_jpim1 184 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 185 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) ) 186 END DO 187 END DO 188 END DO 189 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 190 DO jk = 2, jpkm1 191 DO jj = 1, jpjm1 192 DO ji = 1, fs_jpim1 193 ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 194 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 195 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 196 END DO 197 END DO 198 END DO 199 ENDIF 200 ! 201 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 202 akz(:,:,:) = ah_wslp2(:,:,:) 203 ENDIF 124 204 ENDIF 125 205 ! … … 131 211 !! I - masked horizontal derivative 132 212 !!---------------------------------------------------------------------- 133 !!bug ajout.... why? (1,jpj,:) and (jpi,1,:) should be sufficient....134 zdit (1,:,:) = 0. e0 ; zdit (jpi,:,:) = 0.e0135 zdjt (1,:,:) = 0. e0 ; zdjt (jpi,:,:) = 0.e0213 !!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient.... 214 zdit (1,:,:) = 0._wp ; zdit (jpi,:,:) = 0._wp 215 zdjt (1,:,:) = 0._wp ; zdjt (jpi,:,:) = 0._wp 136 216 !!end 137 217 … … 145 225 END DO 146 226 END DO 147 148 ! partial cell correction 149 IF( ln_zps ) THEN ! partial steps correction at the last ocean level 150 DO jj = 1, jpjm1 227 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 228 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 151 229 DO ji = 1, fs_jpim1 ! vector opt. 152 ! IF useless if zpshde defines pgu everywhere230 !!gm the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere 153 231 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 154 232 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 155 233 END DO 156 234 END DO 235 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 236 DO jj = 1, jpjm1 237 DO ji = 1, fs_jpim1 ! vector opt. 238 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 239 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 240 END DO 241 END DO 242 ENDIF 157 243 ENDIF 158 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the first wet level beneath a cavity 159 DO jj = 1, jpjm1 244 245 !!---------------------------------------------------------------------- 246 !! II - horizontal trend (full) 247 !!---------------------------------------------------------------------- 248 ! 249 DO jk = 1, jpkm1 ! Horizontal slab 250 ! 251 ! !== Vertical tracer gradient 252 zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1 253 ! 254 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 255 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 256 ENDIF 257 !!gm I don't understand why we should need this.... since wmask is used instead of tmask 258 ! IF ( ln_isfcav ) THEN 259 ! DO jj = 1, jpj 260 ! DO ji = 1, jpi ! vector opt. 261 ! ikt = mikt(ji,jj) ! surface level 262 ! zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 263 ! zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 264 ! END DO 265 ! END DO 266 ! END IF 267 !!gm 268 269 DO jj = 1 , jpjm1 !== Horizontal fluxes 160 270 DO ji = 1, fs_jpim1 ! vector opt. 161 IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 162 IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 163 END DO 164 END DO 165 END IF 166 167 !!---------------------------------------------------------------------- 168 !! II - horizontal trend (full) 169 !!---------------------------------------------------------------------- 170 !!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t ) 171 ! 1. Vertical tracer gradient at level jk and jk+1 172 ! ------------------------------------------------ 173 ! 174 ! interior value 175 DO jk = 2, jpkm1 176 DO jj = 1, jpj 177 DO ji = 1, jpi ! vector opt. 178 zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 179 ! 180 zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn ) ) * wmask(ji,jj,jk) 181 END DO 182 END DO 183 END DO 184 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 185 zdk1t(:,:,1) = ( ptb(:,:,1,jn ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 186 zdkt (:,:,1) = zdk1t(:,:,1) 187 IF ( ln_isfcav ) THEN 188 DO jj = 1, jpj 189 DO ji = 1, jpi ! vector opt. 190 ikt = mikt(ji,jj) ! surface level 191 zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 192 zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 193 END DO 194 END DO 195 END IF 196 197 ! 2. Horizontal fluxes 198 ! -------------------- 199 DO jk = 1, jpkm1 200 DO jj = 1 , jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 203 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 271 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk) 272 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) 204 273 ! 205 274 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) & … … 209 278 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. ) 210 279 ! 211 zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku212 zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv280 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 281 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 213 282 ! 214 283 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 215 & + zcof1 * ( zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk) &216 & + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk) ) ) * umask(ji,jj,jk)284 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 285 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 217 286 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 218 & + zcof2 * ( zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk) & 219 & + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk) ) ) * vmask(ji,jj,jk) 220 END DO 221 END DO 222 223 ! II.4 Second derivative (divergence) and add to the general trend 224 ! ---------------------------------------------------------------- 225 DO jj = 2 , jpjm1 287 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 288 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 289 END DO 290 END DO 291 ! 292 DO jj = 2 , jpjm1 !== horizontal divergence and add to pta 226 293 DO ji = fs_2, fs_jpim1 ! vector opt. 227 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 228 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 229 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 230 END DO 231 END DO 232 ! ! =============== 294 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 295 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 296 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 297 END DO 298 END DO 233 299 END DO ! End of slab 234 ! ! =============== 235 ! 236 ! "Poleward" diffusive heat or salt transports (T-S case only) 237 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 238 ! note sign is reversed to give down-gradient diffusive transports (#1043) 239 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 240 IF( jn == jp_sal) str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 241 ENDIF 242 243 IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 244 ! 245 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 246 z2d(:,:) = 0._wp 247 DO jk = 1, jpkm1 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 ! vector opt. 250 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 251 END DO 252 END DO 253 END DO 254 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 255 CALL lbc_lnk( z2d, 'U', -1. ) 256 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 257 ! 258 z2d(:,:) = 0._wp 259 DO jk = 1, jpkm1 260 DO jj = 2, jpjm1 261 DO ji = fs_2, fs_jpim1 ! vector opt. 262 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 263 END DO 264 END DO 265 END DO 266 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 267 CALL lbc_lnk( z2d, 'V', -1. ) 268 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 269 END IF 270 ! 271 ENDIF 272 273 !!---------------------------------------------------------------------- 274 !! III - vertical trend of T & S (extra diagonal terms only) 275 !!---------------------------------------------------------------------- 276 277 ! Local constant initialization 278 ! ----------------------------- 279 ztfw(1,:,:) = 0.e0 ; ztfw(jpi,:,:) = 0.e0 300 301 302 !!---------------------------------------------------------------------- 303 !! III - vertical trend (full) 304 !!---------------------------------------------------------------------- 305 306 ztfw(1,:,:) = 0._wp ; ztfw(jpi,:,:) = 0._wp 280 307 281 308 ! Vertical fluxes … … 283 310 284 311 ! Surface and bottom vertical fluxes set to zero 285 ztfw(:,:, 1 ) = 0. e0 ; ztfw(:,:,jpk) = 0.e0312 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 286 313 287 314 ! interior (2=<jk=<jpk-1) … … 289 316 DO jj = 2, jpjm1 290 317 DO ji = fs_2, fs_jpim1 ! vector opt. 291 zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 292 ! 293 zmsku = 1./MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 294 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk), 1. ) 295 zmskv = 1./MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 296 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk), 1. ) 297 ! 298 zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) 299 zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 318 ! 319 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 320 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 321 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 322 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 323 ! 324 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 325 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 326 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 327 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 328 ! 329 zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked 330 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 300 331 ! 301 332 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & … … 306 337 END DO 307 338 END DO 308 309 310 ! I.5 Divergence of vertical fluxes added to the general tracer trend 311 ! ------------------------------------------------------------------- 312 DO jk = 1, jpkm1 339 ! 340 ! !== add the vertical 33 flux ==! 341 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 342 DO jk = 2, jpkm1 343 DO jj = 1, jpjm1 344 DO ji = fs_2, fs_jpim1 345 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) & 346 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 347 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 348 END DO 349 END DO 350 END DO 351 ! 352 ELSE ! bilaplacian 353 SELECT CASE( kpass ) 354 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 355 DO jk = 2, jpkm1 356 DO jj = 1, jpjm1 357 DO ji = fs_2, fs_jpim1 358 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 359 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 360 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk) 361 END DO 362 END DO 363 END DO 364 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. 365 DO jk = 2, jpkm1 366 DO jj = 1, jpjm1 367 DO ji = fs_2, fs_jpim1 368 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) & 369 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 370 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) 371 END DO 372 END DO 373 END DO 374 END SELECT 375 ENDIF 376 ! 377 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! 313 378 DO jj = 2, jpjm1 314 379 DO ji = fs_2, fs_jpim1 ! vector opt. 315 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 316 ztra = ( ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1) ) * zbtr 317 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 380 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 381 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 318 382 END DO 319 383 END DO 320 384 END DO 321 385 ! 322 END DO 323 ! 324 CALL wrk_dealloc( jpi, jpj, z2d ) 325 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 386 IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==! 387 ( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass (bilaplacian) ==! 388 ! 389 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 390 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 391 ! note sign is reversed to give down-gradient diffusive transports (#1043) 392 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 393 IF( jn == jp_sal) str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 394 ENDIF 395 ! 396 IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 397 ! 398 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 399 z2d(:,:) = zftu(ji,jj,1) 400 DO jk = 2, jpkm1 401 DO jj = 2, jpjm1 402 DO ji = fs_2, fs_jpim1 ! vector opt. 403 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 404 END DO 405 END DO 406 END DO 407 !!gm CAUTION I think there is an error of sign when using BLP operator.... 408 !!gm a multiplication by zsign is required (to be checked twice !) 409 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 410 CALL lbc_lnk( z2d, 'U', -1. ) 411 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 412 ! 413 z2d(:,:) = zftv(ji,jj,1) 414 DO jk = 2, jpkm1 415 DO jj = 2, jpjm1 416 DO ji = fs_2, fs_jpim1 ! vector opt. 417 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 418 END DO 419 END DO 420 END DO 421 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 422 CALL lbc_lnk( z2d, 'V', -1. ) 423 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 424 END IF 425 ! 426 ENDIF 427 ! 428 ENDIF !== end pass selection ==! 429 ! 430 ! ! =============== 431 END DO ! end tracer loop 432 ! ! =============== 433 ! 434 CALL wrk_dealloc( jpi, jpj, zdkt, zdk1t, z2d ) 435 CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw ) 326 436 ! 327 437 IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_iso') 328 438 ! 329 439 END SUBROUTINE tra_ldf_iso 330 331 #else332 !!----------------------------------------------------------------------333 !! default option : Dummy code NO rotation of the diffusive tensor334 !!----------------------------------------------------------------------335 CONTAINS336 SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, pgui, pgvi, ptb, pta, kjpt, pahtb0 ) ! Empty routine337 INTEGER:: kt, kit000338 CHARACTER(len=3) :: cdtype339 REAL, DIMENSION(:,:,:) :: pgu, pgv, pgui, pgvi ! tracer gradient at pstep levels340 REAL, DIMENSION(:,:,:,:) :: ptb, pta341 WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype, &342 & pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0343 END SUBROUTINE tra_ldf_iso344 #endif345 440 346 441 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
