Changeset 791
- Timestamp:
- 2008-01-12T21:33:34+01:00 (16 years ago)
- Location:
- branches/dev_001_GM/NEMO/OPA_SRC/TRA
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv.F90
r786 r791 4 4 !! Ocean active tracers: advection trend 5 5 !!============================================================================== 6 !! History : 9.0 ! 05-11 (G. Madec) Original code 6 !! History : 1.0 ! 2005-11 (G. Madec) Original code 7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 7 8 !!---------------------------------------------------------------------- 8 9 … … 23 24 USE ldftra_oce ! lateral diffusion coefficient on tracers 24 25 USE in_out_manager ! I/O manager 25 ! USE prtctl ! Print control26 26 27 27 IMPLICIT NONE … … 59 59 !! ** Method : - Update (ua,va) with the advection term following nadv 60 60 !!---------------------------------------------------------------------- 61 #if ( defined key_trabbl_adv || defined key_traldf_eiv ) 62 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zun, zvn, zwn ! effective velocity 63 #else 64 USE oce, ONLY : zun => un ! the effective velocity is the 65 USE oce, ONLY : zvn => vn ! Eulerian velocity 66 USE oce, ONLY : zwn => wn ! 67 #endif 61 INTEGER, INTENT( in ) :: kt ! ocean time-step index 68 62 !! 69 INTEGER, INTENT( in ) :: kt ! ocean time-step index 63 INTEGER :: jk ! dummy loop index 64 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zun, zvn, zwn ! effective transports 70 65 !!---------------------------------------------------------------------- 71 66 72 67 IF( kt == nit000 ) CALL tra_adv_ctl ! initialisation & control of options 73 68 69 70 ! ! effective transport 71 DO jk = 1, jpkm1 74 72 #if defined key_trabbl_adv 75 zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:) ! add the bbl velocity 76 zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:) 77 zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:) 73 ! ! eulerian + bbl transport 74 zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * ( un(:,:,jk) - u_bbl(:,:,jk) ) 75 zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * ( vn(:,:,jk) - v_bbl(:,:,jk) ) 76 zwn(:,:,jk) = e1t(:,:) * e2t(:,:) * ( wn(:,:,jk) + w_bbl(:,:,jk) ) 77 #else 78 ! ! eulerian transport only 79 zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 80 zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 81 zwn(:,:,jk) = e1t(:,:) * e2t(:,:) * wn(:,:,jk) 78 82 #endif 79 IF( lk_traldf_eiv ) THEN ! commpute and add the eiv velocity 80 IF( .NOT. lk_trabbl_adv ) THEN 81 zun(:,:,:) = un(:,:,:) 82 zvn(:,:,:) = vn(:,:,:) 83 zwn(:,:,:) = wn(:,:,:) 84 ENDIF 85 CALL tra_adv_eiv( kt, zun, zvn, zwn ) 86 ENDIF 83 END DO 84 zwn(:,:,jpk) = 0.e0 ! no transport trough the bottom 85 86 ! ! add the eiv transport (if necessary) 87 IF( lk_traldf_eiv ) CALL tra_adv_eiv( kt, zun, zvn, zwn ) 88 87 89 88 90 SELECT CASE ( nadv ) ! compute advection trend and add it to general trend 89 CASE ( 0 ) ; CALL tra_adv_cen2 ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! 2nd order centered 90 CALL tra_adv_cen2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! 2nd order centered 91 ! CASE ( 1 ) ; CALL tra_adv_cen2_jki( kt, zun, zvn, zwn ) ! 2nd order centered scheme 91 ! 92 CASE ( 1 ) ; CALL tra_adv_cen2 ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! 2nd order centered 93 CALL tra_adv_cen2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! on T & S 94 ! 92 95 CASE ( 2 ) ; CALL tra_adv_tvd ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! TVD scheme 93 CALL tra_adv_tvd ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! TVD scheme 96 CALL tra_adv_tvd ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! on T & S 97 ! 94 98 CASE ( 3 ) ; CALL tra_adv_muscl ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb , ta ) ! MUSCL scheme 95 CALL tra_adv_muscl ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb , sa ) ! MUSCL scheme 99 CALL tra_adv_muscl ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb , sa ) ! on T & S 100 ! 96 101 CASE ( 4 ) ; CALL tra_adv_muscl2 ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! MUSCL2 scheme 97 CALL tra_adv_muscl2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! MUSCL2 scheme 102 CALL tra_adv_muscl2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! on T & S 103 ! 98 104 CASE ( 5 ) ; CALL tra_adv_ubs ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! UBS scheme 99 CALL tra_adv_ubs ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! UBS scheme 105 CALL tra_adv_ubs ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! on T & S 106 ! 100 107 CASE ( 6 ) ; CALL tra_adv_qck ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! QUICKEST scheme 101 CALL tra_adv_qck ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! QUICKEST scheme108 CALL tra_adv_qck ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! on T & S 102 109 ! 103 CASE (-1 ) ! esopa: test all possibility with control print110 CASE (-1 ) ! NEMO debug : test all the schemes at once 104 111 CALL tra_adv_cen2 ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! 2nd order centered 105 CALL tra_adv_cen2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! 2nd order centered112 CALL tra_adv_cen2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 106 113 CALL tra_adv_tvd ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! TVD scheme 107 CALL tra_adv_tvd ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! TVD scheme114 CALL tra_adv_tvd ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 108 115 CALL tra_adv_muscl ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb , ta ) ! MUSCL scheme 109 CALL tra_adv_muscl ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb , sa ) ! MUSCL scheme116 CALL tra_adv_muscl ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb , sa ) 110 117 CALL tra_adv_muscl2 ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! MUSCL2 scheme 111 CALL tra_adv_muscl2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! MUSCL2 scheme118 CALL tra_adv_muscl2 ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 112 119 CALL tra_adv_ubs ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! UBS scheme 113 CALL tra_adv_ubs ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! UBS scheme120 CALL tra_adv_ubs ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 114 121 CALL tra_adv_qck ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta ) ! QUICKEST scheme 115 CALL tra_adv_qck ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) ! QUICKEST scheme122 CALL tra_adv_qck ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 116 123 END SELECT 117 124 ! … … 127 134 !!---------------------------------------------------------------------- 128 135 INTEGER :: ioptio 129 136 !! 130 137 NAMELIST/nam_traadv/ ln_traadv_cen2 , ln_traadv_tvd, & 131 138 & ln_traadv_muscl, ln_traadv_muscl2, & … … 163 170 & CALL ctl_stop( 'cross-land advection only with 2nd order advection scheme' ) 164 171 165 ! ! Set nadv 166 IF( ln_traadv_cen2 ) nadv = 0 167 #if defined key_mpp_omp 172 ! ! Set nadv 168 173 IF( ln_traadv_cen2 ) nadv = 1 169 #endif170 174 IF( ln_traadv_tvd ) nadv = 2 171 175 IF( ln_traadv_muscl ) nadv = 3 … … 175 179 IF( lk_esopa ) nadv = -1 176 180 177 IF(lwp) THEN ! Print the choice181 IF(lwp) THEN ! Print the choice 178 182 WRITE(numout,*) 179 IF( nadv == 0 ) WRITE(numout,*) ' 2nd order scheme is used' 180 IF( nadv == 1 ) WRITE(numout,*) ' 2nd order scheme is usedi, k-j-i case' 183 IF( nadv == 1 ) WRITE(numout,*) ' 2nd order scheme is used' 181 184 IF( nadv == 2 ) WRITE(numout,*) ' TVD scheme is used' 182 185 IF( nadv == 3 ) WRITE(numout,*) ' MUSCL scheme is used' … … 185 188 IF( nadv == 6 ) WRITE(numout,*) ' PPM scheme is used' 186 189 IF( nadv == 7 ) WRITE(numout,*) ' QUICKEST scheme is used' 187 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection scheme'190 IF( nadv == -1 ) WRITE(numout,*) ' NEMO debug: test all advection scheme at once' 188 191 ENDIF 189 192 ! -
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 8.2 ! 01-08 (G. Madec, E. Durand) trahad+trazad = traadv7 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module8 !! NEMO 1.0 ! 05-11 (V. Garnier) Surface pressure gradient organization9 !! - ! 05-11 (V. Garnier) Surface pressure gradient organization10 !! - ! 06-04 (R. Benshila, G. Madec) Step reorganization11 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : 8.2 ! 2001-08 (G. Madec, E. Durand) trahad+trazad = traadv 7 !! 8.5 ! 2002-06 (G. Madec) F90: Free form and module 8 !! NEMO 1.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 9 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 10 !! - ! 2006-04 (R. Benshila, G. Madec) Step reorganization 11 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 12 12 !!---------------------------------------------------------------------- 13 13 … … 34 34 PUBLIC tra_adv_cen2 ! routine called by step.F90 35 35 36 REAL(wp), DIMENSION(jpi,jpj) :: btr2 ! inverse of T-point surface [1/(e1t*e2t)]37 38 36 !! * Substitutions 39 37 # include "domzgr_substitute.h90" … … 41 39 !!---------------------------------------------------------------------- 42 40 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 43 !! $Id :$41 !! $Id$ 44 42 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 45 43 !!---------------------------------------------------------------------- … … 62 60 !! Part 0 : compute the upstream / centered flag 63 61 !! (3D array, zind, defined at T-point (0<zind<1)) 64 !! Part I : horizontal advection 65 !! * centered flux: 66 !! zcenu = e2u*e3u un mi(ptn) 67 !! zcenv = e1v*e3v vn mj(ptn) 68 !! * upstream flux: 69 !! zupsu = e2u*e3u un (ptb(i) or ptb(i-1) ) [un>0 or <0] 70 !! zupsv = e1v*e3v vn (tb(j) or tb(j-1) ) [vn>0 or <0] 71 !! * mixed upstream / centered horizontal advection scheme 62 !! Part I : advective fluxes 63 !! * centered flux: zcenu = pun mi(ptn) (idem for v and w) 64 !! * upstream flux: zupsu = pun (ptb(i) or ptb(i-1) ) [un>0 or <0] 65 !! * mixed upstream / centered flux: 72 66 !! zcofi = max(zind(i+1), zind(i)) 73 !! zcofj = max(zind(j+1), zind(j))74 67 !! zwx = zcofi * zupsu + (1-zcofi) * zcenu 75 !! zwy = zcofj * zupsv + (1-zcofj) * zcenv68 !! Part III : advective trend 76 69 !! * horizontal advective trend (divergence of the fluxes) 77 !! zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] } 78 !! * Add this trend now to the general trend of tracer (pta): 79 !! pta = pta + zta 80 !! * trend diagnostic (lk_trdtra=T): the trend is 81 !! saved for diagnostics. The trends saved is expressed as 82 !! Uh.gradh(T), i.e. 83 !! save trend = zta + ptn divn 84 !! In addition, the advective trend in the two horizontal direc- 85 !! tion is also re-computed as Uh gradh(T). Indeed hadt+ptn divn is 86 !! equal to (in s-coordinates, and similarly in z-coord.): 87 !! zta+ptn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u un di[ptn] ) 88 !! +mj-1( e1v*e3v vn mj[ptn] ) } 89 !! NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so 90 !! they vanish from the expression of the flux and divergence. 91 !! 92 !! Part II : vertical advection 93 !! the advective trend is computed as follows : 94 !! zta = 1/e3t dk+1[ zwx ] 95 !! where the vertical advective flux, zwx, is given by : 96 !! zwx = zcofk * zupst + (1-zcofk) * zcent 97 !! with 98 !! zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0] 99 !! zcenu = centered flux = wn * mk(ptn) 100 !! The surface boundary condition is : 101 !! rigid-lid (lk_dynspg_frd = T) : zero advective flux 102 !! free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * ptn(:,:,1) 103 !! Add this trend now to the general trend of tracer (ta,sa): 104 !! (ta,sa) = (ta,sa) + ( zta , zsa ) 105 !! Trend diagnostic ('key_trdtra' defined): the trend is 106 !! saved for diagnostics. The trends saved is expressed as : 107 !! save trend = w.gradz(T) = zta - ptn divn. 108 !! 109 !! ** Action : - update (ta,sa) with the now advective tracer trends 70 !! added to the general trend of tracer (pta): 71 !! pta = pta - 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] + dk[zwz] } 72 !! * trend diagnostic (lk_trdtra=T): the advective fluxes are 73 !! where the 3 part of advective trend are recomputed as U.grad[T] 74 !! and then used in on-line diagnostics. 75 !! 76 !! NB: The surface boundary condition is no flux for both 77 !! rigid-lid (lk_dynspg_frd = T) and non-linear free-surface 78 !! but a flux (pwn(1) * ptn(1)) cross the surface in linear 79 !! free-surface case (lk_dynspg_fsc = T) 80 !! 81 !! ** Action : - update pta with the 2nd order centered advective trends 110 82 !! - trend diagnostics (lk_trdtra=T) 111 83 !!---------------------------------------------------------------------- … … 118 90 !! 119 91 INTEGER :: ji, jj, jk ! dummy loop indices 120 REAL(wp) :: zbtr , zta, zhw, ze3tr! temporary scalars121 REAL(wp) :: zcofi, z fui, zcenut, zupsut, zfp_ui, zfm_ui! " "122 REAL(wp) :: zcofj, z fvj, zcenvt, zupsvt, zfp_vj, zfm_vj! " "123 REAL(wp) :: zcofk, zcent , zupst , zfp_w , zfm_w! " "124 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, z ind ! 3D workspace92 REAL(wp) :: zbtr ! temporary scalars 93 REAL(wp) :: zcofi, zcenut, zupsut, zfp_ui, zfm_ui ! " " 94 REAL(wp) :: zcofj, zcenvt, zupsvt, zfp_vj, zfm_vj ! " " 95 REAL(wp) :: zcofk, zcent , zupst , zfp_w , zfm_w ! " " 96 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwz, zind ! 3D workspace 125 97 !!---------------------------------------------------------------------- 126 98 … … 129 101 IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme' 130 102 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 131 !132 btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )133 103 ENDIF 134 104 … … 150 120 END DO 151 121 152 ! I. Horizontal advection 153 ! ----------------------- 154 155 ! ! =============== 156 DO jk = 1, jpkm1 ! Horizontal slab 157 ! ! =============== 158 ! 159 DO jj = 1, jpjm1 ! Horizontal advective fluxes 122 ! I. Advective fluxes 123 ! ------------------- 124 DO jk = 1, jpkm1 ! Horizontal advective fluxes 125 DO jj = 1, jpjm1 160 126 DO ji = 1, fs_jpim1 ! vector opt. 161 127 ! upstream indicator 162 128 zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) 163 129 zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) 164 ! volume fluxes * 1/2165 #if defined key_zco166 zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)167 zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)168 #else169 zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)170 zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)171 #endif172 130 ! upstream scheme 173 zfp_ui = zfui + ABS( zfui)174 zf p_vj = zfvj + ABS( zfvj)175 zf m_ui = zfui - ABS( zfui)176 zfm_vj = zfvj - ABS( zfvj)131 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 132 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 133 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 134 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 177 135 zupsut = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj ,jk) 178 136 zupsvt = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji ,jj+1,jk) 179 137 ! centered scheme 180 zcenut = zfui* ( ptn(ji,jj,jk) + ptn(ji+1,jj ,jk) )181 zcenvt = zfvj* ( ptn(ji,jj,jk) + ptn(ji ,jj+1,jk) )138 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj ,jk) ) 139 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji ,jj+1,jk) ) 182 140 ! mixed centered / upstream scheme 183 zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut 184 zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt 185 END DO 186 END DO 187 188 DO jj = 2, jpjm1 ! horizontal tracer flux divergence added to the general trend 189 DO ji = fs_2, fs_jpim1 ! vector opt. 190 #if defined key_zco 191 zbtr = btr2(ji,jj) 192 #else 193 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 194 #endif 195 ! horizontal advective trends 196 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 197 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 198 ! add it to the general tracer trends 199 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 200 END DO 201 END DO 202 ! ! =============== 203 END DO ! End of slab 204 ! ! =============== 205 206 ! Save the horizontal advective trends for diagnostic 207 ! ----------------------------------------------------- 208 IF( l_trdtra ) THEN 209 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 210 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 211 ENDIF 212 213 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - had: ', mask1=tmask, clinfo3=cdtype ) 214 215 216 ! "Poleward" heat and salt transport 217 ! ---------------------------------- 218 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 219 IF( lk_zco ) THEN 220 DO jk = 1, jpkm1 221 DO jj = 2, jpjm1 222 DO ji = fs_2, fs_jpim1 ! vector opt. 223 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 224 END DO 225 END DO 226 END DO 227 ENDIF 228 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 229 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 230 ENDIF 231 232 233 ! II. Vertical advection 234 ! ---------------------- 235 236 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or non-linear free surface 237 zwx(:,:, 1 ) = 0.e0 ! Surface value : zero flux 238 zwx(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 239 ELSE ! linear free surface 240 zwx(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! Surface : : advection through z=0 241 zwx(:,:,jpk) = 0.e0 ! Bottom : flux set to zero 242 ENDIF 243 244 DO jk = 2, jpk ! Vertical advective fluxes (at w-point) 141 zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut ) 142 zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt ) 143 END DO 144 END DO 145 END DO 146 ! ! Vertical advective fluxes 147 zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 148 ! ! Surface value : 149 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! rigid lid or non-linear fs 150 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! linear free surface 151 ENDIF 152 ! 153 DO jk = 2, jpkm1 ! interior values 245 154 DO jj = 2, jpjm1 246 155 DO ji = fs_2, fs_jpim1 ! vector opt. 247 156 ! upstream indicator 248 157 zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 249 ! velocity * 1/2250 zhw = 0.5 * pwn(ji,jj,jk)251 158 ! upstream scheme 252 zfp_w = zhw + ABS( zhw)253 zfm_w = zhw - ABS( zhw)159 zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 160 zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 254 161 zupst = zfp_w * ptb(ji,jj,jk) + zfm_w * ptb(ji,jj,jk-1) 255 162 ! centered scheme 256 zcent = zhw* ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) )163 zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) 257 164 ! mixed centered / upstream scheme 258 zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 259 END DO 260 END DO 261 END DO 262 263 DO jk = 1, jpkm1 ! Tracer flux divergence at t-point added to the general trend 264 DO jj = 2, jpjm1 165 zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent ) 166 END DO 167 END DO 168 END DO 169 170 ! II. Divergence of advective fluxes 171 ! ---------------------------------- 172 DO jk = 1, jpkm1 173 DO jj = 2, jpjm1 265 174 DO ji = fs_2, fs_jpim1 ! vector opt. 266 ze3tr = 1. / fse3t(ji,jj,jk) 267 ! vertical advective trends 268 zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 269 ! add it to the general tracer trends 270 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 271 END DO 272 END DO 273 END DO 274 275 ! 3. Save the vertical advective trends for diagnostic 276 ! ---------------------------------------------------- 277 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptn ) 278 279 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - zad : ', mask1=tmask, clinfo3=cdtype ) 175 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 176 ! advective trends added to the general tracer trends 177 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 178 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 179 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 180 END DO 181 END DO 182 END DO 183 184 185 ! III. Diagnostics and control print 186 ! ---------------------------------- 187 IF( l_trdtra ) THEN ! trend diagnostics 188 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 189 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 190 CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwz, pwn, ptn ) 191 ENDIF 192 ! ! "Poleward" heat and salt transport 193 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 194 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 195 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 196 ENDIF 197 ! ! control print 198 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - adv: ', mask1=tmask, clinfo3=cdtype ) 280 199 ! 281 200 END SUBROUTINE tra_adv_cen2 -
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_eiv.F90
r719 r791 5 5 !!====================================================================== 6 6 !! History : 9.0 ! 05-11 (G. Madec) Original code, from traldf and zdf _iso 7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_traldf_eiv || defined key_esopa … … 32 33 # include "vectopt_loop_substitute.h90" 33 34 !!---------------------------------------------------------------------- 34 !! OPA 9.0 , LOCEAN-IPSL (2006)35 !! $ Header$35 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 36 !! $Id:$ 36 37 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 37 38 !!---------------------------------------------------------------------- … … 43 44 !! *** ROUTINE tra_adv_eiv *** 44 45 !! 45 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive46 !! trend and add it to the general trend of tracer equation.47 !! 48 !! ** Method : The eddy induced advectionis computed from the slope46 !! ** Purpose : Compute the eddy induced transport and add it to the 47 !! effective transport 48 !! 49 !! ** Method : The eddy induced transport is computed from the slope 49 50 !! of iso-neutral surfaces computed in routine ldf_slp as follows: 50 !! zu_eiv = 1/(e2u e3u) dk[ aeiu e2u mi(wslpi) ] 51 !! zv_eiv = 1/(e1v e3v) dk[ aeiv e1v mj(wslpj) 52 !! zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ] 53 !! + dj[ aeiv e1v mj(wslpj) ] } 51 !! zu_eiv = dk[ aeiu e2u mi(wslpi) ] 52 !! zv_eiv = dk[ aeiv e1v mj(wslpj) ] 53 !! zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] } 54 54 !! add the eiv component to the model velocity: 55 55 !! p.n = p.n + z._eiv 56 !! CAUTION : the horizontal transports not updated along jpi column and jpj row 57 !! the vertical transports not updated along 1 & jpi columns and 1 & jpj rows 56 58 !! 57 59 !! ** Action : - add to p.n the eiv component 58 60 !!---------------------------------------------------------------------- 59 61 INTEGER , INTENT(in ) :: kt ! ocean time-step index 60 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pun ! in : 3 ocean velocitycomponents61 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pvn ! out: 3 ocean velocitycomponents62 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pun ! in : 3 ocean transport components 63 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pvn ! out: 3 ocean transport components 62 64 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pwn ! increased by the eiv 63 65 !! … … 88 90 zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1) 89 91 90 zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) / fse3u(ji,jj,jk)91 zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) / fse3v(ji,jj,jk)92 zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 93 zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 92 94 93 pun(ji,jj,jk) = pun(ji,jj,jk) + zu_eiv94 pvn(ji,jj,jk) = pvn(ji,jj,jk) + zv_eiv95 # if defined key_diaeiv 96 u_eiv(ji,jj,jk) = zu_eiv 97 v_eiv(ji,jj,jk) = zv_eiv 95 pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv 96 pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv 97 # if defined key_diaeiv 98 u_eiv(ji,jj,jk) = zu_eiv / fse3u(ji,jj,jk) 99 v_eiv(ji,jj,jk) = zv_eiv / fse3v(ji,jj,jk) 98 100 # endif 99 101 END DO … … 115 117 zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji ,jj) * vmask(ji ,jj,jk) 116 118 117 zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) & 118 & / ( e1t(ji,jj)*e2t(ji,jj) ) 119 zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 119 120 # endif 120 121 pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv 121 122 122 123 # if defined key_diaeiv 123 w_eiv(ji,jj,jk) = zw_eiv 124 w_eiv(ji,jj,jk) = zw_eiv / ( e1t(ji,jj)*e2t(ji,jj) ) 124 125 # endif 125 126 END DO … … 130 131 ! ! ================= 131 132 END SUBROUTINE tra_adv_eiv 133 134 135 !!gm test tra_adv_eiv better (faster) coded? to be verified 136 137 SUBROUTINE tra_adv_eiv2( kt, pun, pvn, pwn ) 138 !!---------------------------------------------------------------------- 139 !! *** ROUTINE tra_adv_eiv *** 140 !! 141 !! ** Purpose : Compute the eddy induced transport and add it to the 142 !! effective transport 143 !! 144 !! ** Method : The eddy induced transport is computed from the slope 145 !! of iso-neutral surfaces (see ldfslp.F90) as follows: 146 !! zu_eiv = dk[ aeiu e2u mi(wslpi) ] 147 !! zv_eiv = dk[ aeiv e1v mj(wslpj) ] 148 !! zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] } 149 !! add the eiv component to the model velocity: 150 !! p.n = p.n + z._eiv 151 !! CAUTION : the horizontal transports not updated along jpi column and jpj row 152 !! the vertical transports not updated along 1 & jpi columns and 1 & jpj rows 153 !! 154 !! ** Action : - add to p.n the eiv transport component 155 !!---------------------------------------------------------------------- 156 INTEGER , INTENT(in ) :: kt ! ocean time-step index 157 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pun ! in : 3 ocean transport components 158 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pvn ! out: 3 ocean transport components 159 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pwn ! increased by the eiv 160 !! 161 INTEGER :: ji, jj, jk ! dummy loop indices 162 REAL(wp) :: zuwslpi, zvwslpj ! temporary scalar 163 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsfu, zsfv ! eiv stream-function in u and v directions 164 !!---------------------------------------------------------------------- 165 166 IF( kt == nit000 ) THEN 167 IF(lwp) WRITE(numout,*) 168 IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :' 169 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ add to velocity fields the eiv component' 170 ENDIF 171 172 ! eiv stream-function in u- and v-directions 173 ! NB: UW-point mask at level k is umask(:,:,k) idem form VW-point mask 174 zsfu(:,:, 1 ) = 0.e0 ; zsfv(:,:, 1 ) = 0.e0 ! surface value set to zero 175 zsfu(:,:,jpk) = 0.e0 ; zsfv(:,:,jpk) = 0.e0 ! bottom value set to zero 176 DO jk = 2, jpkm1 177 DO jj = 1, jpjm1 178 DO ji = 1, fs_jpim1 ! vector opt. 179 zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 180 zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 181 zsfu(ji,jj,jk) = zuwslpi * e2u(ji,jj) * fsaeiu(ji,jj,jk) * umask(ji,jj,jk) 182 zsfv(ji,jj,jk) = zvwslpj * e1v(ji,jj) * fsaeiv(ji,jj,jk) * vmask(ji,jj,jk) 183 END DO 184 END DO 185 END DO 186 # if defined key_diaeiv 187 ! save eiv stream function in the output 188 !!gm to be done, u_sfeiv and v_sfeiv not defined ==> new IOM.... 189 !!gm and zsfu, zfv notdefined for jpi column and jpj row 190 ! u_sfeiv(:,:,:) = zsfu(:,:,:) 191 ! v_sfeiv(:,:,:) = zsfu(:,:,:) 192 # endif 193 194 ! increase the transport with the eiv transport 195 DO jk = 1, jpkm1 196 DO jj = 1, jpjm1 197 DO ji = 1, fs_jpim1 ! vector opt. 198 pun(ji,jj,jk) = pun(ji,jj,jk) + ( zsfu(ji,jj,jk) - zsfu(ji,jj,jk+1) ) 199 pvn(ji,jj,jk) = pvn(ji,jj,jk) + ( zsfv(ji,jj,jk) - zsfv(ji,jj,jk+1) ) 200 END DO 201 END DO 202 END DO 203 DO jk = 2, jpkm1 204 DO jj = 2, jpjm1 205 DO ji = fs_2, fs_jpim1 ! vector opt. 206 pwn(ji,jj,jk) = pwn(ji,jj,jk) - ( zsfu(ji,jj,jk) - zsfu(ji-1,jj ,jk) & 207 & + zsfv(ji,jj,jk) - zsfv(ji ,jj-1,jk) ) * tmask(ji,jj,jk) 208 END DO 209 END DO 210 END DO 211 ! 212 END SUBROUTINE tra_adv_eiv2 132 213 133 214 #else -
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!====================================================================== 6 !! History : !06-00 (A.Estublier) for passive tracers7 !! !01-08 (E.Durand, G.Madec) adapted for T & S8 !! NEMO 1.0 ! 02-06 (G. Madec) F90: Free form and module9 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : OPA ! 2006-00 (A.Estublier) for passive tracers 7 !! 8.2 ! 2001-08 (E.Durand, G.Madec) adapted for T & S 8 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 9 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 10 10 !!---------------------------------------------------------------------- 11 11 … … 28 28 PRIVATE 29 29 30 PUBLIC tra_adv_muscl ! routine called by step.F9030 PUBLIC tra_adv_muscl ! routine called by step.F90 31 31 32 32 !! * Substitutions … … 35 35 !!---------------------------------------------------------------------- 36 36 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 37 !! $Id :$37 !! $Id$ 38 38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 39 39 !!---------------------------------------------------------------------- … … 66 66 !! 67 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 68 REAL(wp) :: zu, z0u, zzwx 69 REAL(wp) :: zv, z0v, zzwy 70 REAL(wp) :: zw, z0w 71 REAL(wp) :: zbtr, z2, zdt, zalpha 73 72 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwx, zwy, zslpx, zslpy ! 3D workspace 74 73 !!---------------------------------------------------------------------- … … 80 79 ENDIF 81 80 82 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 83 ELSE ; z2 = 2. 81 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ! euler time-stepping 82 ELSE ; z2 = 2. ! leap-frog time-stepping 84 83 ENDIF 85 84 86 85 ! I. Horizontal advective fluxes 87 86 ! ------------------------------ 88 ! first guess of the slopes89 ! interiorvalues90 DO jk = 1, jpkm1 87 ! !-- first guess of the slopes 88 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values 89 DO jk = 1, jpkm1 ! interior values 91 90 DO jj = 1, jpjm1 92 91 DO ji = 1, fs_jpim1 ! vector opt. … … 96 95 END DO 97 96 END DO 98 ! bottom values 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. ) 103 104 ! Slopes 105 ! interior values 106 DO jk = 1, jpkm1 107 DO jj = 2, jpj 108 DO ji = fs_2, jpi ! vector opt. 97 !!gm optimisation (lbc to be suppressed) 98 CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign) 99 CALL lbc_lnk( zwy, 'V', -1. ) 100 !!gm 101 102 ! !-- Slopes of tracer 103 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values 104 DO jk = 1, jpkm1 ! interior values 105 DO jj = 2, jpjm1 106 DO ji = fs_2, jpim1 ! vector opt. 109 107 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) ) )108 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) ) 111 109 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) ) ) 113 END DO 114 END DO 115 END DO 116 ! bottom values 117 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 118 119 ! Slopes limitation 120 DO jk = 1, jpkm1 121 DO jj = 2, jpj 122 DO ji = fs_2, jpi ! vector opt. 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) ) ) 110 & * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) ) 111 END DO 112 END DO 113 END DO 114 !!gm :merge the 2 loop (above & below) 115 DO jk = 1, jpkm1 ! Slopes limitation 116 DO jj = 2, jpjm1 117 DO ji = fs_2, jpim1 ! vector opt. 118 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), & 119 & 2.*ABS( zwx (ji-1,jj,jk) ), & 120 & 2.*ABS( zwx (ji ,jj,jk) ) ) 121 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), & 122 & 2.*ABS( zwy (ji,jj-1,jk) ), & 123 & 2.*ABS( zwy (ji,jj ,jk) ) ) 131 124 END DO 132 125 END DO 133 126 END DO 134 135 ! Advection terms 136 ! interior values 137 DO jk = 1, jpkm1 138 zstep = z2 * rdttra(jk) 139 DO jj = 2, jpjm1 140 DO ji = fs_2, fs_jpim1 ! vector opt. 141 ! volume fluxes 142 #if defined key_zco 143 zeu = e2u(ji,jj) * pun(ji,jj,jk) 144 zev = e1v(ji,jj) * pvn(ji,jj,jk) 145 #else 146 zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 147 zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 148 #endif 149 ! MUSCL fluxes 127 !!gm optimisation: call lbclnk just 1 time here on zslpx & zslpy 128 !! CALL lbc_lnk( zslpx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign) 129 !! CALL lbc_lnk( zslpy, 'V', -1. ) 130 !! 131 !! and loop below from 1 to jpjm1 and to jpim1 132 !!gm 133 134 ! !-- MUSCL horizontal advective fluxes 135 DO jk = 1, jpkm1 ! interior values 136 zdt = z2 * rdttra(jk) 137 DO jj = 2, jpjm1 138 DO ji = fs_2, fs_jpim1 ! vector opt. 150 139 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 151 140 zalpha = 0.5 - z0u 152 zu = z0u - 0.5 * pun(ji,jj,jk) * z step / e1u(ji,jj)141 zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 153 142 zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 154 143 zzwy = ptb(ji ,jj,jk) + zu * zslpx(ji ,jj,jk) 155 zwx(ji,jj,jk) = zeu* ( zalpha * zzwx + (1.-zalpha) * zzwy )144 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 156 145 ! 157 146 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 158 147 zalpha = 0.5 - z0v 159 zv = z0v - 0.5 * pvn(ji,jj,jk) * z step / e2v(ji,jj)148 zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 160 149 zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 161 150 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) 151 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 152 END DO 153 END DO 154 END DO 155 !!gm optimisation (lbc to be suppressed) 156 ! ! lateral boundary conditions on zwx, zwy (changed sign) 169 157 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) 170 158 !!gm 159 160 ! !-- MUSCL advective trend 171 161 ! Tracer flux divergence at t-point added to the general trend 172 162 DO jk = 1, jpkm1 173 163 DO jj = 2, jpjm1 174 164 DO ji = fs_2, fs_jpim1 ! vector opt. 175 #if defined key_zco176 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )177 #else178 165 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 179 #endif 180 ! horizontal advective trends 181 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 182 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 183 ! add it to the general tracer trends 184 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 166 ! horizontal advective trends added to the general tracer trends 167 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 168 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 185 169 END DO 186 170 END DO 187 171 END DO 188 172 189 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 190 191 ! Save the horizontal advective trends for diagnostics 192 IF( l_trdtra ) THEN 173 IF( l_trdtra ) THEN !-- trend diagnostics 193 174 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb ) 194 175 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb ) 195 176 ENDIF 196 177 ! !-- "Poleward" heat or salt transports 197 178 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 198 IF( lk_zco ) THEN199 DO jk = 1, jpkm1200 DO jj = 2, jpjm1201 DO ji = fs_2, fs_jpim1 ! vector opt.202 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)203 END DO204 END DO205 END DO206 ENDIF207 179 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 208 180 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 209 181 ENDIF 182 ! !-- control print 183 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 210 184 211 185 … … 213 187 ! ----------------------------- 214 188 215 ! First guess of the slope216 ! interior values217 DO jk = 2, jpkm1 189 ! !-- first guess of the slopes 190 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions 191 DO jk = 2, jpkm1 ! interior values 218 192 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) ) 219 193 END DO 220 ! surface & bottom boundary conditions 221 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 222 223 ! Slopes 224 DO jk = 2, jpkm1 194 195 ! !-- Slopes of tracer 196 zslpx(:,:,1) = 0.e0 ! surface values 197 DO jk = 2, jpkm1 ! interior value 225 198 DO jj = 1, jpj 226 199 DO ji = 1, jpi … … 230 203 END DO 231 204 END DO 232 233 ! Slopes limitation234 DO jk = 2, jpkm1 ! interior values205 !!gm : merge the above and below loops 206 ! !-- Slopes limitation 207 DO jk = 2, jpkm1 ! interior values 235 208 DO jj = 1, jpj 236 209 DO ji = 1, jpi 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 245 246 ! vertical advective flux 247 DO jk = 1, jpkm1 ! interior values 248 zstep = z2 * rdttra(jk) 249 DO jj = 2, jpjm1 250 DO ji = fs_2, fs_jpim1 ! vector opt. 251 zew = pwn(ji,jj,jk+1) 210 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), & 211 & 2.*ABS( zwx (ji,jj,jk+1) ), & 212 & 2.*ABS( zwx (ji,jj,jk ) ) ) 213 END DO 214 END DO 215 END DO 216 217 ! !-- vertical advective flux 218 ! ! surface values (bottom already set to zero) 219 IF( lk_dynspg_rl .OR. lk_vvl) THEN ; zwx(:,:, 1 ) = 0.e0 ! rigid lid or variable volume 220 ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! linear free surface 221 ENDIF 222 DO jk = 1, jpkm1 ! interior values 223 zdt = z2 * rdttra(jk) 224 DO jj = 2, jpjm1 225 DO ji = fs_2, fs_jpim1 ! vector opt. 252 226 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 253 227 zalpha = 0.5 + z0w 254 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * z step / fse3w(ji,jj,jk+1)228 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt / (e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 255 229 zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 256 230 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 262 IF( lk_dynspg_rl .OR. lk_vvl) THEN ! rigid lid or variable volume: flux set to zero 263 zwx(:,:, 1 ) = 0.e0 ! surface 264 zwx(:,:,jpk) = 0.e0 ! bottom 265 ELSE ! free surface 266 zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! Surface 267 zwx(:,:,jpk) = 0.e0 ! bottom 268 269 ENDIF 270 271 ! Compute & add the vertical advective trend 231 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 232 END DO 233 END DO 234 END DO 235 236 ! !-- vertical advective trend 272 237 DO jk = 1, jpkm1 273 238 DO jj = 2, jpjm1 274 239 DO ji = fs_2, fs_jpim1 ! vector opt. 275 zbtr = 1. / fse3t(ji,jj,jk) 276 ! horizontal advective trends 277 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 278 ! add it to the general tracer trends 279 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 280 END DO 281 END DO 282 END DO 283 284 ! Save the vertical advective trends for diagnostic 285 ! ------------------------------------------------- 240 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 241 ! horizontal advective trends added to the general tracer trends 242 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 243 END DO 244 END DO 245 END DO 246 247 ! !-- trend diagnostic 286 248 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb ) 287 249 ! !-- control print 288 250 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype ) 289 251 ! -
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 1.0 ! 02-06 (G. Madec) from traadv_muscl7 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : 1.0 ! 2002-06 (G. Madec) from traadv_muscl 7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !!---------------------------------------------------------------------- 9 9 … … 26 26 PRIVATE 27 27 28 !! * Accessibility29 28 PUBLIC tra_adv_muscl2 ! routine called by step.F90 30 29 … … 34 33 !!---------------------------------------------------------------------- 35 34 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 36 !! $Id :$35 !! $Id$ 37 36 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 37 !!---------------------------------------------------------------------- … … 66 65 !! 67 66 INTEGER :: ji, jj, jk ! dummy loop indices 68 REAL(wp) :: zu, zv, zw , zeu, zev69 REAL(wp) :: z ew, zbtr, zstep, z267 REAL(wp) :: zu, zv, zw 68 REAL(wp) :: zbtr, zstep, z2 70 69 REAL(wp) :: z0u, z0v, z0w 71 REAL(wp) :: zzwx, zzwy, zalpha , zta70 REAL(wp) :: zzwx, zzwy, zalpha 72 71 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwx, zwy, zslpx, zslpy ! 3D workspace 73 72 !!---------------------------------------------------------------------- … … 98 97 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values 99 98 99 !!gm this lbc_lnk can be omitted 100 100 ! lateral boundary conditions on zwx, zwy (changed sign) 101 101 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) … … 117 117 DO jj = 2, jpj 118 118 DO ji = fs_2, jpi ! vector opt. 119 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 120 & * MIN( ABS( zslpx(ji ,jj,jk) ), & 121 & 2.*ABS( zwx (ji-1,jj,jk) ), & 122 & 2.*ABS( zwx (ji ,jj,jk) ) ) 123 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) & 124 & * MIN( ABS( zslpy(ji,jj ,jk) ), & 125 & 2.*ABS( zwy (ji,jj-1,jk) ), & 126 & 2.*ABS( zwy (ji,jj ,jk) ) ) 119 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), & 120 & 2.*ABS( zwx (ji-1,jj,jk) ), & 121 & 2.*ABS( zwx (ji ,jj,jk) ) ) 122 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), & 123 & 2.*ABS( zwy (ji,jj-1,jk) ), & 124 & 2.*ABS( zwy (ji,jj ,jk) ) ) 127 125 END DO 128 126 END DO … … 134 132 DO jj = 2, jpjm1 135 133 DO ji = fs_2, fs_jpim1 ! vector opt. 136 ! volume fluxes 137 #if defined key_zco 138 zeu = e2u(ji,jj) * pun(ji,jj,jk) 139 zev = e1v(ji,jj) * pvn(ji,jj,jk) 140 #else 141 zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 142 zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 143 #endif 144 ! MUSCL fluxes 145 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 134 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 146 135 zalpha = 0.5 - z0u 147 zu = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)148 zzwx = ptb(ji+1,jj,jk) + zu *zslpx(ji+1,jj,jk)149 zzwy = ptb(ji ,jj,jk) + zu *zslpx(ji ,jj,jk)150 zwx(ji,jj,jk) = zeu* ( zalpha * zzwx + (1.-zalpha) * zzwy )136 zu = z0u - 0.5 * pun(ji,jj,jk) * zstep / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 137 zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 138 zzwy = ptb(ji ,jj,jk) + zu * zslpx(ji ,jj,jk) 139 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 151 140 ! 152 141 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 153 142 zalpha = 0.5 - z0v 154 zv = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj)155 zzwx = ptb(ji,jj+1,jk) + zv *zslpy(ji,jj+1,jk)156 zzwy = ptb(ji,jj ,jk) + zv *zslpy(ji,jj ,jk)157 zwy(ji,jj,jk) = zev* ( zalpha * zzwx + (1.-zalpha) * zzwy )143 zv = z0v - 0.5 * pvn(ji,jj,jk) * zstep / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 144 zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 145 zzwy = ptb(ji,jj ,jk) + zv * zslpy(ji,jj ,jk) 146 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 158 147 END DO 159 148 END DO … … 166 155 DO jj = 2, jpjm1 167 156 DO ji = fs_2, fs_jpim1 ! vector opt. 168 #if defined key_zco169 157 IF( umask(ji,jj,jk) == 0. ) THEN 170 158 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN 171 zwx(ji+1,jj,jk) = e2u(ji+1,jj) *pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5159 zwx(ji+1,jj,jk) = pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5 172 160 ENDIF 173 161 IF( pun(ji-1,jj,jk) < 0. ) THEN 174 zwx(ji-1,jj,jk) = e2u(ji-1,jj) *pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji ,jj,jk) ) * 0.5162 zwx(ji-1,jj,jk) = pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji ,jj,jk) ) * 0.5 175 163 ENDIF 176 164 ENDIF 177 165 IF( vmask(ji,jj,jk) == 0. ) THEN 178 166 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN 179 zwy(ji,jj+1,jk) = e1v(ji,jj+1) *pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5167 zwy(ji,jj+1,jk) = pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5 180 168 ENDIF 181 169 IF( pvn(ji,jj-1,jk) < 0. ) THEN 182 zwy(ji,jj-1,jk) = e1v(ji,jj-1) *pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji ,jj,jk) ) * 0.5170 zwy(ji,jj-1,jk) = pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji ,jj,jk) ) * 0.5 183 171 ENDIF 184 172 ENDIF 185 #else186 IF( umask(ji,jj,jk) == 0. ) THEN187 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN188 zwx(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk) &189 & * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5190 ENDIF191 IF( pun(ji-1,jj,jk) < 0. ) THEN192 zwx(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk) &193 & * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji ,jj,jk) ) * 0.5194 ENDIF195 ENDIF196 IF( vmask(ji,jj,jk) == 0. ) THEN197 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN198 zwy(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk) &199 & * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5200 ENDIF201 IF( pvn(ji,jj-1,jk) < 0. ) THEN202 zwy(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk) &203 & * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji ,jj,jk) ) * 0.5204 ENDIF205 ENDIF206 #endif207 173 END DO 208 174 END DO … … 217 183 DO jj = 2, jpjm1 218 184 DO ji = fs_2, fs_jpim1 ! vector opt. 219 #if defined key_zco220 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )221 #else222 185 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 223 #endif 224 ! horizontal advective trends 225 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 226 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 227 ! add it to the general tracer trends 228 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 186 ! horizontal advective trends added to the general tracer trends 187 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 188 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 229 189 END DO 230 190 END DO … … 270 230 DO ji = 1, jpi 271 231 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) & 272 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )232 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 273 233 END DO 274 234 END DO … … 279 239 DO jj = 1, jpj 280 240 DO ji = 1, jpi 281 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 282 & * MIN( ABS( zslpx(ji,jj,jk ) ), & 283 & 2.*ABS( zwx (ji,jj,jk+1) ), & 284 & 2.*ABS( zwx (ji,jj,jk ) ) ) 241 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), & 242 & 2.*ABS( zwx (ji,jj,jk+1) ), & 243 & 2.*ABS( zwx (ji,jj,jk ) ) ) 285 244 END DO 286 245 END DO … … 293 252 DO jj = 2, jpjm1 294 253 DO ji = fs_2, fs_jpim1 ! vector opt. 295 zew = pwn(ji,jj,jk+1)296 254 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 297 255 zalpha = 0.5 + z0w 298 zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)299 zzwx = ptb(ji,jj,jk+1) + zw *zslpx(ji,jj,jk+1)300 zzwy = ptb(ji,jj,jk ) + zw *zslpx(ji,jj,jk )301 zwx(ji,jj,jk+1) = zew* ( zalpha * zzwx + (1.-zalpha)*zzwy )256 zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / ( e1t(ji,jj) *e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 257 zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 258 zzwy = ptb(ji,jj,jk ) + zw * zslpx(ji,jj,jk ) 259 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha)*zzwy ) 302 260 END DO 303 261 END DO … … 315 273 END DO 316 274 317 ! surface values 318 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or variable volume: flux set to zero 319 zwx(:,:, 1 ) = 0.e0 ! surface 320 zwx(:,:,jpk) = 0.e0 ! bottom 321 ELSE ! free surface 322 zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! surface 323 zwx(:,:,jpk) = 0.e0 ! bottom 324 ENDIF 275 zwx(:,:,jpk) = 0.e0 ! bottom value 276 ! ! surface values 277 IF( lk_dynspg_rl .OR. lk_vvl) THEN ; zwx(:,:, 1 ) = 0.e0 ! rigid lid or variable volume 278 ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! linear free surface 279 ENDIF 325 280 326 281 … … 330 285 DO ji = fs_2, fs_jpim1 ! vector opt. 331 286 zbtr = 1. / fse3t(ji,jj,jk) 332 ! horizontal advective trends 333 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 334 ! add it to the general tracer trends 335 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 287 ! horizontal advective trends added to the general tracer trends 288 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 336 289 END DO 337 290 END DO -
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r786 r791 5 5 !!====================================================================== 6 6 !! History : 1.0 ! 06-09 (G. Reffray) Original code 7 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !!---------------------------------------------------------------------- 9 9 … … 30 30 PUBLIC tra_adv_qck ! routine called by traadv.F90 31 31 32 REAL(wp), DIMENSION(jpi,jpj) :: zbtr233 32 REAL(wp), DIMENSION(jpi,jpj,jpk) :: sl 34 REAL(wp) :: cst1, cst2, dt, coef1 ! temporary scalars35 INTEGER :: ji, jj, jk ! dummy loop indices36 33 37 34 !! * Substitutions … … 40 37 !!---------------------------------------------------------------------- 41 38 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 42 !! $Id :$39 !! $Id$ 43 40 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 41 !!---------------------------------------------------------------------- … … 94 91 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 95 92 !! 96 REAL(wp) :: z2 ! temporary scalar 93 INTEGER :: ji, jj, jk ! dummy loop indices 94 REAL(wp) :: zdt, z2, zcoef1 ! temporary scalars 95 REAL(wp) :: zbtr ! temporary scalars 97 96 !!---------------------------------------------------------------------- 98 97 … … 101 100 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme' 102 101 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 103 104 zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )105 cst1 = 1./12.106 cst2 = 2./3.107 102 ENDIF 108 103 109 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 110 ELSE ; z2 = 2. 104 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ! euler time-stepping 105 ELSE ; z2 = 2. ! leap-frog time-stepping 111 106 ENDIF 112 107 … … 116 111 !--------------------------------------------------------------- 117 112 113 !!gm : optimisation: sl compute twice (for t & s, and even more for trc) 114 !!gm note that sl is in permanent memory be used as workspace in the vertical part ! 118 115 sl(:,:,:) = 100. 119 116 120 ! =============== 121 DO jk = 1, jpkm1 ! Horizontal slab 122 ! ! =============== 123 dt = z2 * rdttra(jk) 124 DO jj = 2, jpjm1 125 DO ji = 2, fs_jpim1 ! vector opt. 126 coef1 = 1.e-12 127 IF (pun(ji-1,jj ,jk ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj ,jk ))*dt/e1t(ji,jj) 128 IF (pun(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pun(ji ,jj ,jk ))*dt/e1t(ji,jj) 129 IF (pvn(ji ,jj-1,jk ).LT.0.) coef1 = coef1 + ABS(pvn(ji ,jj-1,jk ))*dt/e2t(ji,jj) 130 IF (pvn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pvn(ji ,jj ,jk ))*dt/e2t(ji,jj) 131 IF (pwn(ji ,jj ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk+1))*dt/(fse3t(ji,jj,jk)) 132 IF (pwn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk ))*dt/(fse3t(ji,jj,jk)) 133 sl(ji,jj,jk) = 1./coef1 134 sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk)) 135 sl(ji,jj,jk) = MAX(1. ,sl(ji,jj,jk)) 136 ENDDO 137 ENDDO 138 ENDDO 139 !--- Lateral boundary conditions on the limiter slope 140 CALL lbc_lnk( sl(:,:,:), 'T', 1. ) 117 DO jk = 1, jpkm1 118 zdt = z2 * rdttra(jk) 119 DO jj = 2, jpjm1 120 DO ji = 2, fs_jpim1 ! vector opt. 121 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) *fse3t(ji,jj,jk) ) 122 zcoef1 = 1.e-12 123 IF( pun(ji-1,jj ,jk ) < 0.e0 ) zcoef1 = zcoef1 + ABS( pun(ji-1,jj ,jk ) ) * zdt * zbtr 124 IF( pun(ji ,jj ,jk ) > 0.e0 ) zcoef1 = zcoef1 + ABS( pun(ji ,jj ,jk ) ) * zdt * zbtr 125 IF( pvn(ji ,jj-1,jk ) < 0.e0 ) zcoef1 = zcoef1 + ABS( pvn(ji ,jj-1,jk ) ) * zdt * zbtr 126 IF( pvn(ji ,jj ,jk ) > 0.e0 ) zcoef1 = zcoef1 + ABS( pvn(ji ,jj ,jk ) ) * zdt * zbtr 127 IF( pwn(ji ,jj ,jk+1) < 0.e0 ) zcoef1 = zcoef1 + ABS( pwn(ji ,jj ,jk+1) ) * zdt * zbtr 128 IF( pwn(ji ,jj ,jk ) > 0.e0 ) zcoef1 = zcoef1 + ABS( pwn(ji ,jj ,jk ) ) * zdt * zbtr 129 sl(ji,jj,jk) = 1. / zcoef1 130 sl(ji,jj,jk) = MIN( 100., sl(ji,jj,jk) ) 131 sl(ji,jj,jk) = MAX( 1. , sl(ji,jj,jk) ) 132 END DO 133 END DO 134 END DO 135 CALL lbc_lnk( sl(:,:,:), 'T', 1. ) ! Lateral boundary conditions on the limiter slope 141 136 142 137 … … 146 141 CALL tra_adv_qck_hor( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) 147 142 143 ! ! control print 148 144 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - had: ', mask1=tmask, clinfo3=cdtype ) 149 145 … … 152 148 !------------------------------------------------------------------------- 153 149 154 CALL tra_adv_qck_ver( pwn, ptn , pta, z2 ) 155 150 CALL tra_adv_qck_ver( pwn, ptn , pta ) 151 152 ! ! control print 156 153 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - zad: ', mask1=tmask, clinfo3=cdtype ) 157 154 ! … … 159 156 160 157 161 SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) 162 !!---------------------------------------------------------------------- 158 SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, pz2 ) 159 !!---------------------------------------------------------------------- 160 !! *** ROUTINE tra_adv_qck_hor *** 161 !! 162 !! ** Purpose : 163 163 !! 164 164 !!---------------------------------------------------------------------- … … 166 166 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 167 167 INTEGER , INTENT(in ) :: ktra ! tracer index 168 REAL(wp) , INTENT(in ) :: z2 ! ???168 REAL(wp) , INTENT(in ) :: pz2 ! =1 or 2 (euler or leap-frog) 169 169 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! horizontal effective velocity 170 170 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptb ! before tracer field 171 171 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 172 173 REAL(wp) :: za, zbtr, e1, e2, c, dir, fu, fc, fd ! temporary scalars 174 REAL(wp) :: coef2, coef3, fho, mask, dx 175 REAL(wp), DIMENSION(jpi,jpj) :: zee 176 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmask, zlap, dwst, lim 172 !! 173 INTEGER :: ji, jj, jk ! dummy loop indices 174 REAL(wp) :: zdt ! temporary scalars 175 REAL(wp) :: zbtr, zc, zdir, zfu, zfc, zfd ! temporary scalars 176 REAL(wp) :: zcoef1, zcoef2, zcoef3, zfho, zmsk, zdx 177 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmask, zlap, zdwst, zlim 177 178 !!---------------------------------------------------------------------- 178 179 … … 182 183 zmask(:,:,:)= 0.e0 183 184 zlap (:,:,:)= 0.e0 184 dwst(:,:,:)= 0.e0185 lim(:,:,:)= 0.e0185 zdwst(:,:,:)= 0.e0 186 zlim (:,:,:)= 0.e0 186 187 187 188 !---------------------------------------------------------------------- … … 189 190 !---------------------------------------------------------------------- 190 191 191 ! =============== 192 DO jk = 1, jpkm1 ! Horizontal slab 193 ! ! =============== 194 ! Initialization of metric arrays (for z- or s-coordinates) 195 ! --------------------------------------------------------- 196 DO jj = 1, jpjm1 197 DO ji = 1, fs_jpim1 ! vector opt. 198 #if defined key_zco 199 ! z-coordinates, no vertical scale factors 200 zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) 201 #else 202 ! vertical scale factor are used 203 zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 204 #endif 205 END DO 206 END DO 192 DO jk = 1, jpkm1 207 193 208 194 ! Laplacian of tracers (at before time step) … … 211 197 DO jj = 1, jpjm1 212 198 DO ji = 1, fs_jpim1 ! vector opt. 213 zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji+1,jj ,jk) - ptb(ji,jj,jk) ) 214 END DO 215 END DO 216 DO jj = 2, jpjm1 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 #if defined key_zco 219 zee(ji,jj) = e1t(ji,jj) / e2t(ji,jj) 220 #else 221 zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk)) 222 #endif 223 zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji-1,jj,jk) ) 199 zmask(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) & 200 & * ( ptb(ji+1,jj ,jk) - ptb(ji,jj,jk) ) / e1u(ji,jj) 201 END DO 202 END DO 203 DO jj = 2, jpjm1 204 DO ji = fs_2, fs_jpim1 ! vector opt. 205 zlap(ji,jj,jk) = e1t(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji-1,jj,jk) ) & 206 & / ( e2t(ji,jj) * fse3t(ji,jj,jk) ) 224 207 END DO 225 208 END DO … … 231 214 DO ji = 2, fs_jpim1 ! vector opt. 232 215 ! Upstream in the x-direction for the tracer 233 zmask(ji,jj,jk) = ptb(ji-1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji-1,jj,jk))216 zmask(ji,jj,jk) = ptb(ji-1,jj,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji-1,jj,jk) ) 234 217 ! Downstream in the x-direction for the tracer 235 dwst (ji,jj,jk) =ptb(ji+1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji+1,jj,jk)) 236 ENDDO 237 ENDDO 238 END DO 239 !--- Lateral boundary conditions on the laplacian (unchanged sgn) 240 CALL lbc_lnk( zlap(:,:,:), 'T', 1. ) 241 !--- Lateral boundary conditions for the lim function 242 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. ) 243 ! =============== 244 DO jk = 1, jpkm1 ! Horizontal slab 245 ! ! =============== 246 ! --- lim at the U-points in function of the direction of the flow 247 ! ---------------------------------------------------------------- 218 zdwst(ji,jj,jk) = ptb(ji+1,jj,jk) + sl(ji,jj,jk) * ( ptb(ji,jj,jk) - ptb(ji+1,jj,jk) ) 219 END DO 220 END DO 221 END DO 222 223 !--- Lateral boundary conditions on the laplacian and lim functions (unchanged sgn) 224 CALL lbc_lnk( zlap (:,:,:), 'T', 1. ) 225 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zdwst(:,:,:), 'T', 1. ) 226 227 ! --- lim at the U-points in function of the direction of the flow 228 ! ---------------------------------------------------------------- 229 DO jk = 1, jpkm1 248 230 DO jj = 1, jpjm1 249 231 DO ji = 1, fs_jpim1 ! vector opt. 250 dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0251 lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk)232 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun>0 : diru = 1 otherwise diru = 0 233 zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji+1,jj,jk) 252 234 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 253 zmask(ji,jj,jk) =tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2254 END DO 255 END DO 256 END DO 257 !--- Lateral boundary conditions for the mask258 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) 235 zmask(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2 236 END DO 237 END DO 238 END DO 239 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) !--- Lateral boundary conditions for the mask 240 259 241 260 242 ! Horizontal advective fluxes 261 243 ! --------------------------- 262 ! =============== 263 DO jk = 1, jpkm1 ! Horizontal slab 264 ! =============== 265 dt = z2 * rdttra(jk) 244 DO jk = 1, jpkm1 245 zdt = pz2 * rdttra(jk) 266 246 !--- tracer flux at u and v-points 267 247 DO jj = 1, jpjm1 268 248 DO ji = 1, fs_jpim1 ! vector opt. 269 #if defined key_zco 270 e2 = e2u(ji,jj) 271 #else 272 e2 = e2u(ji,jj) * fse3u(ji,jj,jk) 273 #endif 274 dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0 275 276 dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj) 277 c = ABS(pun(ji,jj,jk))*dt/dx ! (0<cx<1 : Courant number on x-direction) 278 279 fu = lim(ji,jj,jk) ! FU + sl(FC-FU) in the x-direction for T 280 fc = dir*ptb(ji ,jj,jk)+(1-dir)*ptb(ji+1,jj,jk) ! FC in the x-direction for T 281 fd = dir*ptb(ji+1,jj,jk)+(1-dir)*ptb(ji ,jj,jk) ! FD in the x-direction for T 249 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun>0 : diru = 1 otherwise diru = 0 250 zdx = ( zdir * e1t(ji,jj) + (1-zdir)* e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk) 251 zc = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<cx<1 : Courant number on x-direction) 252 253 zfu = zlim(ji,jj,jk) ! FU + sl(FC-FU) in the x-direction for T 254 zfc = zdir * ptb(ji ,jj,jk) + (1-zdir) * ptb(ji+1,jj,jk) ! FC in the x-direction for T 255 zfd = zdir * ptb(ji+1,jj,jk) + (1-zdir) * ptb(ji ,jj,jk) ! FD in the x-direction for T 282 256 283 257 !--- QUICKEST scheme 284 258 ! Temperature on the x-direction 285 coef1 = 0.5*(fc+fd)286 coef2 = 0.5*c*(fd-fc)287 coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji+1,jj,jk) )288 fho = coef1-coef2-coef3289 fho = bound(fu,fd,fc,fho)259 zcoef1 = 0.5 * ( zfc + zfd ) 260 zcoef2 = 0.5 * zc * ( zfd - zfc ) 261 zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji+1,jj,jk) ) 262 zfho = zcoef1 - zcoef2 - zcoef3 263 zfho = bound( zfu, zfd, zfc, zfho ) 290 264 !--- If the second ustream point is a land point 291 265 !--- the flux is computed by the 1st order UPWIND scheme 292 mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji+1,jj,jk)293 fho = mask*fho + (1-mask)*fc294 dwst(ji,jj,jk)=e2*pun(ji,jj,jk)*fho266 zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji+1,jj,jk) 267 zfho = zmsk * zfho + (1-zmsk) * zfc 268 zdwst(ji,jj,jk) = pun(ji,jj,jk) * zfho 295 269 END DO 296 270 END DO … … 300 274 DO ji = fs_2, fs_jpim1 ! vector opt. 301 275 !--- horizontal advective trends 302 #if defined key_zco 303 zbtr = zbtr2(ji,jj) 304 #else 305 zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) 306 #endif 307 za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj ,jk) ) 276 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 308 277 !--- add it to the general tracer trends 309 pta(ji,jj,jk) = pta(ji,jj,jk) + za 310 END DO 311 END DO 312 ! ! =============== 313 END DO ! End of slab 314 ! ! =============== 315 316 ! Save the horizontal advective trends for diagnostic 317 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, dwst, pun, ptb ) 278 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zdwst(ji,jj,jk) - zdwst(ji-1,jj ,jk) ) 279 END DO 280 END DO 281 END DO 282 283 ! !-- trend diagnostic 284 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zdwst, pun, ptb ) 318 285 319 286 … … 321 288 ! I. Part 2 : y-direction 322 289 !---------------------------------------------------------------------- 323 ! ============== 324 DO jk = 1, jpkm1 ! Horizontal slab 325 ! ! =============== 326 ! Initialization of metric arrays (for z- or s-coordinates) 327 ! --------------------------------------------------------- 328 DO jj = 1, jpjm1 329 DO ji = 1, fs_jpim1 ! vector opt. 330 #if defined key_zco 331 ! z-coordinates, no vertical scale factors 332 zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) 333 #else 334 ! s-coordinates, vertical scale factor are used 335 zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 336 #endif 337 END DO 338 END DO 290 291 DO jk = 1, jpkm1 339 292 340 293 ! Laplacian of tracers (at before time step) … … 343 296 DO jj = 1, jpjm1 344 297 DO ji = 1, fs_jpim1 ! vector opt. 345 zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji ,jj+1,jk) - ptb(ji,jj,jk) ) 298 zmask(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) & 299 & * ( ptb(ji ,jj+1,jk) - ptb(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk) 346 300 END DO 347 301 END DO … … 349 303 DO jj = 2, jpjm1 350 304 DO ji = fs_2, fs_jpim1 ! vector opt. 351 #if defined key_zco 352 zee(ji,jj) = e2t(ji,jj) / e1t(ji,jj) 353 #else 354 zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk)) 355 #endif 356 zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji,jj-1,jk) ) 305 zlap(ji,jj,jk) = e2t(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji,jj-1,jk) ) & 306 & / ( e1t(ji,jj) * fse3t(ji,jj,jk) ) 357 307 END DO 358 308 END DO … … 362 312 DO ji = 2, fs_jpim1 ! vector opt. 363 313 ! Upstream in the y-direction for the tracer 364 zmask(ji,jj,jk) =ptb(ji,jj-1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj-1,jk))314 zmask(ji,jj,jk) = ptb(ji,jj-1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj-1,jk) ) 365 315 ! Downstream in the y-direction for the tracer 366 dwst (ji,jj,jk)=ptb(ji,jj+1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj+1,jk))316 zdwst(ji,jj,jk) = ptb(ji,jj+1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj+1,jk) ) 367 317 ENDDO 368 318 ENDDO 369 319 END DO 370 !--- Lateral boundary conditions on the laplacian (unchanged sgn) 371 CALL lbc_lnk( zlap(:,:,:), 'T', 1. ) 372 !--- Lateral boundary conditions for the lim function 373 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. ) 374 375 DO jk = 1, jpkm1 ! Horizontal slab 376 ! ! =============== 320 !--- Lateral boundary conditions on the laplacian and lim function (unchanged sgn) 321 CALL lbc_lnk( zlap (:,:,:), 'T', 1. ) 322 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zdwst(:,:,:), 'T', 1. ) 323 324 DO jk = 1, jpkm1 325 377 326 ! --- lim at the V-points in function of the direction of the flow 378 327 ! ---------------------------------------------------------------- 379 328 DO jj = 1, jpjm1 380 329 DO ji = 1, fs_jpim1 ! vector opt. 381 dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0382 lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk)330 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pvn>0 : dirv = 1 otherwise dirv = 0 331 zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji,jj+1,jk) 383 332 ! Mask at the T-points in the y-direction (mask=0 or mask=1) 384 zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 385 END DO 386 END DO 387 END DO 388 !--- Lateral boundary conditions for the mask 389 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) 333 zmask(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 334 END DO 335 END DO 336 END DO 337 CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) !--- Lateral boundary conditions for the mask 390 338 391 339 ! Horizontal advective fluxes 392 ! ------------------------------- 393 ! =============== 394 DO jk = 1, jpkm1 ! Horizontal slab 395 ! =============== 396 dt = z2 * rdttra(jk) 340 ! --------------------------- 341 342 DO jk = 1, jpkm1 343 zdt = pz2 * rdttra(jk) 397 344 !--- tracer flux at u and v-points 398 345 DO jj = 1, jpjm1 399 346 DO ji = 1, fs_jpim1 ! vector opt. 400 #if defined key_zco 401 e1 = e1v(ji,jj) 402 #else 403 e1 = e1v(ji,jj) * fse3v(ji,jj,jk) 404 #endif 405 dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0 406 407 dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1) 408 c = ABS(pvn(ji,jj,jk))*dt/dx ! (0<cy<1 : Courant number on y-direction) 409 410 fu = lim(ji,jj,jk) ! FU + sl(FC-FU) in the y-direction for T 411 fc = dir*ptb(ji,jj ,jk)+(1-dir)*ptb(ji,jj+1,jk) ! FC in the y-direction for T 412 fd = dir*ptb(ji,jj+1,jk)+(1-dir)*ptb(ji,jj ,jk) ! FD in the y-direction for T 347 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pvn>0 : dirv = 1 otherwise dirv = 0 348 zdx = ( zdir * e2t(ji,jj) + (1-zdir)* e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk) 349 zc = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<cy<1 : Courant number on y-direction) 350 351 zfu = zlim(ji,jj,jk) ! FU + sl(FC-FU) in the y-direction for T 352 zfc = zdir * ptb(ji,jj ,jk) + (1-zdir) * ptb(ji,jj+1,jk) ! FC in the y-direction for T 353 zfd = zdir * ptb(ji,jj+1,jk) + (1-zdir) * ptb(ji,jj ,jk) ! FD in the y-direction for T 413 354 414 355 !--- QUICKEST scheme 415 356 ! Temperature on the y-direction 416 coef1 = 0.5*(fc+fd)417 coef2 = 0.5*c*(fd-fc)418 coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji,jj+1,jk) )419 fho = coef1-coef2-coef3420 fho = bound(fu,fd,fc,fho)357 zcoef1 = 0.5 * ( zfc + zfd ) 358 zcoef2 = 0.5 * zc * ( zfd - zfc ) 359 zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji,jj+1,jk) ) 360 zfho = zcoef1 - zcoef2 - zcoef3 361 zfho = bound( zfu, zfd, zfc, zfho ) 421 362 !--- If the second ustream point is a land point 422 363 !--- the flux is computed by the 1st order UPWIND scheme 423 mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji,jj+1,jk)424 fho = mask*fho + (1-mask)*fc425 dwst(ji,jj,jk)=e1*pvn(ji,jj,jk)*fho364 zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji,jj+1,jk) 365 zfho = zmsk * zfho + (1-zmsk) * zfc 366 zdwst(ji,jj,jk) = pvn(ji,jj,jk) * zfho 426 367 END DO 427 368 END DO … … 430 371 DO jj = 2, jpjm1 431 372 DO ji = fs_2, fs_jpim1 ! vector opt. 432 !--- horizontal advective trends 433 #if defined key_zco 434 zbtr = zbtr2(ji,jj) 435 #else 436 zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) 437 #endif 438 za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji ,jj-1,jk) ) 439 !--- add it to the general tracer trends 440 pta(ji,jj,jk) = pta(ji,jj,jk) + za 441 END DO 442 END DO 443 ! ! =============== 444 END DO ! End of slab 445 ! ! =============== 446 447 ! Save the horizontal advective trends for diagnostic 448 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, dwst, pvn, ptb ) 449 450 ! "zonal" mean advective heat and salt transport 373 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 374 ! horizontal advective trends added to the general tracer trends 375 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zdwst(ji,jj,jk) - zdwst(ji,jj-1,jk) ) 376 END DO 377 END DO 378 END DO 379 380 ! !-- trend diagnostic 381 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zdwst, pvn, ptb ) 382 ! ! "Poleward" heat or salt transport 451 383 IF( ln_diaptr .AND. cdtype == 'TRA' .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 452 #if defined key_zco 453 DO jk = 1, jpkm1 454 DO jj = 2, jpjm1 455 DO ji = fs_2, fs_jpim1 ! vector opt. 456 dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk) 457 END DO 458 END DO 459 END DO 460 # endif 461 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( dwst(:,:,:) ) 462 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( dwst(:,:,:) ) 384 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zdwst(:,:,:) ) 385 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zdwst(:,:,:) ) 463 386 ENDIF 464 387 ! … … 466 389 467 390 468 SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta, z2 ) 469 !!---------------------------------------------------------------------- 470 !! 471 !!---------------------------------------------------------------------- 472 REAL(wp), INTENT(in ) :: z2 473 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn 474 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn 475 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta 476 !! 477 REAL(wp) :: za, ze3tr, dt, dir, fc, fd ! temporary scalars 391 SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta ) 392 !!---------------------------------------------------------------------- 393 !! *** ROUTINE tra_adv_qck_ver *** 394 !! 395 !! ** Purpose : 396 !! 397 !!---------------------------------------------------------------------- 398 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical transport 399 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now tracer 400 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 401 !! 402 INTEGER :: ji, jj , jk ! dummy loop indices 403 REAL(wp) :: zbtr, zdir, zfc, zfd ! temporary scalars 404 !!---------------------------------------------------------------------- 478 405 479 406 ! Vertical advection … … 483 410 ! ---------------------------- 484 411 485 sl(:,:,jpk) = 0.e0 !Bottom value : flux set to zero 486 487 ! Surface value 488 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid : flux set to zero 489 sl(:,:, 1 ) = 0.e0 490 ELSE ! free surface-constant volume 491 sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) 412 sl(:,:,jpk) = 0.e0 ! bottom value 413 414 ! ! Surface value 415 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ; sl(:,:, 1 ) = 0.e0 ! rigid lid or non-linear fs 416 ELSE ; sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! linear free surface 492 417 ENDIF 418 419 !!gm bug: code au true 2nd order scheme 420 !!gm : sl used as work array: not good idea for optimisation (sl compute once for all tracers...) 493 421 494 422 ! Second order centered tracer flux at w-point 495 423 DO jk = 2, jpkm1 496 dt = z2 * rdttra(jk) 497 DO jj = 2, jpjm1 498 DO ji = fs_2, fs_jpim1 ! vector opt. 499 dir = 0.5 + sign(0.5,pwn(ji,jj,jk)) ! if pwn>0 : dirw = 1 otherwise dirw = 0 500 fc = dir*ptn(ji,jj,jk )*fse3t(ji,jj,jk-1)+(1-dir)*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk ) ! FC in the z-direction for T 501 fd = dir*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk )+(1-dir)*ptn(ji,jj,jk )*fse3t(ji,jj,jk-1) ! FD in the z-direction for T 424 DO jj = 2, jpjm1 425 DO ji = fs_2, fs_jpim1 ! vector opt. 426 zdir = 0.5 + SIGN( 0.5, pwn(ji,jj,jk) ) ! if pwn>0 : dirw = 1 otherwise dirw = 0 427 zfc = zdir * ptn(ji,jj,jk ) * fse3t(ji,jj,jk-1) + (1-zdir) * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk ) ! FC 428 zfd = zdir * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk ) + (1-zdir) * ptn(ji,jj,jk ) * fse3t(ji,jj,jk-1) ! FD 502 429 !--- Second order centered scheme 503 sl(ji,jj,jk) =pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk))430 sl(ji,jj,jk) = pwn(ji,jj,jk) * ( zfc + zfd ) / ( fse3t(ji,jj,jk-1) + fse3t(ji,jj,jk) ) 504 431 END DO 505 432 END DO … … 511 438 DO jj = 2, jpjm1 512 439 DO ji = fs_2, fs_jpim1 ! vector opt. 513 ze3tr = 1. / fse3t(ji,jj,jk) 514 ! vertical advective trends 515 za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) 516 ! add it to the general tracer trends 517 pta(ji,jj,jk) = pta(ji,jj,jk) + za 440 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 441 ! vertical advective trends added to the general tracer trends 442 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) 518 443 END DO 519 444 END DO … … 523 448 524 449 525 REAL FUNCTION bound(fu,fd,fc,fho) 526 REAL(wp) :: fu, fd, fc, fho, fref1, fref2 527 fref1 = fu 528 fref2 = MAX( MIN( fc , fd ), MIN( MAX( fc , fd ), fref1 ) ) 529 bound = MAX( MIN( fho, fc ), MIN( MAX( fho, fc ), fref2 ) ) 450 FUNCTION bound( pfu, pfd, pfc, pfho ) RESULT( pbound ) 451 !!---------------------------------------------------------------------- 452 !! *** FUNCTION bound *** 453 !! 454 !! ** Purpose : ??? 455 !!---------------------------------------------------------------------- 456 REAL(wp), INTENT(in) :: pfu, pfd, pfc, pfho 457 REAL(wp) :: pbound 458 !! 459 REAL(wp) :: zfref1, zfref2 460 !!---------------------------------------------------------------------- 461 zfref1 = pfu 462 zfref2 = MAX( MIN( pfc , pfd ), MIN( MAX( pfc , pfd ), zfref1 ) ) 463 pbound = MAX( MIN( pfho, pfc ), MIN( MAX( pfho, pfc ), zfref2 ) ) 464 ! 530 465 END FUNCTION 531 466 -
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 ! -------------------- -
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r786 r791 4 4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 1.0 ! 06-08 (L. Debreu, R. Benshila) Original code7 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : 1.0 ! 2006-08 (L. Debreu, R. Benshila) Original code 7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !!---------------------------------------------------------------------- 9 9 … … 27 27 PUBLIC tra_adv_ubs ! routine called by traadv module 28 28 29 REAL(wp), DIMENSION(jpi,jpj) :: e1e2tr ! = 1/(e1t * e2t)30 31 29 !! * Substitutions 32 30 # include "domzgr_substitute.h90" … … 34 32 !!---------------------------------------------------------------------- 35 33 !! NEMO/OPA & TRP 2.4 , LOCEAN-IPSL (2008) 36 !! $Id :$34 !! $Id$ 37 35 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 36 !!---------------------------------------------------------------------- … … 71 69 !! 72 70 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 73 !! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731 –1741.71 !! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731-1741. 74 72 !!---------------------------------------------------------------------- 75 73 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 80 78 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 81 79 !! 82 INTEGER :: ji, jj, jk ! dummy loop indices 83 REAL(wp) :: zta, zbtr, zcoef ! temporary scalars 84 REAL(wp) :: zfui, zfp_ui, zfm_ui, zcenut ! " " 85 REAL(wp) :: zfvj, zfp_vj, zfm_vj, zcenvt ! " " 86 REAL(wp) :: z_hdivn ! " " 87 REAL(wp), DIMENSION(jpi,jpj) :: zeeu, zeev ! temporary 2D workspace 88 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy ! temporary 3D workspace 89 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu , ztv , zltu , zltv, ztrdt ! " " 90 !!---------------------------------------------------------------------- 91 80 INTEGER :: ji, jj, jk ! dummy loop indices 81 REAL(wp) :: zbtr, zcoef, z_hdivn ! temporary scalars 82 REAL(wp) :: zeeu, zfp_ui, zfm_ui, zcenut ! " " 83 REAL(wp) :: zeev, zfp_vj, zfm_vj, zcenvt ! " " 84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, ztu, zltu ! temporary 3D workspace 85 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwy, ztv, zltv ! " " 86 !!---------------------------------------------------------------------- 87 88 !!gm this can be optimized: we don't need zero everywhere 92 89 zltu(:,:,:) = 0.e0 93 90 zltv(:,:,:) = 0.e0 … … 97 94 IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme' 98 95 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 99 ! 100 e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 101 ENDIF 102 103 ! store pta trends 104 ztrdt(:,:,:) = pta(:,:,:) 105 106 zcoef = 1./6. 107 ! ! =============== 108 DO jk = 1, jpkm1 ! Horizontal slab 109 ! ! =============== 110 111 ! Initialization of metric arrays (for z- or s-coordinates) 96 ENDIF 97 98 ! Laplacian 99 DO jk = 1, jpkm1 100 DO jj = 1, jpjm1 ! First derivative (gradient) 101 DO ji = 1, fs_jpim1 ! vector opt. 102 zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 103 zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 104 ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk) - ptb(ji,jj,jk) ) 105 ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk) - ptb(ji,jj,jk) ) 106 END DO 107 END DO 108 DO jj = 2, jpjm1 ! Second derivative (divergence) 109 DO ji = fs_2, fs_jpim1 ! vector opt. 110 zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 111 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 112 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef 113 END DO 114 END DO 115 END DO 116 117 ! Lateral boundary conditions on the laplacian (zlt,zls) (unchanged sgn) 118 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) 119 120 121 ! Horizontal advective fluxes 122 DO jk = 1, jpkm1 112 123 DO jj = 1, jpjm1 113 124 DO ji = 1, fs_jpim1 ! vector opt. 114 #if defined key_zco 115 ! z-coordinates, no vertical scale factors 116 zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) 117 zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) 118 #else 119 ! s-coordinates, vertical scale factor are used 120 zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 121 zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 122 #endif 123 END DO 124 END DO 125 126 ! Laplacian 127 ! First derivative (gradient) 128 DO jj = 1, jpjm1 129 DO ji = 1, fs_jpim1 ! vector opt. 130 ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj ,jk) - ptb(ji,jj,jk) ) 131 ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji ,jj+1,jk) - ptb(ji,jj,jk) ) 132 END DO 133 END DO 134 ! Second derivative (divergence) 135 DO jj = 2, jpjm1 136 DO ji = fs_2, fs_jpim1 ! vector opt. 137 #if ! defined key_zco 138 zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 139 #endif 140 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 141 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef 142 END DO 143 END DO 144 ! ! ================= 145 END DO ! End of slab 146 ! ! ================= 147 148 ! Lateral boundary conditions on the laplacian (zlt,zls) (unchanged sgn) 149 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) 150 151 ! ! =============== 152 DO jk = 1, jpkm1 ! Horizontal slab 153 ! ! =============== 154 ! Horizontal advective fluxes 155 DO jj = 1, jpjm1 156 DO ji = 1, fs_jpim1 ! vector opt. 157 ! volume fluxes * 1/2 158 #if defined key_zco 159 zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) 160 zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) 161 #else 162 zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 163 zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 164 #endif 165 ! upstream scheme 166 zfp_ui = zfui + ABS( zfui ) 167 zfp_vj = zfvj + ABS( zfvj ) 168 zfm_ui = zfui - ABS( zfui ) 169 zfm_vj = zfvj - ABS( zfvj ) 125 ! upstream transport 126 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) ! = 2.pun if pun>0, =0 if pun<0 127 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) ! = 0 if pun>0, = 2.pun if pun<0 128 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) ! idem on pvn 129 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) ! 170 130 ! centered scheme 171 zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj ,jk) ) 172 zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji ,jj+1,jk) ) 173 ! mixed centered / upstream scheme 174 zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk) 175 zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk) 176 END DO 177 END DO 178 179 ! Tracer flux divergence at t-point added to the general trend 180 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 ! horizontal advective trends 183 #if defined key_zco 184 zbtr = e1e2tr(ji,jj) 185 #else 186 zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 187 #endif 188 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 189 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 190 ! add it to the general tracer trends 191 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 192 END DO 193 END DO 194 ! ! =============== 195 END DO ! End of slab 196 ! ! =============== 197 198 ! Horizontal trend used in tra_adv_ztvd subroutine 199 zltu(:,:,:) = pta(:,:,:) - ztrdt(:,:,:) 200 201 ! Save the horizontal advective trends for diagnostic 202 ! ----------------------------------------------------- 203 IF( l_trdtra ) THEN 131 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj ,jk) ) 132 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji ,jj+1,jk) ) 133 ! UBS scheme 134 zwx(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk) ) 135 zwy(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk) ) 136 END DO 137 END DO 138 END DO 139 140 zltu(:,:,:) = pta(:,:,:) ! store pta trends 141 142 ! Horizontal advective trends 143 DO jk = 1, jpkm1 144 DO jj = 2, jpjm1 145 DO ji = fs_2, fs_jpim1 ! vector opt. 146 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 147 ! horizontal advective trends added to the general tracer trends 148 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 149 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 150 END DO 151 END DO 152 END DO 153 154 zltu(:,:,:) = pta(:,:,:) - zltu(:,:,:) ! store the Horizontal advective trend (used in tra_adv_ztvd subroutine) 155 156 157 IF( l_trdtra ) THEN ! trend diagnostics (from advective fluxes) 204 158 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 205 159 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 206 160 ENDIF 207 208 ! "Poleward" heat or salt transport 161 ! ! "Poleward" heat or salt transports 209 162 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 210 IF( lk_zco ) THEN 211 DO jk = 1, jpkm1 212 DO jj = 2, jpjm1 213 DO ji = fs_2, fs_jpim1 ! vector opt. 214 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 215 END DO 216 END DO 217 END DO 218 ENDIF 219 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 220 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 221 ENDIF 222 163 IF( ktra == jp_tem ) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 164 IF( ktra == jp_sal ) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 165 ENDIF 166 ! ! control print 223 167 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - had: ', mask1=tmask, clinfo3=cdtype ) 224 168 … … 226 170 ! II. Vertical advection 227 171 ! ---------------------- 228 IF( l_trdtra ) z trdt(:,:,:) = pta(:,:,:) ! Save ta and sa trends172 IF( l_trdtra ) zltv(:,:,:) = pta(:,:,:) ! store pta if trend diag. 229 173 230 174 ! TVD scheme the vertical direction 231 175 CALL tra_adv_ztvd( kt, pwn, zltu, ptb, ptn, pta ) 232 176 233 IF( l_trdtra ) THEN !vertical advective trend diagnostics234 DO jk = 1, jpkm1 177 IF( l_trdtra ) THEN ! vertical advective trend diagnostics 178 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 235 179 DO jj = 2, jpjm1 236 180 DO ji = fs_2, fs_jpim1 ! vector opt. 237 z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 238 ztrdt(ji,jj,jk) = pta(ji,jj,jk) - ztrdt(ji,jj,jk+1) + ptn(ji,jj,jk) * z_hdivn 181 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 182 z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr 183 zltv(ji,jj,jk) = pta(ji,jj,jk) - zltv(ji,jj,jk+1) + ptn(ji,jj,jk) * z_hdivn 239 184 END DO 240 185 END DO 241 186 END DO 242 CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=ztrdt ) 243 ! 244 ENDIF 245 187 CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=zltv ) 188 ENDIF 189 ! ! control print 246 190 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - zad: ', mask1=tmask, clinfo3=cdtype ) 247 191 ! … … 249 193 250 194 251 SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, ptb, ptn, pta )195 SUBROUTINE tra_adv_ztvd( kt, pwn, phtrd, ptb, ptn, pta ) 252 196 !!---------------------------------------------------------------------- 253 197 !! *** ROUTINE tra_adv_ztvd *** 254 198 !! 255 !! ** Purpose : Compute the now trend due to total advection of256 !! tracers and add it to the general trend of tracer equations199 !! ** Purpose : Compute the vertical advective trend of a tracer 200 !! and add it to its general trend 257 201 !! 258 202 !! ** Method : TVD scheme, i.e. 2nd order centered scheme with 259 203 !! corrected flux (monotonic correction) 260 !! note: - this advection scheme needs a leap-frog time scheme 261 !! 262 !! ** Action : - update (ta,sa) with the now advective tracer trends 263 !! - save the trends in (ztrdt,ztrds) ('key_trdtra') 264 !!---------------------------------------------------------------------- 265 INTEGER , INTENT(in ) :: kt ! ocean time-step 266 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! verical effective velocity 267 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: zttrd ! lateral advective trends on T & S 268 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: ptb, ptn ! before and now tracer fields 269 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 204 !! 205 !! ** Action : - update pta with the vertical advective trend 206 !! - trend diagnostic (lk_trdtra=T) 207 !!---------------------------------------------------------------------- 208 INTEGER , INTENT(in ) :: kt ! ocean time-step 209 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! verical effective velocity 210 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: phtrd ! lateral advective trends on T & S 211 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: ptb, ptn ! before and now tracer fields 212 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer trend 270 213 !! 271 214 INTEGER :: ji, jj, jk ! dummy loop indices 272 REAL(wp) :: z2dtt, zbtr, z ew, z2! temporary scalar273 REAL(wp) :: ztak, zfp_wk, zfm_wk 215 REAL(wp) :: z2dtt, zbtr, z2 ! temporary scalar 216 REAL(wp) :: ztak, zfp_wk, zfm_wk ! " " 274 217 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztw ! temporary 3D workspace 275 218 !!---------------------------------------------------------------------- … … 285 228 ENDIF 286 229 287 ! Bottom value : flux set to zero288 ! --------------289 ztw(:,:,jpk) = 0.e0 ; zti (:,:,:) = 0.e0290 291 230 292 231 ! upstream advection with initial mass fluxes & intermediate update 293 ! ------------------------------------------------------------------- 294 ! Surface value 295 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid : flux set to zero 296 ztw(:,:,1) = 0.e0 297 ELSE ! free surface 298 ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1) 299 ENDIF 300 301 ! Interior value 302 DO jk = 2, jpkm1 232 ztw(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 233 zti(:,:,jpk) = 0.e0 234 ! ! Surface value : 235 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ; ztw(:,:, 1 ) = 0.e0 ! rigid lid or non-linear fs 236 ELSE ; ztw(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! linear free surface 237 ENDIF 238 ! 239 DO jk = 2, jpkm1 ! interior values 303 240 DO jj = 1, jpj 304 241 DO ji = 1, jpi 305 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 306 zfp_wk = zew + ABS( zew ) 307 zfm_wk = zew - ABS( zew ) 308 ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) 242 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 243 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 244 ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) ) 309 245 END DO 310 246 END DO … … 318 254 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 319 255 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 320 pta(ji,jj,jk) = pta(ji,jj,jk) + ztak 321 zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 322 END DO 323 END DO 324 END DO 325 326 ! Lateral boundary conditions on zti, zsi (unchanged sign) 327 CALL lbc_lnk( zti, 'T', 1. ) 256 pta(ji,jj,jk) = pta(ji,jj,jk) + ztak 257 zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + phtrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 258 END DO 259 END DO 260 END DO 261 ! 262 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 328 263 329 264 330 265 ! antidiffusive flux : high order minus low order 331 ! -------------------------------------------------332 266 ztw(:,:,1) = 0.e0 ! Surface value 333 334 267 DO jk = 2, jpkm1 ! Interior value 335 268 DO jj = 1, jpj 336 269 DO ji = 1, jpi 337 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 338 ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 339 END DO 340 END DO 341 END DO 342 343 ! monotonicity algorithm 344 ! ------------------------ 345 CALL nonosc_z( ptb, ztw, zti, z2 ) 270 ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 271 END DO 272 END DO 273 END DO 274 ! 275 CALL nonosc_z( ptb, ztw, zti, z2 ) ! monotonicity algorithm 346 276 347 277 348 278 ! final trend with corrected fluxes 349 ! -----------------------------------350 279 DO jk = 1, jpkm1 351 280 DO jj = 2, jpjm1 352 281 DO ji = fs_2, fs_jpim1 ! vector opt. 353 282 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 354 ! k- vertical advective trends 355 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 356 ! add them to the general tracer trends 357 pta(ji,jj,jk) = pta(ji,jj,jk) + ztak 283 ! k- vertical advective trends added to the general tracer trends 284 pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 358 285 END DO 359 286 END DO
Note: See TracChangeset
for help on using the changeset viewer.