Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
- Property svn:keywords set to Id
r1559 r2528 2 2 !!============================================================================== 3 3 !! *** MODULE traadv_qck *** 4 !! Ocean activetracers: horizontal & vertical advective trend4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 6 !! History : 3.0 ! 2008-07 (G. Reffray) Original code 7 !! 3.3 ! 2010-05 (C.Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 7 8 !!---------------------------------------------------------------------- 8 9 9 10 !!---------------------------------------------------------------------- 10 !! tra_adv_qck 11 !! 12 !! tra_adv_qck_i : 13 !! tra_adv_qck_j : 11 !! tra_adv_qck : update the tracer trend with the horizontal advection 12 !! trends using a 3rd order finite difference scheme 13 !! tra_adv_qck_i : apply QUICK scheme in i-direction 14 !! tra_adv_qck_j : apply QUICK scheme in j-direction 14 15 !! tra_adv_cen2_k : 2nd centered scheme for the vertical advection 15 16 !!---------------------------------------------------------------------- 16 17 USE oce ! ocean dynamics and active tracers 17 18 USE dom_oce ! ocean space and time domain 18 USE trdmod ! ocean active tracers trends19 USE trd mod_oce ! ocean variables trends19 USE trdmod_oce ! ocean space and time domain 20 USE trdtra ! ocean tracers trends 20 21 USE trabbl ! advective term in the BBL 21 22 USE lib_mpp ! distribued memory computing … … 24 25 USE in_out_manager ! I/O manager 25 26 USE diaptr ! poleward transport diagnostics 26 USE prtctl ! Print control27 USE trc_oce ! share passive tracers/Ocean variables 27 28 28 29 IMPLICIT NONE … … 31 32 PUBLIC tra_adv_qck ! routine called by step.F90 32 33 33 REAL(wp), DIMENSION(jpi,jpj) :: btr234 REAL(wp) :: r1_634 LOGICAL :: l_trd ! flag to compute trends 35 REAL(wp) :: r1_6 = 1./ 6. ! 1/6 ratio 35 36 36 37 !! * Substitutions … … 38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3. 2 , LOCEAN-IPSL (2009)41 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 41 42 !! $Id$ 42 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 44 !!---------------------------------------------------------------------- 44 45 45 CONTAINS 46 46 47 SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn ) 47 SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn, & 48 & ptb, ptn, pta, kjpt ) 48 49 !!---------------------------------------------------------------------- 49 50 !! *** ROUTINE tra_adv_qck *** … … 69 70 !! dt = 2*rdtra and the scalar values are tb and sb 70 71 !! 71 !! On the vertical, the simple centered scheme used tn and sn72 !! On the vertical, the simple centered scheme used ptn 72 73 !! 73 74 !! The fluxes are bounded by the ULTIMATE limiter to … … 75 76 !! prevent the appearance of spurious numerical oscillations 76 77 !! 77 !! ** Action : - update ( ta,sa) with the now advective tracer trends78 !! - save the trends ('key_trdtra')78 !! ** Action : - update (pta) with the now advective tracer trends 79 !! - save the trends 79 80 !! 80 81 !! ** Reference : Leonard (1979, 1991) 81 82 !!---------------------------------------------------------------------- 82 USE oce, ONLY : ztrdt => ua ! use ua as workspace 83 USE oce, ONLY : ztrds => va ! use va as workspace 84 !! 85 INTEGER , INTENT(in) :: kt ! ocean time-step index 86 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component 87 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component 88 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component 89 !! 90 INTEGER :: ji, jj, jk ! dummy loop indices 91 REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! temporary scalars 92 REAL(wp) :: zbtr, z2 ! " " 93 !!---------------------------------------------------------------------- 94 95 IF( kt == nit000 ) THEN 83 INTEGER , INTENT(in ) :: kt ! ocean time-step index 84 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 85 INTEGER , INTENT(in ) :: kjpt ! number of tracers 86 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 87 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 88 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 89 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 90 !!---------------------------------------------------------------------- 91 92 IF( kt == nit000 ) THEN 96 93 IF(lwp) WRITE(numout,*) 97 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme '94 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype 98 95 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 99 96 IF(lwp) WRITE(numout,*) 100 btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 101 r1_6 = 1. / 6. 97 ! 98 l_trd = .FALSE. 99 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 102 100 ENDIF 103 101 104 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1.105 ELSE ; z2 = 2.106 ENDIF107 108 102 ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 109 !--------------------------------------------------------------------------- 110 111 CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2) 112 CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2) 113 114 IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 115 116 CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2) 117 CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2) 118 119 IF( l_trdtra ) THEN 120 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 121 ! 122 ztrdt(:,:,:) = ta(:,:,:) ! Save the horizontal up-to-date ta/sa trends 123 ztrds(:,:,:) = sa(:,:,:) 124 END IF 125 126 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck had - Ta: ', mask1=tmask, & 127 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 103 CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 104 CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 128 105 129 106 ! II. The vertical fluxes are computed with the 2nd order centered scheme 130 !------------------------------------------------------------------------- 131 ! 132 CALL tra_adv_cen2_k( pwn, tn, ta ) 133 CALL tra_adv_cen2_k( pwn, sn, sa ) 134 ! 135 !Save the vertical advective trends for diagnostic 136 ! ---------------------------------------------------- 137 IF( l_trdtra ) THEN 138 ! Recompute the vertical advection zta & zsa trends computed 139 ! at the step 2. above in making the difference between the new 140 ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract 141 ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 142 143 DO jk = 1, jpkm1 144 DO jj = 2, jpjm1 145 DO ji = fs_2, fs_jpim1 ! vector opt. 146 #if defined key_zco 147 zbtr = btr2(ji,jj) 148 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 149 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 150 #else 151 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 152 z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk) 153 z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk) 154 #endif 155 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 156 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 157 ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 158 END DO 159 END DO 160 END DO 161 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) 162 ENDIF 163 164 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck zad - Ta: ', mask1=tmask, & 165 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 107 CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt ) 166 108 ! 167 109 END SUBROUTINE tra_adv_qck 168 110 169 111 170 SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 ) 171 !!---------------------------------------------------------------------- 172 !! 173 !!---------------------------------------------------------------------- 174 REAL, INTENT(in) :: z2 175 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran ! horizontal effective velocity 176 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra 177 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa 178 ! 179 INTEGER :: ji, jj, jk 180 REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars 181 REAL(wp) :: z_hdivn_x 182 REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl 183 REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e 184 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux 112 SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, & 113 & ptb, ptn, pta, kjpt ) 114 !!---------------------------------------------------------------------- 115 !! 116 !!---------------------------------------------------------------------- 117 USE oce , zwx => ua ! use ua as workspace 118 !! 119 INTEGER , INTENT(in ) :: kt ! ocean time-step index 120 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 121 INTEGER , INTENT(in ) :: kjpt ! number of tracers 122 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 123 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components 124 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 125 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 126 !! 127 INTEGER :: ji, jj, jk, jn ! dummy loop indices 128 REAL(wp) :: ztra, zbtr ! local scalars 129 REAL(wp) :: zdir, zdx, zdt, zmsk ! local scalars 130 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfc, zfd ! 3D wokspace 185 131 !---------------------------------------------------------------------- 186 132 187 zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0 188 zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0 189 zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0 190 zfho (:,jpj) = 0.e0 191 ! =============== 192 DO jk = 1, jpkm1 ! Horizontal slab 193 ! ! =============== 194 !--- Computation of the ustream and downstream value of the tracer and the mask 195 DO jj = 2, jpjm1 196 DO ji = 2, fs_jpim1 ! vector opt. 197 ! Upstream in the x-direction for the tracer 198 zupst(ji,jj)=tra(ji-1,jj,jk) 199 ! Downstream in the x-direction for the tracer 200 zdwst(ji,jj)=tra(ji+1,jj,jk) 201 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 202 zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2 203 END DO 204 END DO 205 ! 206 !--- Lateral boundary conditions 207 CALL lbc_lnk( zupst(:,:), 'T', 1. ) 208 CALL lbc_lnk( zdwst(:,:), 'T', 1. ) 209 CALL lbc_lnk( zmask(:,:), 'T', 1. ) 133 ! ! =========== 134 DO jn = 1, kjpt ! tracer loop 135 ! ! =========== 136 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 137 zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0 138 ! 139 DO jk = 1, jpkm1 140 ! 141 !--- Computation of the ustream and downstream value of the tracer and the mask 142 DO jj = 2, jpjm1 143 DO ji = fs_2, fs_jpim1 ! vector opt. 144 ! Upstream in the x-direction for the tracer 145 zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) 146 ! Downstream in the x-direction for the tracer 147 zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) 148 END DO 149 END DO 150 END DO 151 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions 152 210 153 ! 211 154 ! Horizontal advective fluxes 212 155 ! --------------------------- 213 156 ! 214 dt = z2 * rdttra(jk) 215 !--- tracer flux at u-points 216 DO jj = 1, jpjm1 217 DO ji = 1, jpi 218 #if defined key_zco 219 zsc_e(ji,jj) = e2u(ji,jj) 220 #else 221 zsc_e(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) 222 #endif 223 dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : dir = 1 otherwise dir = 0 224 dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj) 225 zc_cfl (ji,jj) = ABS(pun(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on x-direction) 226 227 zfu(ji,jj) = dir*zupst(ji ,jj )+(1-dir)*zdwst(ji+1,jj ) ! FU in the x-direction for T 228 zfc(ji,jj) = dir*tra (ji ,jj,jk)+(1-dir)*tra (ji+1,jj,jk) ! FC in the x-direction for T 229 zfd(ji,jj) = dir*tra (ji+1,jj,jk)+(1-dir)*tra (ji ,jj,jk) ! FD in the x-direction for T 230 zmskl(ji,jj) = dir*zmask(ji ,jj) +(1-dir)*zmask(ji+1,jj) 231 END DO 232 END DO 233 ! 157 DO jk = 1, jpkm1 158 DO jj = 2, jpjm1 159 DO ji = fs_2, fs_jpim1 ! vector opt. 160 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 161 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 162 END DO 163 END DO 164 END DO 165 ! 166 DO jk = 1, jpkm1 167 zdt = p2dt(jk) 168 DO jj = 2, jpjm1 169 DO ji = fs_2, fs_jpim1 ! vector opt. 170 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 171 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk) 172 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 173 zfc(ji,jj,jk) = zdir * ptb(ji ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn) ! FC in the x-direction for T 174 zfd(ji,jj,jk) = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji ,jj,jk,jn) ! FD in the x-direction for T 175 END DO 176 END DO 177 END DO 178 !--- Lateral boundary conditions 179 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 180 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) 181 234 182 !--- QUICKEST scheme 183 CALL quickest( zfu, zfd, zfc, zwx ) 184 ! 185 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 186 DO jk = 1, jpkm1 187 DO jj = 2, jpjm1 188 DO ji = fs_2, fs_jpim1 ! vector opt. 189 zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 190 ENDDO 191 END DO 192 END DO 193 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions 194 195 ! 235 196 ! Tracer flux on the x-direction 236 CALL quickest(zfu,zfd,zfc,zfho,zc_cfl) 237 !--- If the second ustream point is a land point 238 !--- the flux is computed by the 1st order UPWIND scheme 239 zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:) 240 !--- Computation of fluxes 241 zflux(:,:,jk) = zsc_e(:,:)*pun(:,:,jk)*zfho(:,:) 242 ! 243 !--- Tracer flux divergence at t-point added to the general trend 244 DO jj = 2, jpjm1 245 DO ji = fs_2, fs_jpim1 ! vector opt. 246 !--- horizontal advective trends 247 #if defined key_zco 248 zbtr = btr2(ji,jj) 249 #else 250 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 251 #endif 252 za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) 253 !--- add it to the general tracer trends 254 traa(ji,jj,jk) = traa(ji,jj,jk) + za 255 END DO 256 END DO 257 ! ! =============== 258 END DO ! End of slab 259 ! ! =============== 197 DO jk = 1, jpkm1 198 ! 199 DO jj = 2, jpjm1 200 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 202 !--- If the second ustream point is a land point 203 !--- the flux is computed by the 1st order UPWIND scheme 204 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 205 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 206 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk) 207 END DO 208 END DO 209 ! 210 ! Computation of the trend 211 DO jj = 2, jpjm1 212 DO ji = fs_2, fs_jpim1 ! vector opt. 213 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 214 ! horizontal advective trends 215 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 216 !--- add it to the general tracer trends 217 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 218 END DO 219 END DO 220 ! 221 END DO 222 ! ! trend diagnostics (contribution of upstream fluxes) 223 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 224 ! 225 END DO 260 226 ! 261 ! Save the horizontal advective trends for diagnostic 262 ! ----------------------------------------------------- 263 IF( l_trdtra ) THEN 264 ! T/S ZONAL advection trends 265 ztrdtra(:,:,:) = 0.e0 266 ! 267 DO jk = 1, jpkm1 227 END SUBROUTINE tra_adv_qck_i 228 229 230 SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, & 231 & ptb, ptn, pta, kjpt ) 232 !!---------------------------------------------------------------------- 233 !! 234 !!---------------------------------------------------------------------- 235 USE oce , zwy => ua ! use ua as workspace 236 !! 237 INTEGER , INTENT(in ) :: kt ! ocean time-step index 238 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 239 INTEGER , INTENT(in ) :: kjpt ! number of tracers 240 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 241 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components 242 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 243 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 244 !! 245 INTEGER :: ji, jj, jk, jn ! dummy loop indices 246 REAL(wp) :: ztra, zbtr ! local scalars 247 REAL(wp) :: zdir, zdx, zdt, zmsk ! local scalars 248 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfc, zfd ! 3D wokspace 249 !---------------------------------------------------------------------- 250 251 ! ! =========== 252 DO jn = 1, kjpt ! tracer loop 253 ! ! =========== 254 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 255 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 256 ! 257 DO jk = 1, jpkm1 258 ! 259 !--- Computation of the ustream and downstream value of the tracer and the mask 268 260 DO jj = 2, jpjm1 269 261 DO ji = fs_2, fs_jpim1 ! vector opt. 270 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 271 ! N.B. This computation is not valid along OBCs (if any) 272 #if defined key_zco 273 zbtr = btr2(ji,jj) 274 z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & 275 & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 276 #else 277 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 278 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 279 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 280 #endif 281 ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) + tran(ji,jj,jk) * z_hdivn_x 282 END DO 283 END DO 284 END DO 285 END IF 286 287 END SUBROUTINE tra_adv_qck_i 288 289 290 SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 ) 291 !!---------------------------------------------------------------------- 292 !! 293 !!---------------------------------------------------------------------- 294 INTEGER, INTENT(in) :: kt ! ocean time-step index 295 REAL, INTENT(in) :: z2 296 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pvn, tra, tran ! horizontal effective velocity 297 REAL(wp), INTENT(out) , DIMENSION(jpj) :: trd_adv 298 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra 299 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa 300 !! 301 INTEGER :: ji, jj, jk 302 REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars 303 REAL(wp) :: z_hdivn_y 304 REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl 305 REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e 306 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux 307 !---------------------------------------------------------------------- 308 ! II. Part 2 : y-direction 309 !---------------------------------------------------------------------- 310 311 zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0 312 zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0 313 zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0 314 zfho (:,jpj) = 0.e0 315 316 ! =============== 317 DO jk = 1, jpkm1 ! Horizontal slab 318 ! ! =============== 319 !--- Computation of the ustream and downstream value of the tracer and the mask 320 DO jj = 2, jpjm1 321 DO ji = 2, fs_jpim1 ! vector opt. 322 ! Upstream in the x-direction for the tracer 323 zupst(ji,jj)=tra(ji,jj-1,jk) 324 ! Downstream in the x-direction for the tracer 325 zdwst(ji,jj)=tra(ji,jj+1,jk) 326 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 327 zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 328 END DO 329 END DO 330 ! 331 !--- Lateral boundary conditions 332 CALL lbc_lnk( zupst(:,:), 'T', 1. ) 333 CALL lbc_lnk( zdwst(:,:), 'T', 1. ) 334 CALL lbc_lnk( zmask(:,:), 'T', 1. ) 262 ! Upstream in the x-direction for the tracer 263 zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn) 264 ! Downstream in the x-direction for the tracer 265 zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn) 266 END DO 267 END DO 268 END DO 269 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions 270 271 335 272 ! 336 273 ! Horizontal advective fluxes 337 274 ! --------------------------- 338 275 ! 339 dt = z2 * rdttra(jk) 340 !--- tracer flux at v-points 341 DO jj = 1, jpjm1 342 DO ji = 1, jpi 343 #if defined key_zco 344 zsc_e(ji,jj) = e1v(ji,jj) 345 #else 346 zsc_e(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) 347 #endif 348 dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dir = 1 otherwise dir = 0 349 dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1) 350 zc_cfl(ji,jj) = ABS(pvn(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on y-direction) 351 352 zfu(ji,jj) = dir*zupst(ji,jj )+(1-dir)*zdwst(ji,jj+1 ) ! FU in the y-direction for T 353 zfc(ji,jj) = dir*tra (ji,jj ,jk)+(1-dir)*tra (ji,jj+1,jk) ! FC in the y-direction for T 354 zfd(ji,jj) = dir*tra (ji,jj+1,jk)+(1-dir)*tra (ji,jj ,jk) ! FD in the y-direction for T 355 zmskl(ji,jj) = dir*zmask(ji,jj )+(1-dir)*zmask(ji,jj+1) 356 END DO 357 END DO 358 ! 276 DO jk = 1, jpkm1 277 DO jj = 2, jpjm1 278 DO ji = fs_2, fs_jpim1 ! vector opt. 279 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 280 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 281 END DO 282 END DO 283 END DO 284 ! 285 DO jk = 1, jpkm1 286 zdt = p2dt(jk) 287 DO jj = 2, jpjm1 288 DO ji = fs_2, fs_jpim1 ! vector opt. 289 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 290 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk) 291 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 292 zfc(ji,jj,jk) = zdir * ptb(ji,jj ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn) ! FC in the x-direction for T 293 zfd(ji,jj,jk) = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj ,jk,jn) ! FD in the x-direction for T 294 END DO 295 END DO 296 END DO 297 298 !--- Lateral boundary conditions 299 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 300 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) 301 359 302 !--- QUICKEST scheme 360 ! Tracer flux on the y-direction 361 CALL quickest(zfu,zfd,zfc,zfho,zc_cfl) 362 !--- If the second ustream point is a land point 363 !--- the flux is computed by the 1st order UPWIND scheme 364 zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:) 365 !--- Computation of fluxes 366 zflux(:,:,jk) = zsc_e(:,:)*pvn(:,:,jk)*zfho(:,:) 367 ! 368 !--- Tracer flux divergence at t-point added to the general trend 369 DO jj = 2, jpjm1 370 DO ji = fs_2, fs_jpim1 ! vector opt. 371 !--- horizontal advective trends 372 #if defined key_zco 373 zbtr = btr2(ji,jj) 374 #else 375 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 376 #endif 377 za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) 378 !--- add it to the general tracer trends 379 traa(ji,jj,jk) = traa(ji,jj,jk) + za 380 END DO 381 END DO 382 ! ! =============== 383 END DO ! End of slab 384 ! ! =============== 385 ! 386 ! Save the horizontal advective trends for diagnostic 387 ! ----------------------------------------------------- 388 IF( l_trdtra ) THEN 389 ! T/S MERIDIONAL advection trends 390 DO jk = 1, jpkm1 391 DO jj = 2, jpjm1 392 DO ji = fs_2, fs_jpim1 ! vector opt. 393 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 394 ! N.B. This computation is not valid along OBCs (if any) 395 #if defined key_zco 396 zbtr = btr2(ji,jj) 397 z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) & 398 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 399 #else 400 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 401 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 402 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 403 #endif 404 ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) + tran(ji,jj,jk) * z_hdivn_y 405 END DO 406 END DO 407 END DO 408 END IF 409 410 ! "zonal" mean advective heat and salt transport 411 ! ---------------------------------------------- 412 413 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 414 IF( lk_zco ) THEN 415 DO jk = 1, jpkm1 416 DO jj = 2, jpjm1 417 DO ji = fs_2, fs_jpim1 ! vector opt. 418 zflux(ji,jj,jk) = zflux(ji,jj,jk) * fse3v(ji,jj,jk) 419 END DO 420 END DO 421 END DO 303 CALL quickest( zfu, zfd, zfc, zwy ) 304 ! 305 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 306 DO jk = 1, jpkm1 307 DO jj = 2, jpjm1 308 DO ji = fs_2, fs_jpim1 ! vector opt. 309 zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 310 END DO 311 END DO 312 END DO 313 !--- Lateral boundary conditions 314 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 315 ! 316 ! Tracer flux on the x-direction 317 DO jk = 1, jpkm1 318 ! 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 322 !--- If the second ustream point is a land point 323 !--- the flux is computed by the 1st order UPWIND scheme 324 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 325 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 326 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk) 327 END DO 328 END DO 329 ! 330 ! Computation of the trend 331 DO jj = 2, jpjm1 332 DO ji = fs_2, fs_jpim1 ! vector opt. 333 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 334 ! horizontal advective trends 335 ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 336 !--- add it to the general tracer trends 337 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 338 END DO 339 END DO 340 ! 341 END DO 342 ! ! trend diagnostics (contribution of upstream fluxes) 343 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 344 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 345 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 346 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 347 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 422 348 ENDIF 423 trd_adv(:) = ptr_vj( zflux(:,:,:) ) 424 ENDIF 425 426 END SUBROUTINE tra_adv_qck_j 427 428 429 SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta ) 430 !!---------------------------------------------------------------------- 431 !! 432 !!---------------------------------------------------------------------- 433 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical effective velocity 434 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now tracer 435 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer general trend 436 !! 437 INTEGER :: ji, jj, jk ! dummy loop indices 438 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux ! 3D workspace 439 !!---------------------------------------------------------------------- 440 ! 441 ! !== Vertical advective fluxes ==! 442 zflux(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 443 ! 444 ! ! Surface value 445 IF( lk_vvl ) THEN ; zflux(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 446 ELSE ; zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! Constant volume : advective flux through the surface 447 ENDIF 448 ! 449 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point 450 DO jj = 2, jpjm1 451 DO ji = fs_2, fs_jpim1 ! vector opt. 452 zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) ) 453 END DO 454 END DO 349 ! 455 350 END DO 456 351 ! 457 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 458 DO jj = 2, jpjm1 459 DO ji = fs_2, fs_jpim1 ! vector opt. 460 pta(ji,jj,jk) = pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) ) & 461 & / fse3t(ji,jj,jk) 462 END DO 463 END DO 352 END SUBROUTINE tra_adv_qck_j 353 354 355 SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, & 356 & ptn, pta, kjpt ) 357 !!---------------------------------------------------------------------- 358 !! 359 !!---------------------------------------------------------------------- 360 USE oce , zwz => ua ! use ua as workspace 361 !! 362 INTEGER , INTENT(in ) :: kt ! ocean time-step index 363 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 364 INTEGER , INTENT(in ) :: kjpt ! number of tracers 365 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity 366 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! before and now tracer fields 367 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 368 !! 369 INTEGER :: ji, jj, jk, jn ! dummy loop indices 370 REAL(wp) :: zbtr , ztra ! temporary scalars 371 !!---------------------------------------------------------------------- 372 373 ! ! =========== 374 DO jn = 1, kjpt ! tracer loop 375 ! ! =========== 376 ! 1. Bottom value : flux set to zero 377 zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 378 ! 379 ! ! Surface value 380 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 381 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface 382 ENDIF 383 ! 384 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point 385 DO jj = 2, jpjm1 386 DO ji = fs_2, fs_jpim1 ! vector opt. 387 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 388 END DO 389 END DO 390 END DO 391 ! 392 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 393 DO jj = 2, jpjm1 394 DO ji = fs_2, fs_jpim1 ! vector opt. 395 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 396 ! k- vertical advective trends 397 ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 398 ! added to the general tracer trends 399 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 400 END DO 401 END DO 402 END DO 403 ! ! Save the vertical advective trends for diagnostic 404 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) ) 405 ! 464 406 END DO 465 407 ! … … 467 409 468 410 469 SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl ) 470 !!---------------------------------------------------------------------- 471 !! 472 !!---------------------------------------------------------------------- 473 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl 474 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho 475 REAL(wp) , DIMENSION(jpi,jpj) :: zcurv, zcoef1, zcoef2, zcoef3 ! temporary scalars 476 ! 477 zcurv (:,:) = fd(:,:) + fu(:,:) - 2.*fc(:,:) 478 zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) ) 479 zcoef2(:,:) = 0.5*fc_cfl(:,:)*( fd(:,:) - fc(:,:) ) 480 zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*r1_6 )*zcurv(:,:) 481 fho (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:) ! phi_f QUICKEST 482 ! 483 zcoef1(:,:) = fd(:,:) - fu(:,:) ! DEL 484 zcoef2(:,:) = ABS( zcoef1(:,:) ) ! ABS(DEL) 485 zcoef3(:,:) = ABS( zcurv(:,:) ) ! ABS(CURV) 486 ! 487 WHERE ( zcoef3(:,:) >= zcoef2(:,:) ) 488 fho(:,:) = fc(:,:) 489 ELSEWHERE 490 zcoef3(:,:) = fu(:,:) + ( ( fc(:,:) - fu(:,:) )/MAX(fc_cfl(:,:),1.e-9) ) ! phi_REF 491 WHERE ( zcoef1(:,:) >= 0.e0 ) 492 fho(:,:) = MAX(fc(:,:),fho(:,:)) 493 fho(:,:) = MIN(fho(:,:),MIN(zcoef3(:,:),fd(:,:))) 494 ELSEWHERE 495 fho(:,:) = MIN(fc(:,:),fho(:,:)) 496 fho(:,:) = MAX(fho(:,:),MAX(zcoef3(:,:),fd(:,:))) 497 ENDWHERE 498 ENDWHERE 411 SUBROUTINE quickest( pfu, pfd, pfc, puc ) 412 !!---------------------------------------------------------------------- 413 !! 414 !! ** Purpose : Computation of advective flux with Quickest scheme 415 !! 416 !! ** Method : 417 !!---------------------------------------------------------------------- 418 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point 419 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point 420 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point) 421 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux 422 !! 423 INTEGER :: ji, jj, jk ! dummy loop indices 424 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars 425 REAL(wp) :: zc, zcurv, zfho ! - - 426 !---------------------------------------------------------------------- 427 428 DO jk = 1, jpkm1 429 DO jj = 1, jpj 430 DO ji = 1, jpi 431 zc = puc(ji,jj,jk) ! Courant number 432 zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) 433 zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) ) 434 zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 435 zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 436 zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST 437 ! 438 zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) 439 zcoef2 = ABS( zcoef1 ) 440 zcoef3 = ABS( zcurv ) 441 IF( zcoef3 >= zcoef2 ) THEN 442 zfho = pfc(ji,jj,jk) 443 ELSE 444 zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF 445 IF( zcoef1 >= 0. ) THEN 446 zfho = MAX( pfc(ji,jj,jk), zfho ) 447 zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 448 ELSE 449 zfho = MIN( pfc(ji,jj,jk), zfho ) 450 zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 451 ENDIF 452 ENDIF 453 puc(ji,jj,jk) = zfho 454 END DO 455 END DO 456 END DO 499 457 ! 500 458 END SUBROUTINE quickest
Note: See TracChangeset
for help on using the changeset viewer.
