Changeset 12377 for NEMO/trunk/src/OCE/TRA/trazdf.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/TRA/trazdf.F90
r10425 r12377 36 36 37 37 !! * Substitutions 38 # include " vectopt_loop_substitute.h90"38 # include "do_loop_substitute.h90" 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 44 44 CONTAINS 45 45 46 SUBROUTINE tra_zdf( kt )46 SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE tra_zdf *** … … 50 50 !! ** Purpose : compute the vertical ocean tracer physics. 51 51 !!--------------------------------------------------------------------- 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 INTEGER , INTENT(in) :: kt ! ocean time-step index 53 INTEGER , INTENT(in) :: Kbb, Kmm, Krhs, Kaa ! time level indices 54 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 53 55 ! 54 56 INTEGER :: jk ! Dummy loop indices … … 70 72 IF( l_trdtra ) THEN !* Save ta and sa trends 71 73 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)74 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 75 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 74 76 ENDIF 75 77 ! 76 78 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )79 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 78 80 79 81 !!gm WHY here ! and I don't like that ! … … 81 83 ! JMM avoid negative salinities near river outlet ! Ugly fix 82 84 ! JMM : restore negative salinities to small salinities: 83 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp85 WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp 84 86 !!gm 85 87 86 88 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 89 DO jk = 1, jpkm1 88 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &89 & / (e3t _n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)90 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &91 & / (e3t _n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)90 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 91 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 92 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 92 94 END DO 93 95 !!gm this should be moved in trdtra.F90 and done on all trends 94 96 CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. ) 95 97 !!gm 96 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )97 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )98 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 99 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 98 100 DEALLOCATE( ztrdt , ztrds ) 99 101 ENDIF 100 102 ! ! print mean trends (used for debugging) 101 IF( ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf - Ta: ', mask1=tmask, &102 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )103 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, & 104 & tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 103 105 ! 104 106 IF( ln_timing ) CALL timing_stop('tra_zdf') … … 107 109 108 110 109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )111 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 110 112 !!---------------------------------------------------------------------- 111 113 !! *** ROUTINE tra_zdf_imp *** … … 125 127 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 126 128 !! 127 !! ** Action : - pt abecomes the after tracer128 !!--------------------------------------------------------------------- 129 INTEGER , INTENT(in ) :: kt ! ocean time-step index130 INTEGER , INTENT(in ) :: kit000 ! first time step index131 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)132 INTEGER , INTENT(in ) :: kjpt ! number of tracers133 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step134 REAL(wp) , DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields135 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt ), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field129 !! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer 130 !!--------------------------------------------------------------------- 131 INTEGER , INTENT(in ) :: kt ! ocean time-step index 132 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 133 INTEGER , INTENT(in ) :: kit000 ! first time step index 134 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 135 INTEGER , INTENT(in ) :: kjpt ! number of tracers 136 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 137 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 136 138 ! 137 139 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 158 160 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 159 161 IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator 160 DO jk = 2, jpkm1 161 DO jj = 2, jpjm1 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 164 END DO 165 END DO 166 END DO 162 DO_3D_00_00( 2, jpkm1 ) 163 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 164 END_3D 167 165 ELSE ! standard or triad iso-neutral operator 168 DO jk = 2, jpkm1 169 DO jj = 2, jpjm1 170 DO ji = fs_2, fs_jpim1 ! vector opt. 171 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 172 END DO 173 END DO 174 END DO 166 DO_3D_00_00( 2, jpkm1 ) 167 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 168 END_3D 175 169 ENDIF 176 170 ENDIF … … 178 172 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 179 173 IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection 180 DO jk = 1, jpkm1 181 DO jj = 2, jpjm1 182 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 183 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) 184 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 185 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws & 186 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 187 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) 188 zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) 189 END DO 190 END DO 191 END DO 174 DO_3D_00_00( 1, jpkm1 ) 175 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 176 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 177 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 178 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 179 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) 180 zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) 181 END_3D 192 182 ELSE 193 DO jk = 1, jpkm1 194 DO jj = 2, jpjm1 195 DO ji = fs_2, fs_jpim1 ! vector opt. 196 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk) 197 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 198 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 199 END DO 200 END DO 201 END DO 183 DO_3D_00_00( 1, jpkm1 ) 184 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 185 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 186 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 187 END_3D 202 188 ENDIF 203 189 ! … … 218 204 ! The solution will be in the 4d array pta. 219 205 ! The 3d array zwt is used as a work space array. 220 ! En route to the solution pt ais used a to evaluate the rhs and then206 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 221 207 ! used as a work space array: its value is modified. 222 208 ! 223 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 224 DO ji = fs_2, fs_jpim1 ! done one for all passive tracers (so included in the IF instruction) 225 zwt(ji,jj,1) = zwd(ji,jj,1) 226 END DO 227 END DO 228 DO jk = 2, jpkm1 229 DO jj = 2, jpjm1 230 DO ji = fs_2, fs_jpim1 231 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 232 END DO 233 END DO 234 END DO 209 DO_2D_00_00 210 zwt(ji,jj,1) = zwd(ji,jj,1) 211 END_2D 212 DO_3D_00_00( 2, jpkm1 ) 213 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 214 END_3D 235 215 ! 236 216 ENDIF 237 217 ! 238 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 239 DO ji = fs_2, fs_jpim1 240 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) 241 END DO 242 END DO 243 DO jk = 2, jpkm1 244 DO jj = 2, jpjm1 245 DO ji = fs_2, fs_jpim1 246 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 247 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 248 END DO 249 END DO 250 END DO 218 DO_2D_00_00 219 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 220 END_2D 221 DO_3D_00_00( 2, jpkm1 ) 222 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side 223 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 224 END_3D 251 225 ! 252 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 253 DO ji = fs_2, fs_jpim1 254 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 255 END DO 256 END DO 257 DO jk = jpk-2, 1, -1 258 DO jj = 2, jpjm1 259 DO ji = fs_2, fs_jpim1 260 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & 261 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 262 END DO 263 END DO 264 END DO 226 DO_2D_00_00 227 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 228 END_2D 229 DO_3DS_00_00( jpk-2, 1, -1 ) 230 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 231 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 232 END_3D 265 233 ! ! ================= ! 266 234 END DO ! end tracer loop !
Note: See TracChangeset
for help on using the changeset viewer.