Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
- Property svn:eol-style deleted
r1528 r2528 2 2 !!============================================================================== 3 3 !! *** MODULE traadv_muscl2 *** 4 !! Ocean activetracers: horizontal & vertical advective trend4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 9.0 ! 02-06 (G. Madec) from traadv_muscl 6 !! History : 1.0 ! 2002-06 (G. Madec) from traadv_muscl 7 !! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 7 8 !!---------------------------------------------------------------------- 8 9 … … 13 14 USE oce ! ocean dynamics and active tracers 14 15 USE dom_oce ! ocean space and time domain 15 USE trdmod ! ocean activetracers trends16 USE trd mod_oce ! ocean variables trends16 USE trdmod_oce ! tracers trends 17 USE trdtra ! tracers trends 17 18 USE in_out_manager ! I/O manager 18 19 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 19 20 USE trabbl ! tracers: bottom boundary layer 20 USE lib_mpp 21 USE lib_mpp ! distribued memory computing 21 22 USE lbclnk ! ocean lateral boundary condition (or mpp link) 22 23 USE diaptr ! poleward transport diagnostics 23 USE prtctl ! Print control 24 USE trc_oce ! share passive tracers/Ocean variables 25 24 26 25 27 IMPLICIT NONE 26 28 PRIVATE 27 29 28 !! * Accessibility 29 PUBLIC tra_adv_muscl2 ! routine called by step.F90 30 PUBLIC tra_adv_muscl2 ! routine called by step.F90 31 32 LOGICAL :: l_trd ! flag to compute trends 30 33 31 34 !! * Substitutions … … 33 36 # include "vectopt_loop_substitute.h90" 34 37 !!---------------------------------------------------------------------- 35 !! OPA 9.0 , LOCEAN-IPSL (2006)38 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 36 39 !! $Id$ 37 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 !!---------------------------------------------------------------------- 39 40 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 41 !!---------------------------------------------------------------------- 40 42 CONTAINS 41 43 42 SUBROUTINE tra_adv_muscl2( kt, pun, pvn, pwn ) 44 SUBROUTINE tra_adv_muscl2( kt, cdtype, p2dt, pun, pvn, pwn, & 45 & ptb, ptn, pta, kjpt ) 43 46 !!---------------------------------------------------------------------- 44 47 !! *** ROUTINE tra_adv_muscl2 *** … … 50 53 !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries 51 54 !! 52 !! ** Action : - update ( ta,sa) with the now advective tracer trends53 !! - save trends in (ztrdt,ztrds) ('key_trdtra')55 !! ** Action : - update (pta) with the now advective tracer trends 56 !! - save trends 54 57 !! 55 58 !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 56 59 !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 57 60 !!---------------------------------------------------------------------- 58 USE oce , ztrdt => ua ! use ua as workspace 59 USE oce , ztrds => va ! use va as workspace 60 !! 61 INTEGER , INTENT(in) :: kt ! ocean time-step index 62 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component 63 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component 64 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component 65 !! 66 INTEGER :: ji, jj, jk ! dummy loop indices 67 REAL(wp) :: & 68 zu, zv, zw, zeu, zev, & 69 zew, zbtr, zstep, & 70 z0u, z0v, z0w, & 71 zzt1, zzt2, zalpha, & 72 zzs1, zzs2, z2, & 73 zta, zsa, & 74 z_hdivn_x, z_hdivn_y, z_hdivn 75 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zt1, zt2, ztp1, ztp2 ! 3D workspace 76 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zs1, zs2, zsp1, zsp2 ! " " 61 USE oce , zwx => ua ! use ua as workspace 62 USE oce , zwy => va ! use va as workspace 63 !! 64 INTEGER , INTENT(in ) :: kt ! ocean time-step index 65 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 66 INTEGER , INTENT(in ) :: kjpt ! number of tracers 67 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 68 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 69 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before & now tracer fields 70 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 71 !! 72 INTEGER :: ji, jj, jk, jn ! dummy loop indices 73 REAL(wp) :: zu, z0u, zzwx ! local scalar 74 REAL(wp) :: zv, z0v, zzwy ! - - 75 REAL(wp) :: zw, z0w ! - - 76 REAL(wp) :: ztra, zbtr, zdt, zalpha 77 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zslpx, zslpy ! 3D workspace 77 78 !!---------------------------------------------------------------------- 78 79 79 IF( kt == nit000 .AND. lwp ) THEN 80 WRITE(numout,*) 81 WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme' 82 WRITE(numout,*) '~~~~~~~~~~~~~~~' 80 IF( kt == nit000 ) THEN 81 IF(lwp) WRITE(numout,*) 82 IF(lwp) WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype 83 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 84 ! 85 l_trd = .FALSE. 86 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 83 87 ENDIF 84 88 85 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 86 ELSE ; z2 = 2. 87 ENDIF 88 89 ! I. Horizontal advective fluxes 90 ! ------------------------------ 91 92 ! first guess of the slopes 93 ! interior values 94 DO jk = 1, jpkm1 95 DO jj = 1, jpjm1 96 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) ) 101 END DO 102 END DO 89 ! ! =========== 90 DO jn = 1, kjpt ! tracer loop 91 ! ! =========== 92 ! I. Horizontal advective fluxes 93 ! ------------------------------ 94 ! first guess of the slopes 95 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values 96 ! interior values 97 DO jk = 1, jpkm1 98 DO jj = 1, jpjm1 99 DO ji = 1, fs_jpim1 ! vector opt. 100 zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) ) 101 zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 102 END DO 103 END DO 104 END DO 105 ! 106 CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign) 107 CALL lbc_lnk( zwy, 'V', -1. ) 108 ! !-- Slopes of tracer 109 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values 110 DO jk = 1, jpkm1 ! interior values 111 DO jj = 2, jpj 112 DO ji = fs_2, jpi ! vector opt. 113 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) & 114 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) ) 115 zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) & 116 & * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) ) 117 END DO 118 END DO 119 END DO 120 ! 121 DO jk = 1, jpkm1 ! Slopes limitation 122 DO jj = 2, jpj 123 DO ji = fs_2, jpi ! vector opt. 124 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * 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) ) * MIN( ABS( zslpy(ji,jj ,jk) ), & 128 & 2.*ABS( zwy (ji,jj-1,jk) ), & 129 & 2.*ABS( zwy (ji,jj ,jk) ) ) 130 END DO 131 END DO 132 END DO ! interior values 133 134 ! !-- MUSCL horizontal advective fluxes 135 DO jk = 1, jpkm1 ! interior values 136 zdt = p2dt(jk) 137 DO jj = 2, jpjm1 138 DO ji = fs_2, fs_jpim1 ! vector opt. 139 ! MUSCL fluxes 140 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 141 zalpha = 0.5 - z0u 142 zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 143 zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk) 144 zzwy = ptb(ji ,jj,jk,jn) + zu * zslpx(ji ,jj,jk) 145 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 146 ! 147 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 148 zalpha = 0.5 - z0v 149 zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 150 zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk) 151 zzwy = ptb(ji,jj ,jk,jn) + zv * zslpy(ji,jj ,jk) 152 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 153 END DO 154 END DO 155 END DO 156 157 !! centered scheme at lateral b.C. if off-shore velocity 158 DO jk = 1, jpkm1 159 DO jj = 2, jpjm1 160 DO ji = fs_2, fs_jpim1 ! vector opt. 161 IF( umask(ji,jj,jk) == 0. ) THEN 162 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN 163 zwx(ji+1,jj,jk) = 0.5 * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk,jn) + ptn(ji+2,jj,jk,jn) ) 164 ENDIF 165 IF( pun(ji-1,jj,jk) < 0. ) THEN 166 zwx(ji-1,jj,jk) = 0.5 * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk,jn) + ptn(ji,jj,jk,jn) ) 167 ENDIF 168 ENDIF 169 IF( vmask(ji,jj,jk) == 0. ) THEN 170 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN 171 zwy(ji,jj+1,jk) = 0.5 * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk,jn) + ptn(ji,jj+2,jk,jn) ) 172 ENDIF 173 IF( pvn(ji,jj-1,jk) < 0. ) THEN 174 zwy(ji,jj-1,jk) = 0.5 * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk,jn) + ptn(ji,jj,jk,jn) ) 175 ENDIF 176 ENDIF 177 END DO 178 END DO 179 END DO 180 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! lateral boundary condition (changed sign) 181 182 ! Tracer flux divergence at t-point added to the general trend 183 DO jk = 1, jpkm1 184 DO jj = 2, jpjm1 185 DO ji = fs_2, fs_jpim1 ! vector opt. 186 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 187 ! horizontal advective trends 188 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 189 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 190 ! added to the general tracer trends 191 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 192 END DO 193 END DO 194 END DO 195 ! ! trend diagnostics (contribution of upstream fluxes) 196 IF( l_trd ) THEN 197 CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) ) 198 CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) ) 199 END IF 200 201 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 202 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 203 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 204 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 205 ENDIF 206 207 ! II. Vertical advective fluxes 208 ! ----------------------------- 209 ! !-- first guess of the slopes 210 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions 211 DO jk = 2, jpkm1 ! interior values 212 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 213 END DO 214 215 ! !-- Slopes of tracer 216 zslpx(:,:,1) = 0.e0 ! surface values 217 DO jk = 2, jpkm1 ! interior value 218 DO jj = 1, jpj 219 DO ji = 1, jpi 220 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) & 221 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 222 END DO 223 END DO 224 END DO 225 ! !-- Slopes limitation 226 DO jk = 2, jpkm1 ! interior values 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), & 230 & 2.*ABS( zwx (ji,jj,jk+1) ), & 231 & 2.*ABS( zwx (ji,jj,jk ) ) ) 232 END DO 233 END DO 234 END DO 235 ! !-- vertical advective flux 236 ! ! surface values (bottom already set to zero) 237 IF( lk_vvl ) THEN ; zwx(:,:, 1 ) = 0.e0 ! variable volume 238 ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 239 ENDIF 240 ! 241 DO jk = 1, jpkm1 ! interior values 242 zdt = p2dt(jk) 243 DO jj = 2, jpjm1 244 DO ji = fs_2, fs_jpim1 ! vector opt. 245 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 246 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 247 zalpha = 0.5 + z0w 248 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 249 zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1) 250 zzwy = ptb(ji,jj,jk ,jn) + zw * zslpx(ji,jj,jk ) 251 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 252 END DO 253 END DO 254 END DO 255 ! 256 DO jk = 2, jpkm1 ! centered near the bottom 257 DO jj = 2, jpjm1 258 DO ji = fs_2, fs_jpim1 ! vector opt. 259 IF( tmask(ji,jj,jk+1) == 0. ) THEN 260 IF( pwn(ji,jj,jk) > 0. ) THEN 261 zwx(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 262 ENDIF 263 ENDIF 264 END DO 265 END DO 266 END DO 267 ! 268 DO jk = 1, jpkm1 ! Compute & add the vertical advective trend 269 DO jj = 2, jpjm1 270 DO ji = fs_2, fs_jpim1 ! vector opt. 271 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 272 ! vertical advective trends 273 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 274 ! added to the general tracer trends 275 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 276 END DO 277 END DO 278 END DO 279 ! ! trend diagnostics (contribution of upstream fluxes) 280 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) ) 281 ! 103 282 END DO 104 ! bottom values105 zt1(:,:,jpk) = 0.e0 ; zt2(:,:,jpk) = 0.e0106 zs1(:,:,jpk) = 0.e0 ; zs2(:,:,jpk) = 0.e0107 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. )111 112 ! Slopes113 ! interior values114 DO jk = 1, jpkm1115 DO jj = 2, jpj116 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) ) )125 END DO126 END DO127 END DO128 ! bottom values129 ztp1(:,:,jpk) = 0.e0 ; ztp2(:,:,jpk) = 0.e0130 zsp1(:,:,jpk) = 0.e0 ; zsp2(:,:,jpk) = 0.e0131 132 ! Slopes limitation133 DO jk = 1, jpkm1134 DO jj = 2, jpj135 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) ) )152 END DO153 END DO154 END DO155 156 ! Advection terms157 ! interior values158 DO jk = 1, jpkm1159 zstep = z2 * rdttra(jk)160 DO jj = 2, jpjm1161 DO ji = fs_2, fs_jpim1 ! vector opt.162 ! volume fluxes163 #if defined key_zco164 zeu = e2u(ji,jj) * pun(ji,jj,jk)165 zev = e1v(ji,jj) * pvn(ji,jj,jk)166 #else167 zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)168 zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)169 #endif170 ! MUSCL fluxes171 z0u = SIGN( 0.5, pun(ji,jj,jk) )172 zalpha = 0.5 - z0u173 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 )180 !181 z0v = SIGN( 0.5, pvn(ji,jj,jk) )182 zalpha = 0.5 - z0v183 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 DO191 END DO192 END DO193 194 !!!! centered scheme at lateral b.C. if off-shore velocity195 DO jk = 1, jpkm1196 DO jj = 2, jpjm1197 DO ji = fs_2, fs_jpim1 ! vector opt.198 #if defined key_zco199 IF( umask(ji,jj,jk) == 0. ) THEN200 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN201 zt1(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5202 zs1(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5203 ENDIF204 IF( pun(ji-1,jj,jk) < 0. ) THEN205 zt1(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji ,jj,jk) ) * 0.5206 zs1(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji ,jj,jk) ) * 0.5207 ENDIF208 ENDIF209 IF( vmask(ji,jj,jk) == 0. ) THEN210 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN211 zt2(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5212 zs2(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5213 ENDIF214 IF( pvn(ji,jj-1,jk) < 0. ) THEN215 zt2(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji ,jj,jk) ) * 0.5216 zs2(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji ,jj,jk) ) * 0.5217 ENDIF218 ENDIF219 #else220 IF( umask(ji,jj,jk) == 0. ) THEN221 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN222 zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk) &223 & * pun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5224 zs1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk) &225 & * pun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5226 ENDIF227 IF( pun(ji-1,jj,jk) < 0. ) THEN228 zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk) &229 & * pun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji ,jj,jk) ) * 0.5230 zs1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk) &231 & * pun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji ,jj,jk) ) * 0.5232 ENDIF233 ENDIF234 IF( vmask(ji,jj,jk) == 0. ) THEN235 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN236 zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk) &237 & * pvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5238 zs2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk) &239 & * pvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5240 ENDIF241 IF( pvn(ji,jj-1,jk) < 0. ) THEN242 zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk) &243 & * pvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji ,jj,jk) ) * 0.5244 zs2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk) &245 & * pvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji ,jj,jk) ) * 0.5246 ENDIF247 ENDIF248 #endif249 END DO250 END DO251 END DO252 253 ! lateral boundary conditions on zt1, zt2 ; zs1, zs2 (changed sign)254 CALL lbc_lnk( zt1, 'U', -1. ) ; CALL lbc_lnk( zs1, 'U', -1. )255 CALL lbc_lnk( zt2, 'V', -1. ) ; CALL lbc_lnk( zs2, 'V', -1. )256 257 ! Compute & add the horizontal advective trend258 259 DO jk = 1, jpkm1260 DO jj = 2, jpjm1261 DO ji = fs_2, fs_jpim1 ! vector opt.262 #if defined key_zco263 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )264 #else265 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )266 #endif267 ! horizontal advective trends268 zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj ,jk ) &269 & + zt2(ji,jj,jk) - zt2(ji ,jj-1,jk ) )270 zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj ,jk ) &271 & + zs2(ji,jj,jk) - zs2(ji ,jj-1,jk ) )272 ! add it to the general tracer trends273 ta(ji,jj,jk) = ta(ji,jj,jk) + zta274 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa275 END DO276 END DO277 END DO278 279 ! Save the horizontal advective trends for diagnostic280 IF( l_trdtra ) THEN281 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0282 !283 ! T/S ZONAL advection trends284 DO jk = 1, jpkm1285 DO jj = 2, jpjm1286 DO ji = fs_2, fs_jpim1 ! vector opt.287 !-- Compute zonal divergence by splitting hdivn (see divcur.F90)288 ! N.B. This computation is not valid along OBCs (if any)289 #if defined key_zco290 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )291 z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) &292 & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr293 #else294 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )295 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) &296 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr297 #endif298 ztrdt(ji,jj,jk) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x299 ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x300 END DO301 END DO302 END DO303 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)304 305 ! T/S MERIDIONAL advection trends306 DO jk = 1, jpkm1307 DO jj = 2, jpjm1308 DO ji = fs_2, fs_jpim1 ! vector opt.309 !-- Compute merid. divergence by splitting hdivn (see divcur.F90)310 ! N.B. This computation is not valid along OBCs (if any)311 #if defined key_zco312 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )313 z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) &314 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr315 #else316 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )317 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) &318 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr319 #endif320 ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y321 ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y322 END DO323 END DO324 END DO325 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)326 327 ! Save the up-to-date ta and sa trends328 ztrdt(:,:,:) = ta(:,:,:)329 ztrds(:,:,:) = sa(:,:,:)330 !331 ENDIF332 333 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl2 had - Ta: ', mask1=tmask, &334 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra')335 336 ! "zonal" mean advective heat and salt transport337 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN338 IF( lk_zco ) THEN339 DO jk = 1, jpkm1340 DO jj = 2, jpjm1341 DO ji = fs_2, fs_jpim1 ! vector opt.342 zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)343 zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)344 END DO345 END DO346 END DO347 ENDIF348 pht_adv(:) = ptr_vj( zt2(:,:,:) )349 pst_adv(:) = ptr_vj( zs2(:,:,:) )350 ENDIF351 352 ! II. Vertical advective fluxes353 ! -----------------------------354 355 ! First guess of the slope356 ! interior values357 DO jk = 2, jpkm1358 zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )359 zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )360 END DO361 ! surface & bottom boundary conditions362 zt1 (:,:, 1 ) = 0.e0 ; zt1 (:,:,jpk) = 0.e0363 zs1 (:,:, 1 ) = 0.e0 ; zs1 (:,:,jpk) = 0.e0364 365 ! Slopes366 DO jk = 2, jpkm1367 DO jj = 1, jpj368 DO ji = 1, jpi369 ztp1(ji,jj,jk) = ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) ) &370 & * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )371 zsp1(ji,jj,jk) = ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) ) &372 & * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )373 END DO374 END DO375 END DO376 377 ! Slopes limitation378 ! interior values379 DO jk = 2, jpkm1380 DO jj = 1, jpj381 DO ji = 1, jpi382 ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) ) &383 & * MIN( ABS( ztp1(ji,jj,jk ) ), &384 & 2.*ABS( zt1 (ji,jj,jk+1) ), &385 & 2.*ABS( zt1 (ji,jj,jk ) ) )386 zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) ) &387 & * MIN( ABS( zsp1(ji,jj,jk ) ), &388 & 2.*ABS( zs1 (ji,jj,jk+1) ), &389 & 2.*ABS( zs1 (ji,jj,jk ) ) )390 END DO391 END DO392 END DO393 ! surface values394 ztp1(:,:,1) = 0.e0395 zsp1(:,:,1) = 0.e0396 397 ! vertical advective flux398 ! interior values399 DO jk = 1, jpkm1400 zstep = z2 * rdttra(jk)401 DO jj = 2, jpjm1402 DO ji = fs_2, fs_jpim1 ! vector opt.403 zew = pwn(ji,jj,jk+1)404 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )405 zalpha = 0.5 + z0w406 zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)407 zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)408 zzt2 = tb(ji,jj,jk ) + zw*ztp1(ji,jj,jk )409 zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)410 zzs2 = sb(ji,jj,jk ) + zw*zsp1(ji,jj,jk )411 zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )412 zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )413 END DO414 END DO415 END DO416 DO jk = 2, jpkm1417 DO jj = 2, jpjm1418 DO ji = fs_2, fs_jpim1 ! vector opt.419 IF( tmask(ji,jj,jk+1) == 0. ) THEN420 IF( pwn(ji,jj,jk) > 0. ) THEN421 zt1(ji,jj,jk) = pwn(ji,jj,jk) * ( tb(ji,jj,jk-1) + tb(ji,jj,jk) ) * 0.5422 zs1(ji,jj,jk) = pwn(ji,jj,jk) * ( sb(ji,jj,jk-1) + sb(ji,jj,jk) ) * 0.5423 ENDIF424 ENDIF425 END DO426 END DO427 END DO428 429 ! surface values430 IF( lk_vvl ) THEN431 ! variable volume : flux set to zero432 zt1(:,:, 1 ) = 0.e0433 zs1(:,:, 1 ) = 0.e0434 ELSE435 ! free surface-constant volume436 zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1)437 zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1)438 ENDIF439 440 ! bottom values441 zt1(:,:,jpk) = 0.e0442 zs1(:,:,jpk) = 0.e0443 444 445 ! Compute & add the vertical advective trend446 447 DO jk = 1, jpkm1448 DO jj = 2, jpjm1449 DO ji = fs_2, fs_jpim1 ! vector opt.450 zbtr = 1. / fse3t(ji,jj,jk)451 ! horizontal advective trends452 zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )453 zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )454 ! add it to the general tracer trends455 ta(ji,jj,jk) = ta(ji,jj,jk) + zta456 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa457 END DO458 END DO459 END DO460 461 ! Save the vertical advective trends for diagnostic462 IF( l_trdtra ) THEN463 ! Recompute the vertical advection zta & zsa trends computed464 ! at the step 2. above in making the difference between the new465 ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract466 ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends467 468 DO jk = 1, jpkm1469 DO jj = 2, jpjm1470 DO ji = fs_2, fs_jpim1 ! vector opt.471 #if defined key_zco472 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )473 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)474 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)475 #else476 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )477 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)478 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)479 #endif480 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr481 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn482 ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn483 END DO484 END DO485 END DO486 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)487 !488 ENDIF489 490 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl2 zad - Ta: ', mask1=tmask, &491 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )492 283 ! 493 284 END SUBROUTINE tra_adv_muscl2
Note: See TracChangeset
for help on using the changeset viewer.
