Changeset 786 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
- Timestamp:
- 2008-01-10T18:11:23+01:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r719 r786 6 6 !! History : ! 06-00 (A.Estublier) for passive tracers 7 7 !! ! 01-08 (E.Durand, G.Madec) adapted for T & S 8 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 8 !! NEMO 1.0 ! 02-06 (G. Madec) F90: Free form and module 9 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC 9 10 !!---------------------------------------------------------------------- 10 11 … … 13 14 !! and vertical advection trends using MUSCL scheme 14 15 !!---------------------------------------------------------------------- 15 USE oce ! ocean dynamics and active tracers16 16 USE dom_oce ! ocean space and time domain 17 17 USE trdmod ! ocean active tracers trends … … 34 34 # include "vectopt_loop_substitute.h90" 35 35 !!---------------------------------------------------------------------- 36 !! OPA 9.0 , LOCEAN-IPSL (2006)37 !! $ Header$36 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 37 !! $Id:$ 38 38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 39 39 !!---------------------------------------------------------------------- … … 41 41 CONTAINS 42 42 43 SUBROUTINE tra_adv_muscl( kt, pun, pvn, pwn ) 43 SUBROUTINE tra_adv_muscl( kt, cdtype, ktra, pun, pvn, pwn, & 44 & ptb , pta ) 44 45 !!---------------------------------------------------------------------- 45 46 !! *** ROUTINE tra_adv_muscl *** … … 51 52 !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries 52 53 !! 53 !! ** Action : - update ( ta,sa) with the now advective tracer trends54 !! ** Action : - update (pta,sa) with the now advective tracer trends 54 55 !! - save trends in (ztrdt,ztrds) ('key_trdtra') 55 56 !! … … 57 58 !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 58 59 !!---------------------------------------------------------------------- 59 USE oce, ONLY : ztrdt => ua ! use ua as workspace 60 USE oce, ONLY : ztrds => va ! use va as workspace 61 !! 62 INTEGER , INTENT(in) :: kt ! ocean time-step index 63 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component 64 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component 65 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component 66 !! 67 INTEGER :: ji, jj, jk ! dummy loop indices 68 REAL(wp) :: & 69 zu, zv, zw, zeu, zev, & 70 zew, zbtr, zstep, & 71 z0u, z0v, z0w, & 72 zzt1, zzt2, zalpha, & 73 zzs1, zzs2, z2, & 74 zta, zsa, & 75 z_hdivn_x, z_hdivn_y, z_hdivn 76 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zt1, zt2, ztp1, ztp2 ! 3D workspace 77 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zs1, zs2, zsp1, zsp2 ! " " 60 INTEGER , INTENT(in ) :: kt ! ocean time-step index 61 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 62 INTEGER , INTENT(in ) :: ktra ! tracer index 63 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn, pwn ! 3 ocean velocity components 64 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptb ! before tracer fields 65 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 66 !! 67 INTEGER :: ji, jj, jk ! dummy loop indices 68 REAL(wp) :: zu, zv, zw, zeu, zev 69 REAL(wp) :: zew, zbtr, z2, zstep 70 REAL(wp) :: z0u, z0v, z0w 71 REAL(wp) :: zzwx, zzwy, zalpha 72 REAL(wp) :: zta 73 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwx, zwy, zslpx, zslpy ! 3D workspace 78 74 !!---------------------------------------------------------------------- 79 75 … … 95 91 DO jj = 1, jpjm1 96 92 DO ji = 1, fs_jpim1 ! vector opt. 97 zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) 98 zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) 99 zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) 100 zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) 93 zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk) - ptb(ji,jj,jk) ) 94 zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk) - ptb(ji,jj,jk) ) 101 95 END DO 102 96 END DO 103 97 END DO 104 98 ! bottom values 105 zt1(:,:,jpk) = 0.e0 ; zt2(:,:,jpk) = 0.e0 106 zs1(:,:,jpk) = 0.e0 ; zs2(:,:,jpk) = 0.e0 107 108 ! lateral boundary conditions on zt1, zt2 ; zs1, zs2 (changed sign) 109 CALL lbc_lnk( zt1, 'U', -1. ) ; CALL lbc_lnk( zs1, 'U', -1. ) 110 CALL lbc_lnk( zt2, 'V', -1. ) ; CALL lbc_lnk( zs2, 'V', -1. ) 99 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 100 101 ! lateral boundary conditions on zwx, zwy (changed sign) 102 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) 111 103 112 104 ! Slopes … … 115 107 DO jj = 2, jpj 116 108 DO ji = fs_2, jpi ! vector opt. 117 ztp1(ji,jj,jk) = ( zt1(ji,jj,jk) + zt1(ji-1,jj ,jk) ) & 118 & * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj ,jk) ) ) 119 zsp1(ji,jj,jk) = ( zs1(ji,jj,jk) + zs1(ji-1,jj ,jk) ) & 120 & * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj ,jk) ) ) 121 ztp2(ji,jj,jk) = ( zt2(ji,jj,jk) + zt2(ji ,jj-1,jk) ) & 122 & * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji ,jj-1,jk) ) ) 123 zsp2(ji,jj,jk) = ( zs2(ji,jj,jk) + zs2(ji ,jj-1,jk) ) & 124 & * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji ,jj-1,jk) ) ) 109 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) & 110 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) ) 111 zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) & 112 & * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) ) 125 113 END DO 126 114 END DO 127 115 END DO 128 116 ! bottom values 129 ztp1(:,:,jpk) = 0.e0 ; ztp2(:,:,jpk) = 0.e0 130 zsp1(:,:,jpk) = 0.e0 ; zsp2(:,:,jpk) = 0.e0 117 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 131 118 132 119 ! Slopes limitation … … 134 121 DO jj = 2, jpj 135 122 DO ji = fs_2, jpi ! vector opt. 136 ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) ) & 137 & * MIN( ABS( ztp1(ji ,jj,jk) ), & 138 & 2.*ABS( zt1 (ji-1,jj,jk) ), & 139 & 2.*ABS( zt1 (ji ,jj,jk) ) ) 140 zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) ) & 141 & * MIN( ABS( zsp1(ji ,jj,jk) ), & 142 & 2.*ABS( zs1 (ji-1,jj,jk) ), & 143 & 2.*ABS( zs1 (ji ,jj,jk) ) ) 144 ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) ) & 145 & * MIN( ABS( ztp2(ji,jj ,jk) ), & 146 & 2.*ABS( zt2 (ji,jj-1,jk) ), & 147 & 2.*ABS( zt2 (ji,jj ,jk) ) ) 148 zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) ) & 149 & * MIN( ABS( zsp2(ji,jj ,jk) ), & 150 & 2.*ABS( zs2 (ji,jj-1,jk) ), & 151 & 2.*ABS( zs2 (ji,jj ,jk) ) ) 123 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 124 & * MIN( ABS( zslpx(ji ,jj,jk) ), & 125 & 2.*ABS( zwx (ji-1,jj,jk) ), & 126 & 2.*ABS( zwx (ji ,jj,jk) ) ) 127 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) & 128 & * MIN( ABS( zslpy(ji,jj ,jk) ), & 129 & 2.*ABS( zwy (ji,jj-1,jk) ), & 130 & 2.*ABS( zwy (ji,jj ,jk) ) ) 152 131 END DO 153 132 END DO … … 172 151 zalpha = 0.5 - z0u 173 152 zu = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj) 174 zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk) 175 zzt2 = tb(ji ,jj,jk) + zu*ztp1(ji ,jj,jk) 176 zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk) 177 zzs2 = sb(ji ,jj,jk) + zu*zsp1(ji ,jj,jk) 178 zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) 179 zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) 153 zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 154 zzwy = ptb(ji ,jj,jk) + zu * zslpx(ji ,jj,jk) 155 zwx(ji,jj,jk) = zeu * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 180 156 ! 181 157 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 182 158 zalpha = 0.5 - z0v 183 159 zv = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj) 184 zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk) 185 zzt2 = tb(ji,jj ,jk) + zv*ztp2(ji,jj ,jk) 186 zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk) 187 zzs2 = sb(ji,jj ,jk) + zv*zsp2(ji,jj ,jk) 188 zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) 189 zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) 190 END DO 191 END DO 192 END DO 193 194 ! lateral boundary conditions on zt1, zt2 ; zs1, zs2 (changed sign) 195 CALL lbc_lnk( zt1, 'U', -1. ) ; CALL lbc_lnk( zs1, 'U', -1. ) 196 CALL lbc_lnk( zt2, 'V', -1. ) ; CALL lbc_lnk( zs2, 'V', -1. ) 160 zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 161 zzwy = ptb(ji,jj ,jk) + zv * zslpy(ji,jj ,jk) 162 zwy(ji,jj,jk) = zev * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 163 END DO 164 END DO 165 END DO 166 167 !!gm bug? there is too many lbc: this have to be checked 168 ! lateral boundary conditions on zwx, zwy (changed sign) 169 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) 197 170 198 171 ! Tracer flux divergence at t-point added to the general trend … … 201 174 DO ji = fs_2, fs_jpim1 ! vector opt. 202 175 #if defined key_zco 203 zbtr = 1. / ( e1t(ji,jj) *e2t(ji,jj) )176 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 204 177 #else 205 zbtr = 1. / ( e1t(ji,jj) *e2t(ji,jj)*fse3t(ji,jj,jk) )178 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 206 179 #endif 207 180 ! horizontal advective trends 208 zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj ,jk ) & 209 & + zt2(ji,jj,jk) - zt2(ji ,jj-1,jk ) ) 210 zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj ,jk ) & 211 & + zs2(ji,jj,jk) - zs2(ji ,jj-1,jk ) ) 181 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 182 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 212 183 ! add it to the general tracer trends 213 ta(ji,jj,jk) = ta(ji,jj,jk) + zta 214 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 184 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 215 185 END DO 216 186 END DO 217 187 END DO 218 188 219 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl had - Ta: ', mask1=tmask , & 220 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 189 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 221 190 222 191 ! Save the horizontal advective trends for diagnostics 223 192 IF( l_trdtra ) THEN 224 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 225 ! 226 ! T/S ZONAL advection trends 227 DO jk = 1, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 231 ! N.B. This computation is not valid along OBCs (if any) 232 #if defined key_zco 233 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 234 z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & 235 & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 236 #else 237 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 238 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 239 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 240 #endif 241 ztrdt(ji,jj,jk) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x 242 ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x 243 END DO 244 END DO 245 END DO 246 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 247 248 ! T/S MERIDIONAL advection trends 249 DO jk = 1, jpkm1 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 253 ! N.B. This computation is not valid along OBCs (if any) 254 #if defined key_zco 255 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 256 z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) & 257 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 258 #else 259 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 260 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 261 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 262 #endif 263 ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y 264 ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y 265 END DO 266 END DO 267 END DO 268 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 269 270 ! Save the up-to-date ta and sa trends 271 ztrdt(:,:,:) = ta(:,:,:) 272 ztrds(:,:,:) = sa(:,:,:) 273 ! 274 ENDIF 275 276 ! "zonal" mean advective heat and salt transport 277 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 193 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb ) 194 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb ) 195 ENDIF 196 197 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 278 198 IF( lk_zco ) THEN 279 199 DO jk = 1, jpkm1 280 200 DO jj = 2, jpjm1 281 201 DO ji = fs_2, fs_jpim1 ! vector opt. 282 zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk) 283 zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk) 202 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 284 203 END DO 285 204 END DO 286 205 END DO 287 206 ENDIF 288 pht_adv(:) = ptr_vj( zt2(:,:,:) ) 289 pst_adv(:) = ptr_vj( zs2(:,:,:) ) 290 ENDIF 207 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 208 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 209 ENDIF 210 291 211 292 212 ! II. Vertical advective fluxes … … 296 216 ! interior values 297 217 DO jk = 2, jpkm1 298 zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) ) 299 zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) ) 218 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) ) 300 219 END DO 301 220 ! surface & bottom boundary conditions 302 zt1 (:,:, 1 ) = 0.e0 ; zt1 (:,:,jpk) = 0.e0 303 zs1 (:,:, 1 ) = 0.e0 ; zs1 (:,:,jpk) = 0.e0 221 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 304 222 305 223 ! Slopes … … 307 225 DO jj = 1, jpj 308 226 DO ji = 1, jpi 309 ztp1(ji,jj,jk) = ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) ) & 310 & * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) ) 311 zsp1(ji,jj,jk) = ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) ) & 312 & * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) ) 227 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) & 228 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 313 229 END DO 314 230 END DO … … 316 232 317 233 ! Slopes limitation 318 ! interior values 319 DO jk = 2, jpkm1 234 DO jk = 2, jpkm1 ! interior values 320 235 DO jj = 1, jpj 321 236 DO ji = 1, jpi 322 ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) ) & 323 & * MIN( ABS( ztp1(ji,jj,jk ) ), & 324 & 2.*ABS( zt1 (ji,jj,jk+1) ), & 325 & 2.*ABS( zt1 (ji,jj,jk ) ) ) 326 zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) ) & 327 & * MIN( ABS( zsp1(ji,jj,jk ) ), & 328 & 2.*ABS( zs1 (ji,jj,jk+1) ), & 329 & 2.*ABS( zs1 (ji,jj,jk ) ) ) 330 END DO 331 END DO 332 END DO 333 ! surface values 334 ztp1(:,:,1) = 0.e0 335 zsp1(:,:,1) = 0.e0 237 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 238 & * MIN( ABS( zslpx(ji,jj,jk ) ), & 239 & 2.*ABS( zwx (ji,jj,jk+1) ), & 240 & 2.*ABS( zwx (ji,jj,jk ) ) ) 241 END DO 242 END DO 243 END DO 244 zslpx(:,:,1) = 0.e0 ! surface values 336 245 337 246 ! vertical advective flux 338 ! interior values 339 DO jk = 1, jpkm1 247 DO jk = 1, jpkm1 ! interior values 340 248 zstep = z2 * rdttra(jk) 341 249 DO jj = 2, jpjm1 … … 344 252 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 345 253 zalpha = 0.5 + z0w 346 zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1) 347 zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1) 348 zzt2 = tb(ji,jj,jk ) + zw*ztp1(ji,jj,jk ) 349 zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1) 350 zzs2 = sb(ji,jj,jk ) + zw*zsp1(ji,jj,jk ) 351 zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 ) 352 zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 ) 353 END DO 354 END DO 355 END DO 356 ! surface values 254 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1) 255 zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 256 zzwy = ptb(ji,jj,jk ) + zw * zslpx(ji,jj,jk ) 257 zwx(ji,jj,jk+1) = zew * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 258 END DO 259 END DO 260 END DO 261 ! ! surface values 357 262 IF( lk_dynspg_rl .OR. lk_vvl) THEN ! rigid lid or variable volume: flux set to zero 358 z t1(:,:, 1 ) = 0.e0359 z s1(:,:, 1 ) = 0.e0263 zwx(:,:, 1 ) = 0.e0 ! surface 264 zwx(:,:,jpk) = 0.e0 ! bottom 360 265 ELSE ! free surface 361 zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1) 362 zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1) 363 ENDIF 364 365 ! bottom values 366 zt1(:,:,jpk) = 0.e0 367 zs1(:,:,jpk) = 0.e0 368 266 zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! Surface 267 zwx(:,:,jpk) = 0.e0 ! bottom 268 269 ENDIF 369 270 370 271 ! Compute & add the vertical advective trend 371 372 272 DO jk = 1, jpkm1 373 273 DO jj = 2, jpjm1 … … 375 275 zbtr = 1. / fse3t(ji,jj,jk) 376 276 ! horizontal advective trends 377 zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) ) 378 zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) ) 277 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 379 278 ! add it to the general tracer trends 380 ta(ji,jj,jk) = ta(ji,jj,jk) + zta 381 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 279 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 382 280 END DO 383 281 END DO … … 386 284 ! Save the vertical advective trends for diagnostic 387 285 ! ------------------------------------------------- 388 IF( l_trdtra ) THEN 389 ! Recompute the vertical advection zta & zsa trends computed 390 ! at the step 2. above in making the difference between the new 391 ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract 392 ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 393 394 DO jk = 1, jpkm1 395 DO jj = 2, jpjm1 396 DO ji = fs_2, fs_jpim1 ! vector opt. 397 #if defined key_zco 398 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 399 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 400 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 401 #else 402 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 403 z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk) 404 z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk) 405 #endif 406 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 407 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 408 ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 409 END DO 410 END DO 411 END DO 412 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) 413 ! 414 ENDIF 415 416 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl zad - Ta: ', mask1=tmask , & 417 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 286 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb ) 287 288 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype ) 418 289 ! 419 290 END SUBROUTINE tra_adv_muscl
Note: See TracChangeset
for help on using the changeset viewer.