- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/trazdf.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/trazdf.F90
r10425 r12928 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 … … 64 66 ENDIF 65 67 ! 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 68 IF( l_trdtra ) THEN !* Save ta and sa trends 71 69 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)70 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 71 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 74 72 ENDIF 75 73 ! 76 74 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r 2dt, tsb, tsa, jpts )75 CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 78 76 79 77 !!gm WHY here ! and I don't like that ! … … 81 79 ! JMM avoid negative salinities near river outlet ! Ugly fix 82 80 ! JMM : restore negative salinities to small salinities: 83 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp81 WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp 84 82 !!gm 85 83 86 84 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 85 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)86 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 87 & / (e3t(:,:,jk,Kmm)*rDt) ) - ztrdt(:,:,jk) 88 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 89 & / (e3t(:,:,jk,Kmm)*rDt) ) - ztrds(:,:,jk) 92 90 END DO 93 91 !!gm this should be moved in trdtra.F90 and done on all trends 94 92 CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. ) 95 93 !!gm 96 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )97 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )94 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 95 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 98 96 DEALLOCATE( ztrdt , ztrds ) 99 97 ENDIF 100 98 ! ! 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' )99 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, & 100 & tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 103 101 ! 104 102 IF( ln_timing ) CALL timing_stop('tra_zdf') … … 107 105 108 106 109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )107 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 110 108 !!---------------------------------------------------------------------- 111 109 !! *** ROUTINE tra_zdf_imp *** … … 125 123 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 126 124 !! 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 field125 !! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer 126 !!--------------------------------------------------------------------- 127 INTEGER , INTENT(in ) :: kt ! ocean time-step index 128 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 129 INTEGER , INTENT(in ) :: kit000 ! first time step index 130 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 131 INTEGER , INTENT(in ) :: kjpt ! number of tracers 132 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 133 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 136 134 ! 137 135 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 158 156 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 159 157 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 158 DO_3D_00_00( 2, jpkm1 ) 159 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 160 END_3D 167 161 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 162 DO_3D_00_00( 2, jpkm1 ) 163 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 164 END_3D 175 165 ENDIF 176 166 ENDIF … … 178 168 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 179 169 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 170 DO_3D_00_00( 1, jpkm1 ) 171 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 172 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 173 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 174 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 175 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) 176 zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) 177 END_3D 192 178 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 179 DO_3D_00_00( 1, jpkm1 ) 180 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 181 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 182 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 183 END_3D 202 184 ENDIF 203 185 ! … … 218 200 ! The solution will be in the 4d array pta. 219 201 ! 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 then202 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 221 203 ! used as a work space array: its value is modified. 222 204 ! 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 205 DO_2D_00_00 206 zwt(ji,jj,1) = zwd(ji,jj,1) 207 END_2D 208 DO_3D_00_00( 2, jpkm1 ) 209 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 210 END_3D 235 211 ! 236 212 ENDIF 237 213 ! 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 214 DO_2D_00_00 215 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) 216 END_2D 217 DO_3D_00_00( 2, jpkm1 ) 218 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 219 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 220 END_3D 251 221 ! 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 222 DO_2D_00_00 223 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 224 END_2D 225 DO_3DS_00_00( jpk-2, 1, -1 ) 226 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 227 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 228 END_3D 265 229 ! ! ================= ! 266 230 END DO ! end tracer loop !
Note: See TracChangeset
for help on using the changeset viewer.
