Changeset 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r4499 r6225 14 14 USE oce ! ocean dynamics and active tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE trdmod_oce ! ocean space and time domain 17 USE trdtra 18 USE lib_mpp 16 USE trc_oce ! share passive tracers/Ocean variables 17 USE trd_oce ! trends: ocean variables 18 USE traadv_fct ! acces to routine interp_4th_cpt 19 USE trdtra ! trends manager: tracers 20 USE diaptr ! poleward transport diagnostics 21 ! 22 USE lib_mpp ! I/O library 19 23 USE lbclnk ! ocean lateral boundary condition (or mpp link) 20 24 USE in_out_manager ! I/O manager 21 USE diaptr ! poleward transport diagnostics22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient23 USE trc_oce ! share passive tracers/Ocean variables24 25 USE wrk_nemo ! Memory Allocation 25 26 USE timing ! Timing … … 34 35 35 36 !! * Substitutions 36 # include "domzgr_substitute.h90"37 37 # include "vectopt_loop_substitute.h90" 38 38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)39 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 40 40 !! $Id$ 41 41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 43 43 CONTAINS 44 44 45 SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,&46 & ptb, ptn, pta, kjpt)45 SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 46 & ptb, ptn, pta, kjpt, kn_ubs_v ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE tra_adv_ubs *** … … 51 51 !! and add it to the general trend of passive tracer equations. 52 52 !! 53 !! ** Method : The upstream biased 3rd order scheme (UBS) is based on an53 !! ** Method : The 3rd order Upstream Biased Scheme (UBS) is based on an 54 54 !! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 55 55 !! It is only used in the horizontal direction. 56 56 !! For example the i-component of the advective fluxes are given by : 57 57 !! ! e2u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0 58 !! z wx= ! or58 !! ztu = ! or 59 59 !! ! e2u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0 60 60 !! where zltu is the second derivative of the before temperature field: 61 61 !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 62 !! This results in a dissipatively dominant (i.e. hyper-diffusive)62 !! This results in a dissipatively dominant (i.e. hyper-diffusive) 63 63 !! truncation error. The overall performance of the advection scheme 64 64 !! is similar to that reported in (Farrow and Stevens, 1995). 65 !! For stability reasons, the first term of the fluxes which corresponds65 !! For stability reasons, the first term of the fluxes which corresponds 66 66 !! to a second order centered scheme is evaluated using the now velocity 67 67 !! (centered in time) while the second term which is the diffusive part 68 68 !! of the scheme, is evaluated using the before velocity (forward in time). 69 69 !! Note that UBS is not positive. Do not use it on passive tracers. 70 !! On the vertical, the advection is evaluated using a TVD scheme, 71 !! as the UBS have been found to be too diffusive. 70 !! On the vertical, the advection is evaluated using a FCT scheme, 71 !! as the UBS have been found to be too diffusive. 72 !! kn_ubs_v argument controles whether the FCT is based on 73 !! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 74 !! scheme (kn_ubs_v=4). 72 75 !! 73 !! ** Action : - update (pta) with the now advective tracer trends 76 !! ** Action : - update pta with the now advective tracer trends 77 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 78 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 74 79 !! 75 80 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 76 81 !! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. 77 82 !!---------------------------------------------------------------------- 78 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace79 !80 83 INTEGER , INTENT(in ) :: kt ! ocean time-step index 81 84 INTEGER , INTENT(in ) :: kit000 ! first time step index 82 85 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 83 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 84 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 87 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 88 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 85 89 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean transport components 86 90 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 88 92 ! 89 93 INTEGER :: ji, jj, jk, jn ! dummy loop indices 90 REAL(wp) :: ztra, zbtr, zcoef , z2dtt! local scalars94 REAL(wp) :: ztra, zbtr, zcoef ! local scalars 91 95 REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - - 92 96 REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - - … … 96 100 IF( nn_timing == 1 ) CALL timing_start('tra_adv_ubs') 97 101 ! 98 CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw ) 99 ! 100 102 CALL wrk_alloc( jpi,jpj,jpk, ztu, ztv, zltu, zltv, zti, ztw ) 103 ! 101 104 IF( kt == kit000 ) THEN 102 105 IF(lwp) WRITE(numout,*) … … 108 111 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 109 112 ! 113 ztw (:,:, 1 ) = 0._wp ! surface & bottom value : set to zero for all tracers 114 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp 115 ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp 116 ! 110 117 ! ! =========== 111 118 DO jn = 1, kjpt ! tracer loop 112 119 ! ! =========== 113 ! 1. Bottom value : flux set to zero114 ! ----------------------------------115 zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0116 120 ! 117 DO jk = 1, jpkm1 ! Horizontal slab 118 ! 119 ! Laplacian 120 DO jj = 1, jpjm1 ! First derivative (gradient) 121 DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==! 122 DO jj = 1, jpjm1 ! First derivative (masked gradient) 121 123 DO ji = 1, fs_jpim1 ! vector opt. 122 zeeu = e2 u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)123 zeev = e1 v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)124 zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 125 zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 124 126 ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 125 127 ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 126 128 END DO 127 129 END DO 128 DO jj = 2, jpjm1 ! Second derivative (divergence)130 DO jj = 2, jpjm1 ! Second derivative (divergence) 129 131 DO ji = fs_2, fs_jpim1 ! vector opt. 130 zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )132 zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) ) 131 133 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 132 134 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef … … 134 136 END DO 135 137 ! 136 END DO ! End of slab138 END DO 137 139 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 138 139 140 ! 140 ! Horizontal advective fluxes 141 DO jk = 1, jpkm1 ! Horizontal slab 141 DO jk = 1, jpkm1 !== Horizontal advective fluxes ==! (UBS) 142 142 DO jj = 1, jpjm1 143 143 DO ji = 1, fs_jpim1 ! vector opt. 144 ! upstream transport (x2) 145 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 144 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) ! upstream transport (x2) 146 145 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 147 146 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 148 147 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 149 ! 2nd order centered advective fluxes (x2)148 ! ! 2nd order centered advective fluxes (x2) 150 149 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) 151 150 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) 152 ! UBS advective fluxes 153 zwx(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 154 zwy(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 155 END DO 156 END DO 157 END DO ! End of slab 158 159 zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends 160 161 ! Horizontal advective trends 162 DO jk = 1, jpkm1 163 ! Tracer flux divergence at t-point added to the general trend 151 ! ! UBS advective fluxes 152 ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 153 ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 154 END DO 155 END DO 156 END DO 157 ! 158 zltu(:,:,:) = pta(:,:,:,jn) ! store the initial trends before its update 159 ! 160 DO jk = 1, jpkm1 !== add the horizontal advective trend ==! 164 161 DO jj = 2, jpjm1 165 162 DO ji = fs_2, fs_jpim1 ! vector opt. 166 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 167 ! horizontal advective 168 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 169 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 170 ! add it to the general tracer trends 171 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 163 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) & 164 & - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) & 165 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 172 166 END DO 173 167 END DO 174 168 ! 175 END DO ! End of slab 176 177 ! Horizontal trend used in tra_adv_ztvd subroutine 178 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 179 180 ! 3. Save the horizontal advective trends for diagnostic 181 ! ------------------------------------------------------ 182 ! ! trend diagnostics (contribution of upstream fluxes) 183 IF( l_trd ) THEN 184 CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 185 CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 169 END DO 170 ! 171 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case 172 ! ! and/or in trend diagnostic (l_trd=T) 173 ! 174 IF( l_trd ) THEN ! trend diagnostics 175 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) ) 176 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) ) 186 177 END IF 187 178 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 188 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN189 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )190 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )179 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 180 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( ztv(:,:,:) ) 181 IF( jn == jp_sal ) str_adv(:) = ptr_sj( ztv(:,:,:) ) 191 182 ENDIF 192 193 ! TVD scheme for the vertical direction 194 ! ---------------------- 195 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 196 197 ! Bottom value : flux set to zero 198 ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0 199 200 ! Surface value 201 IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero 202 ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface 203 ENDIF 204 ! upstream advection with initial mass fluxes & intermediate update 205 ! ------------------------------------------------------------------- 206 ! Interior value 207 DO jk = 2, jpkm1 208 DO jj = 1, jpj 209 DO ji = 1, jpi 210 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 211 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 212 ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 213 END DO 214 END DO 215 END DO 216 ! update and guess with monotonic sheme 217 DO jk = 1, jpkm1 218 z2dtt = p2dt(jk) 219 DO jj = 2, jpjm1 220 DO ji = fs_2, fs_jpim1 ! vector opt. 221 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 222 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 223 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 224 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 225 END DO 226 END DO 227 END DO 228 ! 229 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 230 231 ! antidiffusive flux : high order minus low order 232 ztw(:,:,1) = 0.e0 ! Surface value 233 DO jk = 2, jpkm1 ! Interior value 234 DO jj = 1, jpj 235 DO ji = 1, jpi 236 ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) 237 END DO 238 END DO 239 END DO 240 ! 241 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 242 243 ! final trend with corrected fluxes 244 DO jk = 1, jpkm1 183 ! 184 ! !== vertical advective trend ==! 185 ! 186 SELECT CASE( kn_ubs_v ) ! select the vertical advection scheme 187 ! 188 CASE( 2 ) ! 2nd order FCT 189 ! 190 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 191 ! 192 ! !* upstream advection with initial mass fluxes & intermediate update ==! 193 DO jk = 2, jpkm1 ! Interior value (w-masked) 194 DO jj = 1, jpj 195 DO ji = 1, jpi 196 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 197 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 198 ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 199 END DO 200 END DO 201 END DO 202 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as ztw has been w-masked) 203 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 207 END DO 208 END DO 209 ELSE ! no cavities: only at the ocean surface 210 ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 211 ENDIF 212 ENDIF 213 ! 214 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 215 DO jj = 2, jpjm1 216 DO ji = fs_2, fs_jpim1 ! vector opt. 217 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 218 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 219 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 220 END DO 221 END DO 222 END DO 223 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 224 ! 225 ! !* anti-diffusive flux : high order minus low order 226 DO jk = 2, jpkm1 ! Interior value (w-masked) 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ztw(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 230 & - ztw(ji,jj,jk) ) * wmask(ji,jj,jk) 231 END DO 232 END DO 233 END DO 234 ! ! top ocean value: high order == upstream ==>> zwz=0 235 IF( ln_linssh ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 236 ! 237 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 238 ! 239 CASE( 4 ) ! 4th order COMPACT 240 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! 4th order compact interpolation of T at w-point 241 DO jk = 2, jpkm1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 244 ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 245 END DO 246 END DO 247 END DO 248 IF( ln_linssh ) ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work 249 ! 250 END SELECT 251 ! 252 DO jk = 1, jpkm1 ! final trend with corrected fluxes 245 253 DO jj = 2, jpjm1 246 254 DO ji = fs_2, fs_jpim1 ! vector opt. 247 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 248 ! k- vertical advective trends 249 ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) 250 ! added to the general tracer trends 251 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 252 END DO 253 END DO 254 END DO 255 256 ! Save the final vertical advective trends 257 IF( l_trd ) THEN ! vertical advective trend diagnostics 255 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 256 END DO 257 END DO 258 END DO 259 ! 260 IF( l_trd ) THEN ! vertical advective trend diagnostics 258 261 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 259 262 DO jj = 2, jpjm1 260 263 DO ji = fs_2, fs_jpim1 ! vector opt. 261 z btr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )262 z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr263 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn264 END DO 265 END DO 266 END DO 267 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zltv )264 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) & 265 & + ptn(ji,jj,jk,jn) * ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) & 266 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 267 END DO 268 END DO 269 END DO 270 CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv ) 268 271 ENDIF 269 272 ! 270 END DO271 ! 272 CALL wrk_dealloc( jpi, jpj, jpk,ztu, ztv, zltu, zltv, zti, ztw )273 END DO 274 ! 275 CALL wrk_dealloc( jpi,jpj,jpk, ztu, ztv, zltu, zltv, zti, ztw ) 273 276 ! 274 277 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_ubs') … … 290 293 !! in-space based differencing for fluid 291 294 !!---------------------------------------------------------------------- 292 ! 293 REAL(wp), INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step 295 REAL(wp), INTENT(in ) :: p2dt ! tracer time-step 294 296 REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field 295 297 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field … … 298 300 INTEGER :: ji, jj, jk ! dummy loop indices 299 301 INTEGER :: ikm1 ! local integer 300 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn , z2dtt! local scalars302 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 301 303 REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo 302 304 !!---------------------------------------------------------------------- … … 304 306 IF( nn_timing == 1 ) CALL timing_start('nonosc_z') 305 307 ! 306 CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo ) 307 ! 308 308 CALL wrk_alloc( jpi,jpj,jpk, zbetup, zbetdo ) 309 ! 309 310 zbig = 1.e+40_wp 310 311 zrtrn = 1.e-15_wp 311 312 zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp 312 313 ! 313 314 ! Search local extrema 314 315 ! -------------------- 315 ! large negative value (-zbig) inside land316 ! ! large negative value (-zbig) inside land 316 317 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 317 318 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 318 ! search maximum in neighbourhood319 DO jk = 1, jpkm1 319 ! 320 DO jk = 1, jpkm1 ! search maximum in neighbourhood 320 321 ikm1 = MAX(jk-1,1) 321 322 DO jj = 2, jpjm1 … … 327 328 END DO 328 329 END DO 329 ! large positive value (+zbig) inside land330 ! ! large positive value (+zbig) inside land 330 331 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 331 332 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 332 ! search minimum in neighbourhood333 DO jk = 1, jpkm1 333 ! 334 DO jk = 1, jpkm1 ! search minimum in neighbourhood 334 335 ikm1 = MAX(jk-1,1) 335 336 DO jj = 2, jpjm1 … … 341 342 END DO 342 343 END DO 343 344 ! restore masked values to zero 344 ! ! restore masked values to zero 345 345 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 346 346 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 347 348 349 ! 2. Positive and negative part of fluxes and beta terms 350 ! ------------------------------------------------------ 351 347 ! 348 ! Positive and negative part of fluxes and beta terms 349 ! --------------------------------------------------- 352 350 DO jk = 1, jpkm1 353 z2dtt = p2dt(jk)354 351 DO jj = 2, jpjm1 355 352 DO ji = fs_2, fs_jpim1 ! vector opt. … … 358 355 zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 359 356 ! up & down beta terms 360 zbt = e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt357 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 361 358 zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 362 359 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt … … 364 361 END DO 365 362 END DO 363 ! 366 364 ! monotonic flux in the k direction, i.e. pcc 367 365 ! ------------------------------------------- … … 377 375 END DO 378 376 ! 379 CALL wrk_dealloc( jpi, jpj, jpk,zbetup, zbetdo )377 CALL wrk_dealloc( jpi,jpj,jpk, zbetup, zbetdo ) 380 378 ! 381 379 IF( nn_timing == 1 ) CALL timing_stop('nonosc_z')
Note: See TracChangeset
for help on using the changeset viewer.
