Changeset 791 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
- Timestamp:
- 2008-01-12T21:33:34+01:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 7.0 ! 95-12 (L. Mortier) Original code7 !! 8.0 ! 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 !! NEMO 1.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl13 !! - ! 08-04 (S. Cravatte) add the i-, j- & k- trends computation14 !! - ! 05-11 (V. Garnier) Surface pressure gradient organization15 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : 7.0 ! 1995-12 (L. Mortier) Original code 7 !! 8.0 ! 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 !! - ! 2008-04 (S. Cravatte) add the i-, j- & k- trends computation 14 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 15 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 16 16 !!---------------------------------------------------------------------- 17 17 … … 42 42 !!---------------------------------------------------------------------- 43 43 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 44 !! $Id :$44 !! $Id$ 45 45 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 46 46 !!---------------------------------------------------------------------- … … 70 70 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 71 71 !! 72 INTEGER :: ji, jj, jk ! dummy loop indices 73 REAL(wp) :: ztai, ztaj, ztak 74 REAL(wp) :: z2dtt, zbtr, zeu, zev ! temporary scalar 75 REAL(wp) :: zew, z2 ! temporary scalar 76 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! " " 77 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! " " 78 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztu, ztv, ztw ! temporary workspace 72 INTEGER :: ji, jj, jk ! dummy loop indices 73 REAL(wp) :: z2dtt, zbtr, z2, zzti ! temporary scalar 74 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! " " 75 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! " " 76 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztu, ztv, ztw ! 3D workspace 79 77 !!---------------------------------------------------------------------- 80 78 … … 87 85 ENDIF 88 86 89 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 90 ELSE ; z2 = 2. 87 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ! euler time-stepping 88 ELSE ; z2 = 2. ! leap-frog time-stepping 91 89 ENDIF 92 90 … … 103 101 DO jj = 1, jpjm1 104 102 DO ji = 1, fs_jpim1 ! vector opt. 105 zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 106 zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 107 ! upstream scheme 108 zfp_ui = zeu + ABS( zeu ) 109 zfm_ui = zeu - ABS( zeu ) 110 zfp_vj = zev + ABS( zev ) 111 zfm_vj = zev - ABS( zev ) 112 ztu(ji,jj,jk) = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj ,jk) 113 ztv(ji,jj,jk) = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji ,jj+1,jk) 114 END DO 115 END DO 116 END DO 117 103 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 104 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 105 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 106 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 107 ztu(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj ,jk) ) 108 ztv(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji ,jj+1,jk) ) 109 END DO 110 END DO 111 END DO 112 ! 118 113 ! upstream tracer flux in the k direction 119 ! Surface value 120 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or variable volume: flux set to zero 121 ztw(:,:,1) = 0.e0 122 ELSE ! free surface 123 ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1) 124 ENDIF 125 126 ! Interior value 127 DO jk = 2, jpkm1 114 ! ! Surface value 115 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! rigid lid or non-linear fs 116 ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1) ! free surface 117 ENDIF 118 ! 119 DO jk = 2, jpkm1 ! Interior value 128 120 DO jj = 1, jpj 129 121 DO ji = 1, jpi 130 z ew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)131 zf p_wk = zew + ABS( zew)132 z fm_wk = zew - ABS( zew)133 ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)134 135 136 END DO137 138 ! total advective trend139 DO jk = 1, jpkm1122 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 123 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 124 ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) ) 125 END DO 126 END DO 127 END DO 128 ! 129 ! upstream advective trend 130 DO jk = 1, jpkm1 131 z2dtt = z2 * rdttra(jk) 140 132 DO jj = 2, jpjm1 141 133 DO ji = fs_2, fs_jpim1 ! vector opt. 142 134 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 143 ! i- j- horizontal & k- vertical advective trends144 ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) ) * zbtr145 ztaj = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) ) * zbtr146 ztak = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr147 135 ! total intermediate advective trends 148 zti(ji,jj,jk) = ztai + ztaj + ztak 149 END DO 150 END DO 151 END DO 152 153 ! Save the horizontal advective trends for diagnostic 154 ! ----------------------------------------------------- 136 zzti = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) & 137 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) & 138 & + ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr 139 ! update and guess with monotonic sheme 140 pta(ji,jj,jk) = pta(ji,jj,jk) + zzti 141 zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zzti ) * tmask(ji,jj,jk) 142 END DO 143 END DO 144 END DO 145 ! 146 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti (unchanged sign) 147 148 149 ! ! trend diagnostics (contribution of upstream fluxes) 155 150 IF( l_trdtra ) THEN 156 151 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn ) … … 158 153 CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn ) 159 154 ENDIF 160 161 ! update and guess with monotonic sheme 162 DO jk = 1, jpkm1 163 z2dtt = z2 * rdttra(jk) 164 DO jj = 2, jpjm1 165 DO ji = fs_2, fs_jpim1 ! vector opt. 166 pta(ji,jj,jk) = pta(ji,jj,jk) + zti(ji,jj,jk) 167 zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk) 168 END DO 169 END DO 170 END DO 171 172 ! Lateral boundary conditions on zti (unchanged sign) 173 CALL lbc_lnk( zti, 'T', 1. ) 155 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 156 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 157 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( ztv(:,:,:) ) 158 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( ztv(:,:,:) ) 159 ENDIF 174 160 175 161 176 162 ! 3. antidiffusive flux : high order minus low order 177 163 ! -------------------------------------------------- 178 ! antidiffusive flux on i and j164 ! ! anti-diffusive flux on i and j 179 165 DO jk = 1, jpkm1 180 166 DO jj = 1, jpjm1 181 167 DO ji = 1, fs_jpim1 ! vector opt. 182 zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 183 zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 184 ztu(ji,jj,jk) = zeu * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk) 185 ztv(ji,jj,jk) = zev * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk) 186 END DO 187 END DO 188 END DO 189 190 ! antidiffusive flux on k 191 ! Surface value 192 ztw(:,:,1) = 0.e0 193 194 ! Interior value 195 DO jk = 2, jpkm1 168 ztu(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk) 169 ztv(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk) 170 END DO 171 END DO 172 END DO 173 ! ! antidiffusive flux on k 174 ztw(:,:,1) = 0.e0 ! Surface value 175 ! 176 DO jk = 2, jpkm1 ! Interior value 196 177 DO jj = 1, jpj 197 178 DO ji = 1, jpi 198 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 199 ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 200 END DO 201 END DO 202 END DO 203 204 ! Lateral bondary conditions 205 CALL lbc_lnk( ztu, 'U', -1. ) 179 ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 180 END DO 181 END DO 182 END DO 183 ! 184 CALL lbc_lnk( ztu, 'U', -1. ) ! Lateral bondary conditions on upstream fluxes 206 185 CALL lbc_lnk( ztv, 'V', -1. ) 207 186 CALL lbc_lnk( ztw, 'W', 1. ) 208 187 209 ! 4. monotonicity algorithm 210 ! ------------------------- 188 ! ! monotonicity algorithm 211 189 CALL nonosc( ptb, ztu, ztv, ztw, zti, z2 ) 212 190 213 191 214 ! 5. final trend with correctedfluxes215 ! ------------------------------------ 192 ! 4. final trend with anti-diffusive fluxes 193 ! ----------------------------------------- 216 194 DO jk = 1, jpkm1 217 195 DO jj = 2, jpjm1 218 196 DO ji = fs_2, fs_jpim1 ! vector opt. 219 197 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 220 ! i- j- horizontal & k- vertical advective trends 221 ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk )) * zbtr 222 ztaj = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk )) * zbtr 223 ztak = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1)) * zbtr 224 225 ! add them to the general tracer trends 226 pta(ji,jj,jk) = pta(ji,jj,jk) + ztai + ztaj + ztak 227 END DO 228 END DO 229 END DO 230 231 !!gm the transport computation is wrong, the upstream part is missing ! 232 ! "zonal" mean advective heat and salt transport 233 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 234 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( ztv(:,:,:) ) 235 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( ztv(:,:,:) ) 236 ENDIF 237 238 ! Save the horizontal advective trends for diagnostic 239 ! ----------------------------------------------------- 240 IF( l_trdtra ) THEN 198 ! anti-diffusive trends added to the general tracer trends 199 pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) & 200 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) & 201 & + ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr 202 END DO 203 END DO 204 END DO 205 206 IF( l_trdtra ) THEN ! trend diagnostic (contribution of anti-diffusive fluxes) 241 207 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn, cnbpas='bis' ) ! <<< Add to iad trend 242 208 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn, cnbpas='bis' ) ! <<< Add to jad trend 243 209 CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn, cnbpas='bis' ) ! <<< Add to zad trend 244 210 ENDIF 245 211 ! ! "Poleward" heat and salt transports (contribution of anti-diffusive fluxes) 212 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 213 IF( ktra == jp_tem) pht_adv(:) = pht_adv(:) + ptr_vj( ztv(:,:,:) ) 214 IF( ktra == jp_sal) pst_adv(:) = pst_adv(:) + ptr_vj( ztv(:,:,:) ) 215 ENDIF 216 ! ! control print 246 217 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' tvd - adv: ', mask1=tmask, clinfo3=cdtype ) 247 218 ! … … 276 247 zbetup(:,:,:) = 0.e0 ; zbetdo(:,:,:) = 0.e0 277 248 249 !!gm optimisation : add the optimal version I wrote 1 year ago 278 250 ! Search local extrema 279 251 ! --------------------
Note: See TracChangeset
for help on using the changeset viewer.
