- Timestamp:
- 2020-03-20T23:26:56+01:00 (4 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90
r12481 r12581 1 MODULE traatf 1 MODULE traatfLF 2 2 !!====================================================================== 3 !! *** MODULE traatf ***3 !! *** MODULE traatfLF *** 4 4 !! Ocean active tracers: Asselin time filtering for temperature and salinity 5 5 !!====================================================================== … … 17 17 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter + modified LF-RA 18 18 !! - ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 19 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename tranxt.F90 -> traatf .F90. Now only does time filtering.19 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename tranxt.F90 -> traatfLF.F90. Now only does time filtering. 20 20 !!---------------------------------------------------------------------- 21 21 … … 26 26 !!---------------------------------------------------------------------- 27 27 USE oce ! ocean dynamics and tracers variables 28 USE dom_oce ! ocean space and time domain variables 28 USE dom_oce ! ocean space and time domain variables 29 29 USE sbc_oce ! surface boundary condition: ocean 30 30 USE sbcrnf ! river runoffs … … 33 33 USE domvvl ! variable volume 34 34 USE trd_oce ! trends: ocean variables 35 USE trdtra ! trends manager: tracers 35 USE trdtra ! trends manager: tracers 36 36 USE traqsr ! penetrative solar radiation (needed for nksr) 37 37 USE phycst ! physical constant … … 52 52 PRIVATE 53 53 54 PUBLIC tra_atf ! routine called by step.F9055 PUBLIC tra_atf_fix ! to be used in trcnxt56 PUBLIC tra_atf_ vvl ! to be used in trcnxt54 PUBLIC tra_atf_lf ! routine called by step.F90 55 PUBLIC tra_atf_fix_lf ! to be used in trcnxt !!st WARNING discrepancy here interpol is used 56 PUBLIC tra_atf_qe_lf ! to be used in trcnxt !!st WARNING discrepancy here interpol is used 57 57 58 58 !! * Substitutions … … 65 65 CONTAINS 66 66 67 SUBROUTINE tra_atf ( kt, Kbb, Kmm, Kaa, pts )68 !!---------------------------------------------------------------------- 69 !! *** ROUTINE traatf ***70 !! 71 !! ** Purpose : Apply the boundary condition on the after temperature 67 SUBROUTINE tra_atf_lf( kt, Kbb, Kmm, Kaa, pts ) 68 !!---------------------------------------------------------------------- 69 !! *** ROUTINE traatfLF *** 70 !! 71 !! ** Purpose : Apply the boundary condition on the after temperature 72 72 !! and salinity fields and add the Asselin time filter on now fields. 73 !! 74 !! ** Method : At this stage of the computation, ta and sa are the 73 !! 74 !! ** Method : At this stage of the computation, ta and sa are the 75 75 !! after temperature and salinity as the time stepping has 76 76 !! been performed in trazdf_imp or trazdf_exp module. 77 77 !! 78 !! - Apply lateral boundary conditions on (ta,sa) 79 !! at the local domain boundaries through lbc_lnk call, 80 !! at the one-way open boundaries (ln_bdy=T), 78 !! - Apply lateral boundary conditions on (ta,sa) 79 !! at the local domain boundaries through lbc_lnk call, 80 !! at the one-way open boundaries (ln_bdy=T), 81 81 !! at the AGRIF zoom boundaries (lk_agrif=T) 82 82 !! … … 88 88 INTEGER , INTENT(in ) :: kt ! ocean time-step index 89 89 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices 90 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 90 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 91 91 !! 92 92 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 95 95 !!---------------------------------------------------------------------- 96 96 ! 97 IF( ln_timing ) CALL timing_start( 'tra_atf ')97 IF( ln_timing ) CALL timing_start( 'tra_atf_lf') 98 98 ! 99 99 IF( kt == nit000 ) THEN 100 100 IF(lwp) WRITE(numout,*) 101 IF(lwp) WRITE(numout,*) 'tra_atf : apply Asselin time filter to "now" fields'101 IF(lwp) WRITE(numout,*) 'tra_atf_lf : apply Asselin time filter to "now" fields' 102 102 IF(lwp) WRITE(numout,*) '~~~~~~~' 103 103 ENDIF 104 104 105 105 ! Update after tracer on domain lateral boundaries 106 ! 106 ! 107 107 #if defined key_agrif 108 108 CALL Agrif_tra ! AGRIF zoom boundaries 109 109 #endif 110 110 ! ! local domain boundaries (T-point, unchanged sign) 111 CALL lbc_lnk_multi( 'traatf ', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )111 CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 112 112 ! 113 113 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 114 114 115 115 ! set time step size (Euler/Leapfrog) 116 116 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler) … … 119 119 120 120 ! trends computation initialisation 121 IF( l_trdtra ) THEN 121 IF( l_trdtra ) THEN 122 122 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 123 123 ztrdt(:,:,jpk) = 0._wp 124 124 ztrds(:,:,jpk) = 0._wp 125 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 125 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 126 126 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 127 127 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds ) 128 128 ENDIF 129 ! total trend for the non-time-filtered variables. 129 ! total trend for the non-time-filtered variables. 130 130 zfact = 1.0 / rdt 131 131 ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms … … 137 137 CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds ) 138 138 IF( ln_linssh ) THEN ! linear sea surface height only 139 ! Store now fields before applying the Asselin filter 139 ! Store now fields before applying the Asselin filter 140 140 ! in order to calculate Asselin filter trend later. 141 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 141 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 142 142 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm) 143 143 ENDIF 144 144 ENDIF 145 145 146 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping 146 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping 147 147 ! 148 148 IF (l_trdtra .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl … … 156 156 ELSE ! Leap-Frog + Asselin filter time stepping 157 157 ! 158 IF( ln_linssh ) THEN ; CALL tra_atf_fix ( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface159 ELSE ; CALL tra_atf_ vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface160 ENDIF 161 ! 162 CALL lbc_lnk_multi( 'traatf ', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., &158 IF( ln_linssh ) THEN ; CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 159 ELSE ; CALL tra_atf_qe_lf( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 160 ENDIF 161 ! 162 CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., & 163 163 & pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., & 164 164 & pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 165 165 ! 166 ENDIF 167 ! 168 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 169 zfact = 1._wp / r2dt 166 ENDIF 167 ! 168 IF( l_trdtra .AND. ln_linssh ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 169 zfact = 1._wp / r2dt 170 170 DO jk = 1, jpkm1 171 171 ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact … … 181 181 & tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2= ' Sn: ', mask2=tmask ) 182 182 ! 183 IF( ln_timing ) CALL timing_stop('tra_atf ')184 ! 185 END SUBROUTINE tra_atf 186 187 188 SUBROUTINE tra_atf_fix ( kt, Kbb, Kmm, Kaa, kit000, cdtype, pt, kjpt )183 IF( ln_timing ) CALL timing_stop('tra_atf_lf') 184 ! 185 END SUBROUTINE tra_atf_lf 186 187 188 SUBROUTINE tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, kit000, cdtype, pt, kjpt ) 189 189 !!---------------------------------------------------------------------- 190 190 !! *** ROUTINE tra_atf_fix *** 191 191 !! 192 192 !! ** Purpose : fixed volume: apply the Asselin time filter to the "now" field 193 !! 193 !! 194 194 !! ** Method : - Apply a Asselin time filter on now fields. 195 195 !! … … 209 209 IF( kt == kit000 ) THEN 210 210 IF(lwp) WRITE(numout,*) 211 IF(lwp) WRITE(numout,*) 'tra_atf_fix : time filtering', cdtype211 IF(lwp) WRITE(numout,*) 'tra_atf_fix_lf : time filtering', cdtype 212 212 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 213 213 ENDIF … … 216 216 ! 217 217 DO_3D_00_00( 1, jpkm1 ) 218 ztn = pt(ji,jj,jk,jn,Kmm) 218 ztn = pt(ji,jj,jk,jn,Kmm) 219 219 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers 220 220 ! … … 224 224 END DO 225 225 ! 226 END SUBROUTINE tra_atf_fix 227 228 229 SUBROUTINE tra_atf_ vvl( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt )226 END SUBROUTINE tra_atf_fix_lf 227 228 229 SUBROUTINE tra_atf_qe_lf( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt ) 230 230 !!---------------------------------------------------------------------- 231 231 !! *** ROUTINE tra_atf_vvl *** 232 232 !! 233 !! ** Purpose : Time varying volume: apply the Asselin time filter 234 !! 233 !! ** Purpose : Time varying volume: apply the Asselin time filter 234 !! 235 235 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 236 236 !! pt(Kmm) = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] ) … … 252 252 INTEGER :: ji, jj, jk, jn ! dummy loop indices 253 253 REAL(wp) :: zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar 254 REAL(wp) :: zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d , zscale! - -254 REAL(wp) :: zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d ! - - 255 255 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrd_atf 256 256 !!---------------------------------------------------------------------- … … 262 262 ENDIF 263 263 ! 264 IF( cdtype == 'TRA' ) THEN 264 IF( cdtype == 'TRA' ) THEN 265 265 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 266 266 ll_rnf = ln_rnf ! active tracers case and river runoffs … … 268 268 ELSE ! passive tracers case 269 269 ll_traqsr = .FALSE. ! NO solar penetration 270 ll_rnf = .FALSE. ! NO river runoffs ???? !!gm BUG ? 271 ll_isf = .FALSE. ! NO ice shelf melting/freezing !!gm BUG ?? 270 ll_rnf = .FALSE. ! NO river runoffs ???? !!gm BUG ? 271 ll_isf = .FALSE. ! NO ice shelf melting/freezing !!gm BUG ?? 272 272 ENDIF 273 273 ! … … 279 279 zfact1 = atfp * p2dt 280 280 zfact2 = zfact1 * r1_rau0 281 DO jn = 1, kjpt 281 DO jn = 1, kjpt 282 282 DO_3D_00_00( 1, jpkm1 ) 283 283 ze3t_b = e3t(ji,jj,jk,Kbb) … … 292 292 ztc_d = ztc_a - 2. * ztc_n + ztc_b 293 293 ! 294 ze3t_f = ze3t_n + atfp * ze3t_d295 294 ztc_f = ztc_n + atfp * ztc_d 296 295 ! 297 ! Add asselin correction on scale factors: 298 zscale = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 299 ze3t_f = ze3t_f - zfact2 * zscale * ( emp_b(ji,jj) - emp(ji,jj) ) 300 IF ( ll_rnf ) ze3t_f = ze3t_f + zfact2 * zscale * ( rnf_b(ji,jj) - rnf(ji,jj) ) 301 IF ( ll_isf ) THEN 302 IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) 303 IF ( ln_isfpar_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) ) 304 ENDIF 305 ! 306 IF( jk == mikt(ji,jj) ) THEN ! first level 296 ! Asselin correction on scale factors is done via ssh in r3t_f 297 ze3t_f = e3t_0(ji,jj,jk) * ( 1._wp + r3t_f(ji,jj) * tmask(ji,jj,jk) ) 298 299 ! 300 IF( jk == mikt(ji,jj) ) THEN ! first level 307 301 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 308 302 ENDIF 309 303 ! 310 304 ! solar penetration (temperature only) 311 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & 312 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 305 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & 306 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 313 307 ! 314 308 ! 315 309 IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & 316 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 310 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 317 311 & * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) 318 312 … … 328 322 & * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) 329 323 END IF 330 ! level partially include in Losch_2008 ice shelf boundary layer 324 ! level partially include in Losch_2008 ice shelf boundary layer 331 325 IF ( jk == misfkb_cav(ji,jj) ) THEN 332 326 ztc_f = ztc_f - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) & … … 342 336 & * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) 343 337 END IF 344 ! level partially include in Losch_2008 ice shelf boundary layer 338 ! level partially include in Losch_2008 ice shelf boundary layer 345 339 IF ( jk == misfkb_par(ji,jj) ) THEN 346 340 ztc_f = ztc_f - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) ) & … … 356 350 ztc_f = ztc_f + zfact1 * risfcpl_tsc(ji,jj,jk,jn) * r1_e1e2t(ji,jj) 357 351 ! Shouldn't volume increment be spread according thanks to zscale ? 358 ze3t_f = ze3t_f - zfact1 * risfcpl_vol(ji,jj,jk ) * r1_e1e2t(ji,jj)359 352 END IF 360 353 ! … … 371 364 ! 372 365 END_3D 373 ! 366 ! 374 367 END DO 375 368 ! 376 369 IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) ) THEN 377 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 370 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 378 371 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 379 372 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) … … 387 380 ENDIF 388 381 ! 389 END SUBROUTINE tra_atf_ vvl382 END SUBROUTINE tra_atf_qe_lf 390 383 391 384 !!====================================================================== 392 END MODULE traatf 385 END MODULE traatfLF
Note: See TracChangeset
for help on using the changeset viewer.
