- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traldf_iso.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traldf_iso.F90
r12178 r12928 40 40 41 41 !! * Substitutions 42 # include " vectopt_loop_substitute.h90"42 # include "do_loop_substitute.h90" 43 43 !!---------------------------------------------------------------------- 44 44 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 48 48 CONTAINS 49 49 50 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,&51 & 52 & pt b , ptbb, pta , kjpt, kpass )50 SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv, & 51 & pgu , pgv , pgui, pgvi, & 52 & pt , pt2 , pt_rhs , kjpt , kpass ) 53 53 !!---------------------------------------------------------------------- 54 54 !! *** ROUTINE tra_ldf_iso *** … … 87 87 !! difft = 1/(e1e2t*e3t) dk[ zftw ] 88 88 !! Add this trend to the general trend (ta,sa): 89 !! pt a = pta+ difft90 !! 91 !! ** Action : Update pt aarrays with the before rotated diffusion89 !! pt_rhs = pt_rhs + difft 90 !! 91 !! ** Action : Update pt_rhs arrays with the before rotated diffusion 92 92 !!---------------------------------------------------------------------- 93 93 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 96 96 INTEGER , INTENT(in ) :: kjpt ! number of tracers 97 97 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 98 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 98 99 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 99 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 100 101 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt b! tracer (kpass=1) or laplacian of tracer (kpass=2)102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt bb! tracer (only used in kpass=2)103 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt a! tracer trend102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 103 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 104 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 104 105 ! 105 106 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 108 109 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 109 110 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 110 REAL(wp) :: zcoef0, ze3w_2, zsign , z2dt, z1_2dt! - -111 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 111 112 REAL(wp), DIMENSION(jpi,jpj) :: zdkt, zdk1t, z2d 112 113 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw … … 124 125 l_hst = .FALSE. 125 126 l_ptr = .FALSE. 126 IF( cdtype == 'TRA' .AND. ln_diaptr )l_ptr = .TRUE.127 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 127 128 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 128 129 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 129 130 ! 130 ! ! set time step size (Euler/Leapfrog)131 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler)132 ELSE ; z2dt = 2.* rdt ! (Leapfrog)133 ENDIF134 z1_2dt = 1._wp / z2dt135 131 ! 136 132 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 144 140 IF( kpass == 1 ) THEN !== first pass only ==! 145 141 ! 146 DO jk = 2, jpkm1 147 DO jj = 2, jpjm1 148 DO ji = fs_2, fs_jpim1 ! vector opt. 149 ! 150 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 151 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 152 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 153 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 154 ! 155 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 156 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 157 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 158 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 159 ! 160 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 161 & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 162 END DO 163 END DO 164 END DO 142 DO_3D_00_00( 2, jpkm1 ) 143 ! 144 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 145 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 146 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 147 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 148 ! 149 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 150 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 151 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 152 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 153 ! 154 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 155 & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 156 END_3D 165 157 ! 166 158 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient 167 DO jk = 2, jpkm1 168 DO jj = 2, jpjm1 169 DO ji = fs_2, fs_jpim1 170 akz(ji,jj,jk) = 0.25_wp * ( & 171 & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 172 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 173 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 174 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 175 END DO 176 END DO 177 END DO 159 DO_3D_00_00( 2, jpkm1 ) 160 akz(ji,jj,jk) = 0.25_wp * ( & 161 & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 162 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 163 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 164 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 165 END_3D 178 166 ! 179 167 IF( ln_traldf_blp ) THEN ! bilaplacian operator 180 DO jk = 2, jpkm1 181 DO jj = 1, jpjm1 182 DO ji = 1, fs_jpim1 183 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 184 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 185 END DO 186 END DO 187 END DO 168 DO_3D_10_10( 2, jpkm1 ) 169 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 170 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) ) ) 171 END_3D 188 172 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 189 DO jk = 2, jpkm1 190 DO jj = 1, jpjm1 191 DO ji = 1, fs_jpim1 192 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 193 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 194 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 195 END DO 196 END DO 197 END DO 173 DO_3D_10_10( 2, jpkm1 ) 174 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 175 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 176 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 177 END_3D 198 178 ENDIF 199 179 ! … … 216 196 217 197 ! Horizontal tracer gradient 218 DO jk = 1, jpkm1 219 DO jj = 1, jpjm1 220 DO ji = 1, fs_jpim1 ! vector opt. 221 zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 222 zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 223 END DO 224 END DO 225 END DO 198 DO_3D_10_10( 1, jpkm1 ) 199 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 200 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 201 END_3D 226 202 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 227 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 228 DO ji = 1, fs_jpim1 ! vector opt. 229 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 230 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 231 END DO 232 END DO 203 DO_2D_10_10 204 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 205 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 206 END_2D 233 207 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 234 DO jj = 1, jpjm1 235 DO ji = 1, fs_jpim1 ! vector opt. 236 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 237 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 238 END DO 239 END DO 208 DO_2D_10_10 209 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 210 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 211 END_2D 240 212 ENDIF 241 213 ENDIF … … 248 220 ! 249 221 ! !== Vertical tracer gradient 250 zdk1t(:,:) = ( pt b(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1222 zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1 251 223 ! 252 224 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 253 ELSE ; zdkt(:,:) = ( pt b(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)225 ELSE ; zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 254 226 ENDIF 255 DO jj = 1 , jpjm1 !== Horizontal fluxes 256 DO ji = 1, fs_jpim1 ! vector opt. 257 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 258 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 259 ! 260 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 261 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 262 ! 263 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 264 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 265 ! 266 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 267 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 268 ! 269 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 270 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 271 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 272 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 273 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 274 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 275 END DO 276 END DO 277 ! 278 DO jj = 2 , jpjm1 !== horizontal divergence and add to pta 279 DO ji = fs_2, fs_jpim1 ! vector opt. 280 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 281 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 282 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 283 END DO 284 END DO 227 DO_2D_10_10 228 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 229 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 230 ! 231 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 232 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 233 ! 234 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 235 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 236 ! 237 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 238 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 239 ! 240 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 241 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 242 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 243 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 244 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 245 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 246 END_2D 247 ! 248 DO_2D_00_00 249 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 250 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 251 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 252 END_2D 285 253 END DO ! End of slab 286 254 … … 288 256 !! III - vertical trend (full) 289 257 !!---------------------------------------------------------------------- 290 !291 ztfw(fs_2:1,:,:) = 0._wp ; ztfw(jpi:fs_jpim1,:,:) = 0._wp ! avoid to potentially manipulate NaN values292 258 ! 293 259 ! Vertical fluxes … … 296 262 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 297 263 298 DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) 299 DO jj = 2, jpjm1 300 DO ji = fs_2, fs_jpim1 ! vector opt. 301 ! 302 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 303 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 304 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 305 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 306 ! 307 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 308 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 309 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 310 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 311 ! 312 zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked 313 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 314 ! 315 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & 316 & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & 317 & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & 318 & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) 319 END DO 320 END DO 321 END DO 264 DO_3D_00_00( 2, jpkm1 ) 265 ! 266 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 267 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 268 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 269 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 270 ! 271 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 272 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 273 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 274 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 275 ! 276 zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked 277 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 278 ! 279 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & 280 & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & 281 & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & 282 & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) 283 END_3D 322 284 ! !== add the vertical 33 flux ==! 323 285 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 324 DO jk = 2, jpkm1 325 DO jj = 2, jpjm1 326 DO ji = fs_2, fs_jpim1 327 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 328 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 329 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 330 END DO 331 END DO 332 END DO 286 DO_3D_00_00( 2, jpkm1 ) 287 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 288 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 289 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 290 END_3D 333 291 ! 334 292 ELSE ! bilaplacian 335 293 SELECT CASE( kpass ) 336 294 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 337 DO jk = 2, jpkm1 338 DO jj = 2, jpjm1 339 DO ji = fs_2, fs_jpim1 340 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 341 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 342 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 343 END DO 344 END DO 345 END DO 346 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. 347 DO jk = 2, jpkm1 348 DO jj = 2, jpjm1 349 DO ji = fs_2, fs_jpim1 350 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 351 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 352 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) 353 END DO 354 END DO 355 END DO 295 DO_3D_00_00( 2, jpkm1 ) 296 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 297 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 298 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 299 END_3D 300 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp. 301 DO_3D_00_00( 2, jpkm1 ) 302 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 303 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) & 304 & + akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) ) 305 END_3D 356 306 END SELECT 357 307 ENDIF 358 308 ! 359 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! 360 DO jj = 2, jpjm1 361 DO ji = fs_2, fs_jpim1 ! vector opt. 362 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 363 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 364 END DO 365 END DO 366 END DO 309 DO_3D_00_00( 1, jpkm1 ) 310 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 311 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 312 END_3D 367 313 ! 368 314 IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==!
Note: See TracChangeset
for help on using the changeset viewer.
