- Timestamp:
- 2017-12-01T18:44:09+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r8698 r8882 4 4 !! Ocean active tracers: vertical component of the tracer mixing trend 5 5 !!============================================================================== 6 !! History : 1.0 ! 2005-11 (G. Madec) Original code 7 !! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA 6 !! History : 1.0 ! 2005-11 (G. Madec) Original code 7 !! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA 8 !! 4.0 ! 2017-06 (G. Madec) remove explict time-stepping option 8 9 !!---------------------------------------------------------------------- 9 10 10 11 !!---------------------------------------------------------------------- 11 12 !! tra_zdf : Update the tracer trend with the vertical diffusion 12 !! tra_zdf_init : initialisation of the computation13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and tracers variables … … 20 20 USE ldftra ! lateral diffusion: eddy diffusivity 21 21 USE ldfslp ! lateral diffusion: iso-neutral slope 22 USE trazdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine)23 USE trazdf_imp ! vertical diffusion: implicit (tra_zdf_imp routine)24 22 USE trd_oce ! trends: ocean variables 25 23 USE trdtra ! trends: tracer trend manager … … 29 27 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 30 28 USE lib_mpp ! MPP library 31 USE wrk_nemo ! Memory allocation32 29 USE timing ! Timing 33 30 … … 35 32 PRIVATE 36 33 37 PUBLIC tra_zdf ! routine called by step.F90 38 PUBLIC tra_zdf_init ! routine called by nemogcm.F90 39 40 INTEGER :: nzdf = 0 ! type vertical diffusion algorithm used (defined from ln_zdf... namlist logicals) 34 PUBLIC tra_zdf ! called by step.F90 35 PUBLIC tra_zdf_imp ! called by trczdf.F90 41 36 42 37 !! * Substitutions 43 # include "zdfddm_substitute.h90"44 38 # include "vectopt_loop_substitute.h90" 45 39 !!---------------------------------------------------------------------- 46 !! NEMO/OPA 3.7 , NEMO Consortium (2015)40 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 47 41 !! $Id$ 48 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 56 50 !! ** Purpose : compute the vertical ocean tracer physics. 57 51 !!--------------------------------------------------------------------- 58 INTEGER, INTENT( in ) :: kt ! ocean time-step index 59 ! 60 INTEGER :: jk ! Dummy loop indices 61 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds ! 3D workspace 62 !!--------------------------------------------------------------------- 63 ! 64 IF( nn_timing == 1 ) CALL timing_start('tra_zdf') 65 ! 66 IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000 67 r2dt = rdt ! = rdt (restarting with Euler time stepping) 68 ELSEIF( kt <= nit000 + 1) THEN ! at nit000 or nit000+1 69 r2dt = 2. * rdt ! = 2 rdt (leapfrog) 70 ENDIF 71 ! 72 IF( l_trdtra ) THEN !* Save ta and sa trends 73 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 53 ! 54 INTEGER :: jk ! Dummy loop indices 55 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace 56 !!--------------------------------------------------------------------- 57 ! 58 IF( ln_timing ) CALL timing_start('tra_zdf') 59 ! 60 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000, = rdt (restarting with Euler time stepping) 61 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! otherwise, = 2 rdt (leapfrog) 62 ENDIF 63 ! 64 IF( l_trdtra ) THEN !* Save ta and sa trends 65 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 74 66 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 75 67 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 76 68 ENDIF 77 69 ! 78 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend 79 CASE ( 0 ) ; CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts ) ! explicit scheme 80 CASE ( 1 ) ; CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) ! implicit scheme 81 END SELECT 70 ! !* compute lateral mixing trend and add it to the general trend 71 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) 72 82 73 !!gm WHY here ! and I don't like that ! 83 74 ! DRAKKAR SSS control { … … 90 81 DO jk = 1, jpkm1 91 82 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) & 92 &/ (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)83 & / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 93 84 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) & 94 &/ (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)85 & / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 95 86 END DO 96 87 !!gm this should be moved in trdtra.F90 and done on all trends … … 100 91 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 101 92 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 102 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )93 DEALLOCATE( ztrdt , ztrds ) 103 94 ENDIF 104 95 ! ! print mean trends (used for debugging) … … 106 97 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 107 98 ! 108 IF( nn_timing == 1 )CALL timing_stop('tra_zdf')99 IF( ln_timing ) CALL timing_stop('tra_zdf') 109 100 ! 110 101 END SUBROUTINE tra_zdf 111 102 112 113 SUBROUTINE tra_zdf_i nit103 104 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 114 105 !!---------------------------------------------------------------------- 115 !! *** ROUTINE tra_zdf_init *** 116 !! 117 !! ** Purpose : Choose the vertical mixing scheme 118 !! 119 !! ** Method : Set nzdf from ln_zdfexp 120 !! nzdf = 0 explicit (time-splitting) scheme (ln_zdfexp=T) 121 !! = 1 implicit (euler backward) scheme (ln_zdfexp=F) 122 !! NB: rotation of lateral mixing operator or TKE & GLS schemes, 123 !! an implicit scheme is required. 124 !!---------------------------------------------------------------------- 125 USE zdftke 126 USE zdfgls 127 !!---------------------------------------------------------------------- 128 ! 129 ! Choice from ln_zdfexp already read in namelist in zdfini module 130 IF( ln_zdfexp ) THEN ; nzdf = 0 ! use explicit scheme 131 ELSE ; nzdf = 1 ! use implicit scheme 132 ENDIF 133 ! 134 ! Force implicit schemes 135 IF( lk_zdftke .OR. lk_zdfgls ) nzdf = 1 ! TKE, or GLS physics 136 IF( ln_traldf_iso ) nzdf = 1 ! iso-neutral lateral physics 137 IF( ln_traldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 138 IF( ln_zdfexp .AND. nzdf == 1 ) CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator', & 139 & ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 140 ! 141 IF(lwp) THEN 142 WRITE(numout,*) 143 WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme' 144 WRITE(numout,*) '~~~~~~~~~~~' 145 IF( nzdf == 0 ) WRITE(numout,*) ' ===>> Explicit time-splitting scheme' 146 IF( nzdf == 1 ) WRITE(numout,*) ' ===>> Implicit (euler backward) scheme' 147 ENDIF 148 ! 149 END SUBROUTINE tra_zdf_init 106 !! *** ROUTINE tra_zdf_imp *** 107 !! 108 !! ** Purpose : Compute the after tracer through a implicit computation 109 !! of the vertical tracer diffusion (including the vertical component 110 !! of lateral mixing (only for 2nd order operator, for fourth order 111 !! it is already computed and add to the general trend in traldf) 112 !! 113 !! ** Method : The vertical diffusion of a tracer ,t , is given by: 114 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 115 !! It is computed using a backward time scheme (t=after field) 116 !! which provide directly the after tracer field. 117 !! If ln_zdfddm=T, use avs for salinity or for passive tracers 118 !! Surface and bottom boundary conditions: no diffusive flux on 119 !! both tracers (bottom, applied through the masked field avt). 120 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 121 !! 122 !! ** Action : - pta becomes the after tracer 123 !!--------------------------------------------------------------------- 124 INTEGER , INTENT(in ) :: kt ! ocean time-step index 125 INTEGER , INTENT(in ) :: kit000 ! first time step index 126 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 127 INTEGER , INTENT(in ) :: kjpt ! number of tracers 128 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 129 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 130 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field 131 ! 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 REAL(wp) :: zrhs ! local scalars 134 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt, zwd, zws 135 !!--------------------------------------------------------------------- 136 ! 137 IF( ln_timing ) CALL timing_start('tra_zdf_imp') 138 ! 139 IF( kt == kit000 ) THEN 140 IF(lwp)WRITE(numout,*) 141 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 142 IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 143 ENDIF 144 ! ! ============= ! 145 DO jn = 1, kjpt ! tracer loop ! 146 ! ! ============= ! 147 ! Matrix construction 148 ! -------------------- 149 ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer 150 ! 151 IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. & 152 & ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN 153 ! 154 ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 155 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt(:,:,2:jpk) 156 ELSE ; zwt(:,:,2:jpk) = avs(:,:,2:jpk) 157 ENDIF 158 zwt(:,:,1) = 0._wp 159 ! 160 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 161 IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator 162 DO jk = 2, jpkm1 163 DO jj = 2, jpjm1 164 DO ji = fs_2, fs_jpim1 ! vector opt. 165 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 166 END DO 167 END DO 168 END DO 169 ELSE ! standard or triad iso-neutral operator 170 DO jk = 2, jpkm1 171 DO jj = 2, jpjm1 172 DO ji = fs_2, fs_jpim1 ! vector opt. 173 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 174 END DO 175 END DO 176 END DO 177 ENDIF 178 ENDIF 179 ! 180 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 181 DO jk = 1, jpkm1 182 DO jj = 2, jpjm1 183 DO ji = fs_2, fs_jpim1 ! vector opt. 184 !!gm BUG I think, use e3w_a instead of e3w_n, not sure of that 185 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) 186 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 187 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 188 END DO 189 END DO 190 END DO 191 ! 192 !! Matrix inversion from the first level 193 !!---------------------------------------------------------------------- 194 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 195 ! 196 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 197 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 198 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 199 ! ( ... )( ... ) ( ... ) 200 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 201 ! 202 ! m is decomposed in the product of an upper and lower triangular matrix. 203 ! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi. 204 ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal 205 ! and "superior" (above diagonal) components of the tridiagonal system. 206 ! The solution will be in the 4d array pta. 207 ! The 3d array zwt is used as a work space array. 208 ! En route to the solution pta is used a to evaluate the rhs and then 209 ! used as a work space array: its value is modified. 210 ! 211 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 212 DO ji = fs_2, fs_jpim1 ! done one for all passive tracers (so included in the IF instruction) 213 zwt(ji,jj,1) = zwd(ji,jj,1) 214 END DO 215 END DO 216 DO jk = 2, jpkm1 217 DO jj = 2, jpjm1 218 DO ji = fs_2, fs_jpim1 219 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 220 END DO 221 END DO 222 END DO 223 ! 224 ENDIF 225 ! 226 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 227 DO ji = fs_2, fs_jpim1 228 pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn) 229 END DO 230 END DO 231 DO jk = 2, jpkm1 232 DO jj = 2, jpjm1 233 DO ji = fs_2, fs_jpim1 234 zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn) ! zrhs=right hand side 235 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 236 END DO 237 END DO 238 END DO 239 ! 240 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 241 DO ji = fs_2, fs_jpim1 242 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 243 END DO 244 END DO 245 DO jk = jpk-2, 1, -1 246 DO jj = 2, jpjm1 247 DO ji = fs_2, fs_jpim1 248 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & 249 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 250 END DO 251 END DO 252 END DO 253 ! ! ================= ! 254 END DO ! end tracer loop ! 255 ! ! ================= ! 256 ! 257 IF( ln_timing ) CALL timing_stop('tra_zdf_imp') 258 ! 259 END SUBROUTINE tra_zdf_imp 150 260 151 261 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.