Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.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_tvd.F90
- Property svn:eol-style deleted
r1877 r2528 2 2 !!============================================================================== 3 3 !! *** MODULE traadv_tvd *** 4 !! Ocean activetracers: horizontal & vertical advective trend4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : !95-12 (L. Mortier) Original code7 !! ! 00-01 (H. Loukos) adapted to ORCA8 !! ! 00-10 (MA Foujols E.Kestenare) include file not routine9 !! ! 00-12 (E. Kestenare M. Levy) fix bug in trtrd indexes10 !! ! 01-07 (E. Durand G. Madec) adaptation to ORCA config11 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module12 !! 9.0 !04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl13 !! 9.0 !08-04 (S. Cravatte) add the i-, j- & k- trends computation14 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization15 !! ----------------------------------------------------------------------16 6 !! History : OPA ! 1995-12 (L. Mortier) Original code 7 !! ! 2000-01 (H. Loukos) adapted to ORCA 8 !! ! 2000-10 (MA Foujols E.Kestenare) include file not routine 9 !! ! 2000-12 (E. Kestenare M. Levy) fix bug in trtrd indexes 10 !! ! 2001-07 (E. Durand G. Madec) adaptation to ORCA config 11 !! 8.5 ! 2002-06 (G. Madec) F90: Free form and module 12 !! NEMO 1.0 ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 13 !! 2.0 ! 2008-04 (S. Cravatte) add the i-, j- & k- trends computation 14 !! - ! 2009-11 (V. Garnier) Surface pressure gradient organization 15 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 16 !!---------------------------------------------------------------------- 17 17 18 18 !!---------------------------------------------------------------------- … … 24 24 USE oce ! ocean dynamics and active tracers 25 25 USE dom_oce ! ocean space and time domain 26 USE trdmod ! ocean active tracers trends27 USE trd mod_oce ! ocean variables trends26 USE trdmod_oce ! tracers trends 27 USE trdtra ! tracers trends 28 28 USE in_out_manager ! I/O manager 29 29 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 30 USE trabbl ! Advective term of BBL31 30 USE lib_mpp 32 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 33 32 USE diaptr ! poleward transport diagnostics 34 USE prtctl ! Print control33 USE trc_oce ! share passive tracers/Ocean variables 35 34 36 35 … … 39 38 40 39 PUBLIC tra_adv_tvd ! routine called by step.F90 40 41 LOGICAL :: l_trd ! flag to compute trends 41 42 42 43 !! * Substitutions … … 44 45 # include "vectopt_loop_substitute.h90" 45 46 !!---------------------------------------------------------------------- 46 !! OPA 9.0 , LOCEAN-IPSL (2006)47 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 47 48 !! $Id$ 48 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 49 !!---------------------------------------------------------------------- 50 49 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 50 !!---------------------------------------------------------------------- 51 51 CONTAINS 52 52 53 SUBROUTINE tra_adv_tvd( kt, pun, pvn, pwn ) 53 SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn, & 54 & ptb, ptn, pta, kjpt ) 54 55 !!---------------------------------------------------------------------- 55 56 !! *** ROUTINE tra_adv_tvd *** … … 62 63 !! note: - this advection scheme needs a leap-frog time scheme 63 64 !! 64 !! ** Action : - update (ta,sa) with the now advective tracer trends 65 !! - save the trends in (ztrdt,ztrds) ('key_trdtra') 66 !!---------------------------------------------------------------------- 67 USE oce , ztrdt => ua ! use ua as workspace 68 USE oce , ztrds => va ! use va as workspace 69 !! 70 INTEGER , INTENT(in) :: kt ! ocean time-step index 71 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component 72 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component 73 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component 74 !! 75 INTEGER :: ji, jj, jk ! dummy loop indices 76 REAL(wp) :: & ! temporary scalar 77 ztat, zsat, & ! " " 78 z_hdivn_x, z_hdivn_y, z_hdivn 79 REAL(wp) :: & 80 z2dtt, zbtr, zeu, zev, & ! temporary scalar 81 zew, z2, zbtr1, & ! temporary scalar 82 zfp_ui, zfp_vj, zfp_wk, & ! " " 83 zfm_ui, zfm_vj, zfm_wk ! " " 84 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztu, ztv, ztw ! temporary workspace 85 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zsi, zsu, zsv, zsw ! " " 86 !!---------------------------------------------------------------------- 87 88 zti(:,:,:) = 0.e0 ; zsi(:,:,:) = 0.e0 89 90 IF( kt == nit000 .AND. lwp ) THEN 91 WRITE(numout,*) 92 WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme' 93 WRITE(numout,*) '~~~~~~~~~~~' 65 !! ** Action : - update (pta) with the now advective tracer trends 66 !! - save the trends 67 !!---------------------------------------------------------------------- 68 USE oce , zwx => ua ! use ua as workspace 69 USE oce , zwy => va ! use va as workspace 70 !! 71 INTEGER , INTENT(in ) :: kt ! ocean time-step index 72 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 73 INTEGER , INTENT(in ) :: kjpt ! number of tracers 74 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 75 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 76 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 77 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 78 !! 79 INTEGER :: ji, jj, jk, jn ! dummy loop indices 80 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 81 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 82 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 83 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwi, zwz ! 3D workspace 84 REAL(wp), DIMENSION (:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz 85 !!---------------------------------------------------------------------- 86 87 IF( kt == nit000 ) THEN 88 IF(lwp) WRITE(numout,*) 89 IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype 90 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 91 ! 92 l_trd = .FALSE. 93 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 94 94 ENDIF 95 96 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 97 ELSE ; z2 = 2. 98 ENDIF 99 100 ! 1. Bottom value : flux set to zero 101 ! --------------- 102 ztu(:,:,jpk) = 0.e0 ; zsu(:,:,jpk) = 0.e0 103 ztv(:,:,jpk) = 0.e0 ; zsv(:,:,jpk) = 0.e0 104 ztw(:,:,jpk) = 0.e0 ; zsw(:,:,jpk) = 0.e0 105 zti(:,:,jpk) = 0.e0 ; zsi(:,:,jpk) = 0.e0 106 107 108 ! 2. upstream advection with initial mass fluxes & intermediate update 109 ! -------------------------------------------------------------------- 110 ! upstream tracer flux in the i and j direction 111 DO jk = 1, jpkm1 112 DO jj = 1, jpjm1 113 DO ji = 1, fs_jpim1 ! vector opt. 114 zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 115 zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 116 ! upstream scheme 117 zfp_ui = zeu + ABS( zeu ) 118 zfm_ui = zeu - ABS( zeu ) 119 zfp_vj = zev + ABS( zev ) 120 zfm_vj = zev - ABS( zev ) 121 ztu(ji,jj,jk) = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj ,jk) 122 ztv(ji,jj,jk) = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji ,jj+1,jk) 123 zsu(ji,jj,jk) = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj ,jk) 124 zsv(ji,jj,jk) = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji ,jj+1,jk) 125 END DO 126 END DO 127 END DO 128 129 ! upstream tracer flux in the k direction 130 ! Surface value 131 IF( lk_vvl ) THEN 132 ! variable volume : flux set to zero 133 ztw(:,:,1) = 0.e0 134 zsw(:,:,1) = 0.e0 135 ELSE 136 ! free surface-constant volume 137 DO jj = 1, jpj 138 DO ji = 1, jpi 139 zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1) 140 ztw(ji,jj,1) = zew * tb(ji,jj,1) 141 zsw(ji,jj,1) = zew * sb(ji,jj,1) 142 END DO 143 END DO 144 ENDIF 145 146 ! Interior value 147 DO jk = 2, jpkm1 148 DO jj = 1, jpj 149 DO ji = 1, jpi 150 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 151 zfp_wk = zew + ABS( zew ) 152 zfm_wk = zew - ABS( zew ) 153 ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1) 154 zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1) 155 END DO 156 END DO 157 END DO 158 159 ! total advective trend 160 DO jk = 1, jpkm1 161 z2dtt = z2 * rdttra(jk) 162 DO jj = 2, jpjm1 163 DO ji = fs_2, fs_jpim1 ! vector opt. 164 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 165 ! total intermediate advective trends 166 ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) & 167 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) & 168 & + ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr 169 zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj ,jk ) & 170 & + zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) & 171 & + zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr 172 ! update and guess with monotonic sheme 173 ta(ji,jj,jk) = ta(ji,jj,jk) + ztat 174 sa(ji,jj,jk) = sa(ji,jj,jk) + zsat 175 zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ztat ) * tmask(ji,jj,jk) 176 zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsat ) * tmask(ji,jj,jk) 177 END DO 178 END DO 179 END DO 180 181 ! "zonal" mean advective heat and salt transport 182 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 183 pht_adv(:) = ptr_vj( ztv(:,:,:) ) 184 pst_adv(:) = ptr_vj( zsv(:,:,:) ) 185 ENDIF 186 187 ! Save the intermediate i / j / k advective trends for diagnostics 188 ! ------------------------------------------------------------------- 189 ! Warning : We should use zun instead of un in the computations below, but we 190 ! also use hdivn which is computed with un, vn (check ???). So we use un, vn 191 ! for consistency. Results are therefore approximate with key_trabbl_adv. 192 193 IF( l_trdtra ) THEN 194 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 195 ! 196 ! T/S ZONAL advection trends 95 ! 96 IF( l_trd ) THEN 97 ALLOCATE( ztrdx(jpi,jpj,jpk) ) ; ztrdx(:,:,:) = 0.e0 98 ALLOCATE( ztrdy(jpi,jpj,jpk) ) ; ztrdy(:,:,:) = 0.e0 99 ALLOCATE( ztrdz(jpi,jpj,jpk) ) ; ztrdz(:,:,:) = 0.e0 100 END IF 101 ! 102 zwi(:,:,:) = 0.e0 103 ! 104 ! ! =========== 105 DO jn = 1, kjpt ! tracer loop 106 ! ! =========== 107 ! 1. Bottom value : flux set to zero 108 ! ---------------------------------- 109 zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0 110 zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0 111 112 ! 2. upstream advection with initial mass fluxes & intermediate update 113 ! -------------------------------------------------------------------- 114 ! upstream tracer flux in the i and j direction 197 115 DO jk = 1, jpkm1 116 DO jj = 1, jpjm1 117 DO ji = 1, fs_jpim1 ! vector opt. 118 ! upstream scheme 119 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 120 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 121 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 122 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 123 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 124 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 125 END DO 126 END DO 127 END DO 128 129 ! upstream tracer flux in the k direction 130 ! Surface value 131 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable 132 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 133 ENDIF 134 ! Interior value 135 DO jk = 2, jpkm1 136 DO jj = 1, jpj 137 DO ji = 1, jpi 138 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 139 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 140 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 141 END DO 142 END DO 143 END DO 144 145 ! total advective trend 146 DO jk = 1, jpkm1 147 z2dtt = p2dt(jk) 198 148 DO jj = 2, jpjm1 199 149 DO ji = fs_2, fs_jpim1 ! vector opt. 200 150 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 201 ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr 202 ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr 203 END DO 204 END DO 205 END DO 206 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) ! save the trends 151 ! total intermediate advective trends 152 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 153 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 154 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 155 ! update and guess with monotonic sheme 156 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 157 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 158 END DO 159 END DO 160 END DO 161 ! ! Lateral boundary conditions on zwi (unchanged sign) 162 CALL lbc_lnk( zwi, 'T', 1. ) 163 164 ! ! trend diagnostics (contribution of upstream fluxes) 165 IF( l_trd ) THEN 166 ! store intermediate advective trends 167 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 168 END IF 169 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 170 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 171 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 172 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 173 ENDIF 174 175 ! 3. antidiffusive flux : high order minus low order 176 ! -------------------------------------------------- 177 ! antidiffusive flux on i and j 178 DO jk = 1, jpkm1 179 DO jj = 1, jpjm1 180 DO ji = 1, fs_jpim1 ! vector opt. 181 zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 182 zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 183 END DO 184 END DO 185 END DO 186 187 ! antidiffusive flux on k 188 zwz(:,:,1) = 0.e0 ! Surface value 207 189 ! 208 ! T/S MERIDIONAL advection trends 190 DO jk = 2, jpkm1 ! Interior value 191 DO jj = 1, jpj 192 DO ji = 1, jpi 193 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 194 END DO 195 END DO 196 END DO 197 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 198 CALL lbc_lnk( zwz, 'W', 1. ) 199 200 ! 4. monotonicity algorithm 201 ! ------------------------- 202 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 203 204 205 ! 5. final trend with corrected fluxes 206 ! ------------------------------------ 209 207 DO jk = 1, jpkm1 210 208 DO jj = 2, jpjm1 211 DO ji = fs_2, fs_jpim1 ! vector opt. 212 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 213 ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr 214 ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr 215 END DO 216 END DO 217 END DO 218 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) ! save the trends 209 DO ji = fs_2, fs_jpim1 ! vector opt. 210 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 211 ! total advective trends 212 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 213 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 214 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 215 ! add them to the general tracer trends 216 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 217 END DO 218 END DO 219 END DO 220 221 ! ! trend diagnostics (contribution of upstream fluxes) 222 IF( l_trd ) THEN 223 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 224 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 225 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 226 227 CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) ) 228 CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 229 CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 230 END IF 231 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 232 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 233 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 234 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 235 ENDIF 219 236 ! 220 ! T/S VERTICAL advection trends 221 DO jk = 1, jpkm1 222 DO jj = 2, jpjm1 223 DO ji = fs_2, fs_jpim1 ! vector opt. 224 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 225 ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 226 ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr 227 END DO 228 END DO 229 END DO 230 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) ! save the trends 231 ! 232 ENDIF 233 234 ! Lateral boundary conditions on zti, zsi (unchanged sign) 235 CALL lbc_lnk( zti, 'T', 1. ) 236 CALL lbc_lnk( zsi, 'T', 1. ) 237 238 239 ! 3. antidiffusive flux : high order minus low order 240 ! -------------------------------------------------- 241 ! antidiffusive flux on i and j 242 DO jk = 1, jpkm1 243 DO jj = 1, jpjm1 244 DO ji = 1, fs_jpim1 ! vector opt. 245 zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 246 zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 247 ztu(ji,jj,jk) = zeu * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) - ztu(ji,jj,jk) 248 zsu(ji,jj,jk) = zeu * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) ) - zsu(ji,jj,jk) 249 ztv(ji,jj,jk) = zev * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) - ztv(ji,jj,jk) 250 zsv(ji,jj,jk) = zev * ( sn(ji,jj,jk) + sn(ji,jj+1,jk) ) - zsv(ji,jj,jk) 251 END DO 252 END DO 253 END DO 254 255 ! antidiffusive flux on k 256 ! Surface value 257 ztw(:,:,1) = 0.e0 258 zsw(:,:,1) = 0.e0 259 260 ! Interior value 261 DO jk = 2, jpkm1 262 DO jj = 1, jpj 263 DO ji = 1, jpi 264 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 265 ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 266 zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk) 267 END DO 268 END DO 269 END DO 270 271 ! Lateral bondary conditions 272 CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( zsu, 'U', -1. ) 273 CALL lbc_lnk( ztv, 'V', -1. ) ; CALL lbc_lnk( zsv, 'V', -1. ) 274 CALL lbc_lnk( ztw, 'W', 1. ) ; CALL lbc_lnk( zsw, 'W', 1. ) 275 276 ! 4. monotonicity algorithm 277 ! ------------------------- 278 CALL nonosc( tb, ztu, ztv, ztw, zti, z2 ) 279 CALL nonosc( sb, zsu, zsv, zsw, zsi, z2 ) 280 281 282 ! 5. final trend with corrected fluxes 283 ! ------------------------------------ 284 DO jk = 1, jpkm1 285 DO jj = 2, jpjm1 286 DO ji = fs_2, fs_jpim1 ! vector opt. 287 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 288 ! total advective trends 289 ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) & 290 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) & 291 & + ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr 292 zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj ,jk ) & 293 & + zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) & 294 & + zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr 295 ! add them to the general tracer trends 296 ta(ji,jj,jk) = ta(ji,jj,jk) + ztat 297 sa(ji,jj,jk) = sa(ji,jj,jk) + zsat 298 END DO 299 END DO 300 END DO 301 302 303 ! Save the advective trends for diagnostics 304 ! -------------------------------------------- 305 306 IF( l_trdtra ) THEN 307 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 308 ! 309 ! T/S ZONAL advection trends 310 DO jk = 1, jpkm1 311 DO jj = 2, jpjm1 312 DO ji = fs_2, fs_jpim1 ! vector opt. 313 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 314 ! N.B. This computation is not valid along OBCs (if any) 315 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 316 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 317 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 318 !-- Compute T/S zonal advection trends 319 ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x 320 ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x 321 END DO 322 END DO 323 END DO 324 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED 325 ! 326 ! T/S MERIDIONAL advection trends 327 DO jk = 1, jpkm1 328 DO jj = 2, jpjm1 329 DO ji = fs_2, fs_jpim1 ! vector opt. 330 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 331 ! N.B. This computation is not valid along OBCs (if any) 332 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 333 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 334 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 335 !-- Compute T/S meridional advection trends 336 ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y 337 ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y 338 END DO 339 END DO 340 END DO 341 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED 342 ! 343 ! T/S VERTICAL advection trends 344 DO jk = 1, jpkm1 345 DO jj = 2, jpjm1 346 DO ji = fs_2, fs_jpim1 ! vector opt. 347 zbtr1 = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 348 #if defined key_zco 349 zbtr = zbtr1 350 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 351 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 352 #else 353 zbtr = zbtr1 / fse3t(ji,jj,jk) 354 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) 355 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) 356 #endif 357 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 358 zbtr = zbtr1 / fse3t(ji,jj,jk) 359 ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr - tn(ji,jj,jk) * z_hdivn 360 ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr - sn(ji,jj,jk) * z_hdivn 361 END DO 362 END DO 363 END DO 364 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED 365 ! 366 ENDIF 367 368 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' tvd adv - Ta: ', mask1=tmask, & 369 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 370 371 ! "zonal" mean advective heat and salt transport 372 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 373 pht_adv(:) = ptr_vj( ztv(:,:,:) ) + pht_adv(:) 374 pst_adv(:) = ptr_vj( zsv(:,:,:) ) + pst_adv(:) 375 ENDIF 237 ENDDO 238 ! 239 IF( l_trd ) THEN 240 DEALLOCATE( ztrdx ) ; DEALLOCATE( ztrdy ) ; DEALLOCATE( ztrdz ) 241 END IF 376 242 ! 377 243 END SUBROUTINE tra_adv_tvd 378 244 379 245 380 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p rdt )246 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) 381 247 !!--------------------------------------------------------------------- 382 248 !! *** ROUTINE nonosc *** … … 391 257 !! in-space based differencing for fluid 392 258 !!---------------------------------------------------------------------- 393 REAL(wp), INTENT( in ) :: prdt ! ??? 394 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) :: & 395 pbef, & ! before field 396 paft, & ! after field 397 paa, & ! monotonic flux in the i direction 398 pbb, & ! monotonic flux in the j direction 399 pcc ! monotonic flux in the k direction 259 REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step 260 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 261 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 400 262 !! 401 263 INTEGER :: ji, jj, jk ! dummy loop indices … … 423 285 DO jk = 1, jpkm1 424 286 ikm1 = MAX(jk-1,1) 425 z2dtt = p rdt * rdttra(jk)287 z2dtt = p2dt(jk) 426 288 DO jj = 2, jpjm1 427 289 DO ji = fs_2, fs_jpim1 ! vector opt. … … 456 318 END DO 457 319 END DO 458 459 ! lateral boundary condition on zbetup & zbetdo (unchanged sign) 460 CALL lbc_lnk( zbetup, 'T', 1. ) 461 CALL lbc_lnk( zbetdo, 'T', 1. ) 320 CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 321 462 322 463 323 … … 486 346 END DO 487 347 END DO 488 489 ! lateral boundary condition on paa, pbb, pcc 490 CALL lbc_lnk( paa, 'U', -1. ) ! changed sign 491 CALL lbc_lnk( pbb, 'V', -1. ) ! changed sign 348 CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign) 492 349 ! 493 350 END SUBROUTINE nonosc
Note: See TracChangeset
for help on using the changeset viewer.
