- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/trazdf.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/trazdf.F90
r10425 r13463 36 36 37 37 !! * Substitutions 38 # include "vectopt_loop_substitute.h90" 38 # include "do_loop_substitute.h90" 39 # include "domzgr_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 41 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 44 45 CONTAINS 45 46 46 SUBROUTINE tra_zdf( kt )47 SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa ) 47 48 !!---------------------------------------------------------------------- 48 49 !! *** ROUTINE tra_zdf *** … … 50 51 !! ** Purpose : compute the vertical ocean tracer physics. 51 52 !!--------------------------------------------------------------------- 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 53 INTEGER , INTENT(in) :: kt ! ocean time-step index 54 INTEGER , INTENT(in) :: Kbb, Kmm, Krhs, Kaa ! time level indices 55 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 53 56 ! 54 57 INTEGER :: jk ! Dummy loop indices … … 64 67 ENDIF 65 68 ! 66 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000, = rdt (restarting with Euler time stepping)67 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! otherwise, = 2 rdt (leapfrog)68 ENDIF69 !70 69 IF( l_trdtra ) THEN !* Save ta and sa trends 71 70 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)71 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 72 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 74 73 ENDIF 75 74 ! 76 75 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r 2dt, tsb, tsa, jpts )76 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 78 77 79 78 !!gm WHY here ! and I don't like that ! … … 81 80 ! JMM avoid negative salinities near river outlet ! Ugly fix 82 81 ! JMM : restore negative salinities to small salinities: 83 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp82 WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp 84 83 !!gm 85 84 86 85 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 86 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) 87 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) & 88 & - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 89 & / ( e3t(:,:,jk,Kmm)*rDt ) ) & 90 & - ztrdt(:,:,jk) 91 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) & 92 & - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / ( e3t(:,:,jk,Kmm)*rDt ) ) & 94 & - ztrds(:,:,jk) 92 95 END DO 93 96 !!gm this should be moved in trdtra.F90 and done on all trends 94 CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1.)97 CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1.0_wp , ztrds, 'T', 1.0_wp ) 95 98 !!gm 96 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )97 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )99 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 100 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 98 101 DEALLOCATE( ztrdt , ztrds ) 99 102 ENDIF 100 103 ! ! 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' )104 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, & 105 & tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 103 106 ! 104 107 IF( ln_timing ) CALL timing_stop('tra_zdf') … … 107 110 108 111 109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )112 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 110 113 !!---------------------------------------------------------------------- 111 114 !! *** ROUTINE tra_zdf_imp *** … … 125 128 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 126 129 !! 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 field130 !! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer 131 !!--------------------------------------------------------------------- 132 INTEGER , INTENT(in ) :: kt ! ocean time-step index 133 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 134 INTEGER , INTENT(in ) :: kit000 ! first time step index 135 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 136 INTEGER , INTENT(in ) :: kjpt ! number of tracers 137 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 138 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 136 139 ! 137 140 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 158 161 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 159 162 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 163 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 164 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 165 END_3D 167 166 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 167 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 168 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 169 END_3D 175 170 ENDIF 176 171 ENDIF … … 178 173 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 179 174 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 175 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 176 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 177 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 178 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 179 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 180 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) 181 zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) 182 END_3D 192 183 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 184 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 185 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 186 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 187 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 188 END_3D 202 189 ENDIF 203 190 ! … … 218 205 ! The solution will be in the 4d array pta. 219 206 ! 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 then207 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 221 208 ! used as a work space array: its value is modified. 222 209 ! 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 210 DO_2D( 0, 0, 0, 0 ) 211 zwt(ji,jj,1) = zwd(ji,jj,1) 212 END_2D 213 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 214 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 215 END_3D 235 216 ! 236 217 ENDIF 237 218 ! 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 219 DO_2D( 0, 0, 0, 0 ) 220 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) & 221 & + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 222 END_2D 223 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 224 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) & 225 & + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side 226 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 227 END_3D 251 228 ! 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 229 DO_2D( 0, 0, 0, 0 ) 230 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 231 END_2D 232 DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 233 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 234 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 235 END_3D 265 236 ! ! ================= ! 266 237 END DO ! end tracer loop !
Note: See TracChangeset
for help on using the changeset viewer.
