Changeset 11057 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90
- Timestamp:
- 2019-05-24T17:36:26+02:00 (5 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90
r11056 r11057 1 MODULE tra nxt1 MODULE traatf 2 2 !!====================================================================== 3 !! *** MODULE tra nxt***3 !! *** MODULE traatf *** 4 4 !! Ocean active tracers: time stepping on temperature and salinity 5 5 !!====================================================================== … … 20 20 21 21 !!---------------------------------------------------------------------- 22 !! tra_ nxt: time stepping on tracers23 !! tra_ nxt_fix : time stepping on tracers : fixed volume case24 !! tra_ nxt_vvl : time stepping on tracers : variable volume case22 !! tra_atf : time stepping on tracers 23 !! tra_atf_fix : time stepping on tracers : fixed volume case 24 !! tra_atf_vvl : time stepping on tracers : variable volume case 25 25 !!---------------------------------------------------------------------- 26 26 USE oce ! ocean dynamics and tracers variables … … 51 51 PRIVATE 52 52 53 PUBLIC tra_ nxt! routine called by step.F9054 PUBLIC tra_ nxt_fix ! to be used in trcnxt55 PUBLIC tra_ nxt_vvl ! to be used in trcnxt53 PUBLIC tra_atf ! routine called by step.F90 54 PUBLIC tra_atf_fix ! to be used in trcnxt 55 PUBLIC tra_atf_vvl ! to be used in trcnxt 56 56 57 57 !! * Substitutions … … 64 64 CONTAINS 65 65 66 SUBROUTINE tra_nxt( kt, Kbb, Kmm, Krhs, Kaa ) 67 !!---------------------------------------------------------------------- 68 !! *** ROUTINE tranxt *** 66 SUBROUTINE tra_atf( kt, Kbb, Kmm, Kaa, pts ) 67 !!---------------------------------------------------------------------- 68 !! *** ROUTINE traatf *** 69 !! 70 !! !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!! 69 71 !! 70 72 !! ** Purpose : Apply the boundary condition on the after temperature … … 86 88 !! ** Action : - ts(Kbb) & ts(Kmm) ready for the next time step 87 89 !!---------------------------------------------------------------------- 88 INTEGER, INTENT(in) :: kt ! ocean time-step index 89 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs, Kaa ! time level indices 90 INTEGER , INTENT(in ) :: kt ! ocean time-step index 91 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices 92 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 90 93 !! 91 94 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 94 97 !!---------------------------------------------------------------------- 95 98 ! 96 IF( ln_timing ) CALL timing_start( 'tra_ nxt')99 IF( ln_timing ) CALL timing_start( 'tra_atf') 97 100 ! 98 101 IF( kt == nit000 ) THEN 99 102 IF(lwp) WRITE(numout,*) 100 IF(lwp) WRITE(numout,*) 'tra_ nxt: achieve the time stepping by Asselin filter and array swap'103 IF(lwp) WRITE(numout,*) 'tra_atf : achieve the time stepping by Asselin filter and array swap' 101 104 IF(lwp) WRITE(numout,*) '~~~~~~~' 102 105 ENDIF … … 108 111 #endif 109 112 ! ! local domain boundaries (T-point, unchanged sign) 110 CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. ) 111 ! 112 !! IMMERSE development : Following the general pattern for the new code we want to pass in the 113 !! velocities to bdy_dyn as arguments so here we use "ts" even 114 !! though we haven't converted the tracer names in the rest of tranxt.F90 115 !! because it will be completely rewritten. DS. 116 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, ts, Kaa ) ! BDY open boundaries 113 CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 114 ! 115 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 117 116 118 117 ! set time step size (Euler/Leapfrog) … … 127 126 ztrds(:,:,jpk) = 0._wp 128 127 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 129 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_tem, jptra_zdfp, ztrdt )130 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_sal, jptra_zdfp, ztrds )128 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 129 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds ) 131 130 ENDIF 132 131 ! total trend for the non-time-filtered variables. 133 132 zfact = 1.0 / rdt 134 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms133 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 135 134 DO jk = 1, jpkm1 136 ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) -ts(:,:,jk,jp_tem,Kmm)) * zfact137 ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Krhs)*e3t(:,:,jk,Krhs) / e3t(:,:,jk,Kmm) -ts(:,:,jk,jp_sal,Kmm)) * zfact138 END DO 139 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_tem, jptra_tot, ztrdt )140 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_sal, jptra_tot, ztrds )135 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_tem,Kmm)) * zfact 136 ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_sal,Kmm)) * zfact 137 END DO 138 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_tot, ztrdt ) 139 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds ) 141 140 IF( ln_linssh ) THEN ! linear sea surface height only 142 141 ! Store now fields before applying the Asselin filter 143 142 ! in order to calculate Asselin filter trend later. 144 ztrdt(:,:,:) = ts(:,:,:,jp_tem,Kmm)145 ztrds(:,:,:) = ts(:,:,:,jp_sal,Kmm)143 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 144 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm) 146 145 ENDIF 147 146 ENDIF … … 150 149 DO jn = 1, jpts 151 150 DO jk = 1, jpkm1 152 ts(:,:,jk,jn,Kmm) = ts(:,:,jk,jn,Krhs)151 pts(:,:,jk,jn,Kmm) = pts(:,:,jk,jn,Kaa) 153 152 END DO 154 153 END DO 155 154 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl 156 ! ! Asselin filter is output by tra_ nxt_vvl that is not called on this time step155 ! ! Asselin filter is output by tra_atf_vvl that is not called on this time step 157 156 ztrdt(:,:,:) = 0._wp 158 157 ztrds(:,:,:) = 0._wp 159 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_tem, jptra_atf, ztrdt )160 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_sal, jptra_atf, ztrds )158 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt ) 159 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds ) 161 160 END IF 162 161 ! 163 162 ELSE ! Leap-Frog + Asselin filter time stepping 164 163 ! 165 IF( ln_linssh ) THEN ; CALL tra_nxt_fix( kt, Kmm, nit000, 'TRA', & 166 & ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs), jpts ) ! linear free surface 167 ELSE ; CALL tra_nxt_vvl( kt, Kbb, Kmm, Krhs, nit000, rdt, 'TRA', & 168 & ts(:,:,:,:,Kbb), ts(:,:,:,:,Kmm), ts(:,:,:,:,Krhs), & 169 & sbc_tsc , sbc_tsc_b , jpts ) ! non-linear free surface 170 ENDIF 171 ! 172 CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Kbb) , 'T', 1., ts(:,:,:,jp_sal,Kbb) , 'T', 1., & 173 & ts(:,:,:,jp_tem,Kmm) , 'T', 1., ts(:,:,:,jp_sal,Kmm) , 'T', 1., & 174 & ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. ) 164 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 165 ELSE ; CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 166 ENDIF 167 ! 168 CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., & 169 & pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., & 170 & pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 175 171 ! 176 172 ENDIF … … 179 175 zfact = 1._wp / r2dt 180 176 DO jk = 1, jpkm1 181 ztrdt(:,:,jk) = ( ts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact182 ztrds(:,:,jk) = ( ts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact183 END DO 184 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_tem, jptra_atf, ztrdt )185 CALL trd_tra( kt, Kmm, K rhs, 'TRA', jp_sal, jptra_atf, ztrds )177 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact 178 ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact 179 END DO 180 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt ) 181 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds ) 186 182 END IF 187 183 IF( l_trdtra ) DEALLOCATE( ztrdt , ztrds ) 188 184 ! 189 185 ! ! control print 190 IF(ln_ctl) CALL prt_ctl( tab3d_1= ts(:,:,:,jp_tem,Kmm), clinfo1=' nxt - Tn: ', mask1=tmask, &191 & tab3d_2= ts(:,:,:,jp_sal,Kmm), clinfo2= ' Sn: ', mask2=tmask )192 ! 193 IF( ln_timing ) CALL timing_stop('tra_ nxt')194 ! 195 END SUBROUTINE tra_ nxt196 197 198 SUBROUTINE tra_ nxt_fix( kt, Kmm, kit000, cdtype, ptb, ptn, pta, kjpt )199 !!---------------------------------------------------------------------- 200 !! *** ROUTINE tra_ nxt_fix ***186 IF(ln_ctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kmm), clinfo1=' nxt - Tn: ', mask1=tmask, & 187 & tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2= ' Sn: ', mask2=tmask ) 188 ! 189 IF( ln_timing ) CALL timing_stop('tra_atf') 190 ! 191 END SUBROUTINE tra_atf 192 193 194 SUBROUTINE tra_atf_fix( kt, kit000, Kbb, Kmm, Kaa, cdtype, pt, kjpt ) 195 !!---------------------------------------------------------------------- 196 !! *** ROUTINE tra_atf_fix *** 201 197 !! 202 198 !! ** Purpose : fixed volume: apply the Asselin time filter and … … 208 204 !! ** Action : - ptb & ptn ready for the next time step 209 205 !!---------------------------------------------------------------------- 210 INTEGER , INTENT(in ) :: kt ! ocean time-step index 211 INTEGER , INTENT(in ) :: Kmm ! time level index 212 INTEGER , INTENT(in ) :: kit000 ! first time step index 213 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 214 INTEGER , INTENT(in ) :: kjpt ! number of tracers 215 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb ! before tracer fields 216 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptn ! now tracer fields 217 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 206 INTEGER , INTENT(in ) :: kt ! ocean time-step index 207 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices 208 INTEGER , INTENT(in ) :: kit000 ! first time step index 209 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 210 INTEGER , INTENT(in ) :: kjpt ! number of tracers 211 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracer fields 218 212 ! 219 213 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 223 217 IF( kt == kit000 ) THEN 224 218 IF(lwp) WRITE(numout,*) 225 IF(lwp) WRITE(numout,*) 'tra_ nxt_fix : time stepping', cdtype219 IF(lwp) WRITE(numout,*) 'tra_atf_fix : time stepping', cdtype 226 220 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 227 221 ENDIF … … 232 226 DO jj = 2, jpjm1 233 227 DO ji = fs_2, fs_jpim1 234 ztn = ptn(ji,jj,jk,jn) 235 ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 236 ! 237 ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn 238 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 228 ztn = pt(ji,jj,jk,jn,Kmm) 229 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers 230 ! 231 pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd ! ptb <-- filtered ptn 239 232 END DO 240 233 END DO … … 243 236 END DO 244 237 ! 245 END SUBROUTINE tra_ nxt_fix246 247 248 SUBROUTINE tra_ nxt_vvl( kt, Kbb, Kmm, Krhs, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt )249 !!---------------------------------------------------------------------- 250 !! *** ROUTINE tra_ nxt_vvl ***238 END SUBROUTINE tra_atf_fix 239 240 241 SUBROUTINE tra_atf_vvl( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt ) 242 !!---------------------------------------------------------------------- 243 !! *** ROUTINE tra_atf_vvl *** 251 244 !! 252 245 !! ** Purpose : Time varying volume: apply the Asselin time filter … … 256 249 !! - swap tracer fields to prepare the next time_step. 257 250 !! tb = ( e3t(Kmm)*tn + atfp*[ e3t(Kbb)*tb - 2 e3t(Kmm)*tn + e3t_a*ta ] ) 258 !! /( e3t(Kmm) + atfp*[ e3t(Kbb) - 2 e3t(Kmm) + e3t(K rhs) ] )251 !! /( e3t(Kmm) + atfp*[ e3t(Kbb) - 2 e3t(Kmm) + e3t(Kaa) ] ) 259 252 !! tn = ta 260 253 !! 261 254 !! ** Action : - ptb & ptn ready for the next time step 262 255 !!---------------------------------------------------------------------- 263 INTEGER , INTENT(in ) :: kt ! ocean time-step index 264 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! time level indices 265 INTEGER , INTENT(in ) :: kit000 ! first time step index 266 REAL(wp) , INTENT(in ) :: p2dt ! time-step 267 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 268 INTEGER , INTENT(in ) :: kjpt ! number of tracers 269 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb ! before tracer fields 270 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptn ! now tracer fields 271 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 272 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: psbc_tc ! surface tracer content 273 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: psbc_tc_b ! before surface tracer content 256 INTEGER , INTENT(in ) :: kt ! ocean time-step index 257 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices 258 INTEGER , INTENT(in ) :: kit000 ! first time step index 259 REAL(wp) , INTENT(in ) :: p2dt ! time-step 260 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 261 INTEGER , INTENT(in ) :: kjpt ! number of tracers 262 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracer fields 263 REAL(wp), DIMENSION(jpi,jpj ,kjpt) , INTENT(in ) :: psbc_tc ! surface tracer content 264 REAL(wp), DIMENSION(jpi,jpj ,kjpt) , INTENT(in ) :: psbc_tc_b ! before surface tracer content 274 265 ! 275 266 LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical … … 282 273 IF( kt == kit000 ) THEN 283 274 IF(lwp) WRITE(numout,*) 284 IF(lwp) WRITE(numout,*) 'tra_ nxt_vvl : time stepping', cdtype275 IF(lwp) WRITE(numout,*) 'tra_atf_vvl : time stepping', cdtype 285 276 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 286 277 ENDIF … … 309 300 ze3t_b = e3t(ji,jj,jk,Kbb) 310 301 ze3t_n = e3t(ji,jj,jk,Kmm) 311 ze3t_a = e3t(ji,jj,jk,K rhs)302 ze3t_a = e3t(ji,jj,jk,Kaa) 312 303 ! ! tracer content at Before, now and after 313 ztc_b = pt b(ji,jj,jk,jn) * ze3t_b314 ztc_n = pt n(ji,jj,jk,jn) * ze3t_n315 ztc_a = pt a(ji,jj,jk,jn) * ze3t_a304 ztc_b = pt(ji,jj,jk,jn,Kbb) * ze3t_b 305 ztc_n = pt(ji,jj,jk,jn,Kmm) * ze3t_n 306 ztc_a = pt(ji,jj,jk,jn,Kaa) * ze3t_a 316 307 ! 317 308 ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b … … 361 352 ! 362 353 ze3t_f = 1.e0 / ze3t_f 363 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 364 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 354 pt(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f ! time filtered "now" field 365 355 ! 366 356 IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN … … 376 366 IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) ) THEN 377 367 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 378 CALL trd_tra( kt, Kmm, K rhs, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )379 CALL trd_tra( kt, Kmm, K rhs, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )368 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 369 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) 380 370 ENDIF 381 371 IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN 382 372 DO jn = 1, kjpt 383 CALL trd_tra( kt, Kmm, K rhs, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )373 CALL trd_tra( kt, Kmm, Kaa, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) ) 384 374 END DO 385 375 ENDIF … … 387 377 ENDIF 388 378 ! 389 END SUBROUTINE tra_ nxt_vvl379 END SUBROUTINE tra_atf_vvl 390 380 391 381 !!====================================================================== 392 END MODULE tra nxt382 END MODULE traatf
Note: See TracChangeset
for help on using the changeset viewer.
