Changeset 12581
- Timestamp:
- 2020-03-20T23:26:56+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE
- Files:
-
- 4 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
r12579 r12581 45 45 PUBLIC dom_qe_sf_update ! called by step.F90 46 46 PUBLIC dom_qe_interpol ! called by dynnxt.F90 47 PUBLIC dom_qe_r3c ! called by steplf.F90 47 48 48 49 ! !!* Namelist nam_vvl -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatf.F90
r12377 r12581 13 13 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 14 14 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 16 16 !! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module … … 22 22 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering. 23 23 !!------------------------------------------------------------------------- 24 24 25 25 !!---------------------------------------------------------------------------------------------- 26 26 !! dyn_atf : apply Asselin time filtering to "now" velocities and vertical scale factors … … 42 42 USE trdken ! trend manager: kinetic energy 43 43 USE isf_oce , ONLY: ln_isf ! ice shelf 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 45 45 ! 46 46 USE in_out_manager ! I/O manager … … 63 63 !!---------------------------------------------------------------------- 64 64 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 65 !! $Id$ 65 !! $Id$ 66 66 !! Software governed by the CeCILL license (see ./LICENSE) 67 67 !!---------------------------------------------------------------------- … … 71 71 !!---------------------------------------------------------------------- 72 72 !! *** ROUTINE dyn_atf *** 73 !! 74 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 73 !! 74 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 75 75 !! condition on the after velocity and apply the Asselin time 76 76 !! filter to the now fields. … … 79 79 !! estimate (ln_dynspg_ts=T) 80 80 !! 81 !! * Apply lateral boundary conditions on after velocity 81 !! * Apply lateral boundary conditions on after velocity 82 82 !! at the local domain boundaries through lbc_lnk call, 83 83 !! at the one-way open boundaries (ln_bdy=T), … … 86 86 !! * Apply the Asselin time filter to the now fields 87 87 !! arrays to start the next time step: 88 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 88 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 89 89 !! + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 90 90 !! Note that with flux form advection and non linear free surface, … … 92 92 !! As a result, dyn_atf MUST be called after tra_atf. 93 93 !! 94 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 94 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 95 95 !!---------------------------------------------------------------------- 96 96 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 103 103 REAL(wp) :: zve3a, zve3n, zve3b, z1_2dt ! - - 104 104 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve, zwfld 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 106 106 !!---------------------------------------------------------------------- 107 107 ! … … 131 131 ! 132 132 IF( .NOT.ln_bt_fw ) THEN 133 ! Remove advective velocity from "now velocities" 134 ! prior to asselin filtering 135 ! In the forward case, this is done below after asselin filtering 136 ! so that asselin contribution is removed at the same time 133 ! Remove advective velocity from "now velocities" 134 ! prior to asselin filtering 135 ! In the forward case, this is done below after asselin filtering 136 ! so that asselin contribution is removed at the same time 137 137 DO jk = 1, jpkm1 138 138 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 139 139 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 140 END DO 140 END DO 141 141 ENDIF 142 142 ENDIF 143 143 144 144 ! Update after velocity on domain lateral boundaries 145 ! -------------------------------------------------- 145 ! -------------------------------------------------- 146 146 # if defined key_agrif 147 147 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries … … 157 157 ! 158 158 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 160 160 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt 161 161 ! … … 177 177 ! Time filter and swap of dynamics arrays 178 178 ! ------------------------------------------ 179 180 IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN !* Leap-Frog : Asselin time filter 179 180 IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN !* Leap-Frog : Asselin time filter 181 181 ! ! =============! 182 182 IF( ln_linssh ) THEN ! Fixed volume ! … … 202 202 DO jk = 1, jpkm1 203 203 ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 204 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 204 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 205 205 END DO 206 206 ! … … 239 239 pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 240 240 END_3D 241 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 241 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 242 242 pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 243 243 ! … … 250 250 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 251 251 ! Revert filtered "now" velocities to time split estimate 252 ! Doing it here also means that asselin filter contribution is removed 252 ! Doing it here also means that asselin filter contribution is removed 253 253 zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 254 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 254 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 255 255 DO jk = 2, jpkm1 256 256 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 257 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 257 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 258 258 END DO 259 259 DO jk = 1, jpkm1 … … 307 307 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, & 308 308 & tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask ) 309 ! 309 ! 310 310 IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) 311 311 IF( l_trddyn ) DEALLOCATE( zua, zva ) -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90
r12481 r12581 1 MODULE dynatf 1 MODULE dynatfLF 2 2 !!========================================================================= 3 !! *** MODULE dynatf ***3 !! *** MODULE dynatfLF *** 4 4 !! Ocean dynamics: time filtering 5 5 !!========================================================================= … … 13 13 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 14 14 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 16 16 !! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module … … 20 20 !! 3.6 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 22 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf .F90. Now just does time filtering.22 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatfLF.F90. Now just does time filtering. 23 23 !!------------------------------------------------------------------------- 24 24 25 25 !!---------------------------------------------------------------------------------------------- 26 !! dyn_atf : apply Asselin time filtering to "now" velocities and vertical scale factors26 !! dyn_atfLF : apply Asselin time filtering to "now" velocities and vertical scale factors 27 27 !!---------------------------------------------------------------------------------------------- 28 28 USE oce ! ocean dynamics and tracers … … 42 42 USE trdken ! trend manager: kinetic energy 43 43 USE isf_oce , ONLY: ln_isf ! ice shelf 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 45 45 ! 46 46 USE in_out_manager ! I/O manager … … 57 57 PRIVATE 58 58 59 PUBLIC dyn_atf ! routine called by step.F9059 PUBLIC dyn_atf_lf ! routine called by step.F90 60 60 61 61 !! * Substitutions … … 63 63 !!---------------------------------------------------------------------- 64 64 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 65 !! $Id$ 65 !! $Id$ 66 66 !! Software governed by the CeCILL license (see ./LICENSE) 67 67 !!---------------------------------------------------------------------- 68 68 CONTAINS 69 69 70 SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )70 SUBROUTINE dyn_atf_lf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 71 71 !!---------------------------------------------------------------------- 72 !! *** ROUTINE dyn_atf ***73 !! 74 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 72 !! *** ROUTINE dyn_atf_lf *** 73 !! 74 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 75 75 !! condition on the after velocity and apply the Asselin time 76 76 !! filter to the now fields. … … 79 79 !! estimate (ln_dynspg_ts=T) 80 80 !! 81 !! * Apply lateral boundary conditions on after velocity 81 !! * Apply lateral boundary conditions on after velocity 82 82 !! at the local domain boundaries through lbc_lnk call, 83 83 !! at the one-way open boundaries (ln_bdy=T), … … 86 86 !! * Apply the Asselin time filter to the now fields 87 87 !! arrays to start the next time step: 88 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 88 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 89 89 !! + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 90 90 !! Note that with flux form advection and non linear free surface, 91 91 !! the time filter is applied on thickness weighted velocity. 92 !! As a result, dyn_atf MUST be called after tra_atf.93 !! 94 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 92 !! As a result, dyn_atf_lf MUST be called after tra_atf. 93 !! 94 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 95 95 !!---------------------------------------------------------------------- 96 96 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 102 102 REAL(wp) :: zue3a, zue3n, zue3b, zcoef ! local scalars 103 103 REAL(wp) :: zve3a, zve3n, zve3b, z1_2dt ! - - 104 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve , zwfld105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z e3t_f, ze3u_f, ze3v_f, zua, zva104 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zua, zva 106 106 !!---------------------------------------------------------------------- 107 107 ! 108 IF( ln_timing ) CALL timing_start('dyn_atf ')108 IF( ln_timing ) CALL timing_start('dyn_atf_lf') 109 109 IF( ln_dynspg_ts ) ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) ) 110 110 IF( l_trddyn ) ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) ) … … 112 112 IF( kt == nit000 ) THEN 113 113 IF(lwp) WRITE(numout,*) 114 IF(lwp) WRITE(numout,*) 'dyn_atf : Asselin time filtering'114 IF(lwp) WRITE(numout,*) 'dyn_atf_lf : Asselin time filtering' 115 115 IF(lwp) WRITE(numout,*) '~~~~~~~' 116 116 ENDIF … … 131 131 ! 132 132 IF( .NOT.ln_bt_fw ) THEN 133 ! Remove advective velocity from "now velocities" 134 ! prior to asselin filtering 135 ! In the forward case, this is done below after asselin filtering 136 ! so that asselin contribution is removed at the same time 133 ! Remove advective velocity from "now velocities" 134 ! prior to asselin filtering 135 ! In the forward case, this is done below after asselin filtering 136 ! so that asselin contribution is removed at the same time 137 137 DO jk = 1, jpkm1 138 138 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 139 139 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 140 END DO 140 END DO 141 141 ENDIF 142 142 ENDIF 143 143 144 144 ! Update after velocity on domain lateral boundaries 145 ! -------------------------------------------------- 145 ! -------------------------------------------------- 146 146 # if defined key_agrif 147 147 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 148 148 # endif 149 149 ! 150 CALL lbc_lnk_multi( 'dynatf ', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. ) !* local domain boundaries150 CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. ) !* local domain boundaries 151 151 ! 152 152 ! !* BDY open boundaries … … 157 157 ! 158 158 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 160 160 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt 161 161 ! … … 177 177 ! Time filter and swap of dynamics arrays 178 178 ! ------------------------------------------ 179 180 IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN !* Leap-Frog : Asselin time filter 179 180 IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN !* Leap-Frog : Asselin time filter 181 181 ! ! =============! 182 182 IF( ln_linssh ) THEN ! Fixed volume ! … … 191 191 ! Time-filtered scale factor at t-points 192 192 ! ---------------------------------------------------- 193 ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) ) 194 DO jk = 1, jpkm1 195 ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) ) 196 END DO 197 ! Add volume filter correction: compatibility with tracer advection scheme 198 ! => time filter + conservation correction 199 zcoef = atfp * rdt * r1_rau0 200 zwfld(:,:) = emp_b(:,:) - emp(:,:) 201 IF ( ln_rnf ) zwfld(:,:) = zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) ) 202 DO jk = 1, jpkm1 203 ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 204 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 193 DO jk = 1, jpk ! filtered scale factor at T-points 194 pe3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t_f(:,:) * tmask(:,:,jk) ) 205 195 END DO 206 196 ! 207 ! ice shelf melting (deal separately as it can be in depth)208 ! PM: we could probably define a generic subroutine to do the in depth correction209 ! to manage rnf, isf and possibly in the futur icb, tide water glacier (...)210 ! ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)211 IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rdt )212 !213 pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1) ! filtered scale factor at T-points214 197 ! 215 198 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 216 199 ! Before filtered scale factor at (u/v)-points 217 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 218 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 200 DO jk = 1, jpk 201 pe3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) ) 202 pe3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) ) 203 END DO 204 ! 219 205 DO_3D_11_11( 1, jpkm1 ) 220 206 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) … … 224 210 ELSE ! Asselin filter applied on thickness weighted velocity 225 211 ! 226 ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )227 ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f228 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' )229 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' )230 212 DO_3D_11_11( 1, jpkm1 ) 231 213 zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) … … 235 217 zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb) 236 218 zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb) 219 ! ! filtered scale factor at U-,V-points 220 pe3u(ji,jj,jk,Kmm) = e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) * umask(ji,jj,jk) ) 221 pe3v(ji,jj,jk,Kmm) = e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) * vmask(ji,jj,jk) ) 237 222 ! 238 puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)239 pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)223 puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / pe3u(ji,jj,jk,Kmm) !!st ze3u_f(ji,jj,jk) 224 pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / pe3v(ji,jj,jk,Kmm) !!st ze3v_f(ji,jj,jk) 240 225 END_3D 241 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 242 pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 243 ! 244 DEALLOCATE( ze3u_f , ze3v_f ) 226 ! 245 227 ENDIF 246 228 ! 247 DEALLOCATE( ze3t_f, zwfld )248 229 ENDIF 249 230 ! 250 231 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 251 232 ! Revert filtered "now" velocities to time split estimate 252 ! Doing it here also means that asselin filter contribution is removed 233 ! Doing it here also means that asselin filter contribution is removed 253 234 zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 254 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 235 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 255 236 DO jk = 2, jpkm1 256 237 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 257 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 238 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 258 239 END DO 259 240 DO jk = 1, jpkm1 … … 307 288 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, & 308 289 & tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask ) 309 ! 290 ! 310 291 IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) 311 292 IF( l_trddyn ) DEALLOCATE( zua, zva ) 312 IF( ln_timing ) CALL timing_stop('dyn_atf ')313 ! 314 END SUBROUTINE dyn_atf 293 IF( ln_timing ) CALL timing_stop('dyn_atf_lf') 294 ! 295 END SUBROUTINE dyn_atf_lf 315 296 316 297 !!========================================================================= 317 END MODULE dynatf 298 END MODULE dynatfLF -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatf.F90
r12377 r12581 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 … … 69 69 !! *** ROUTINE traatf *** 70 70 !! 71 !! ** Purpose : Apply the boundary condition on the after temperature 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 … … 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 … … 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 surface 158 IF( ln_linssh ) THEN ; CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000, 'TRA', pts, jpts ) ! linear free surface 159 159 ELSE ; CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 160 160 ENDIF … … 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 … … 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 !! … … 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 ! … … 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) ] ) … … 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) … … 296 296 ! 297 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) ) 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 301 IF ( ll_isf ) THEN 302 302 IF ( ln_isfcav_mlt ) ze3t_f = ze3t_f - zfact2 * zscale * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) ) … … 304 304 ENDIF 305 305 ! 306 IF( jk == mikt(ji,jj) ) THEN ! first level 306 IF( jk == mikt(ji,jj) ) THEN ! first level 307 307 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 308 308 ENDIF 309 309 ! 310 310 ! 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) ) 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) ) 313 313 ! 314 314 ! 315 315 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) ) & 316 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 317 317 & * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) 318 318 … … 328 328 & * e3t(ji,jj,jk,Kmm) / rhisf_tbl_cav(ji,jj) 329 329 END IF 330 ! level partially include in Losch_2008 ice shelf boundary layer 330 ! level partially include in Losch_2008 ice shelf boundary layer 331 331 IF ( jk == misfkb_cav(ji,jj) ) THEN 332 332 ztc_f = ztc_f - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) & … … 342 342 & * e3t(ji,jj,jk,Kmm) / rhisf_tbl_par(ji,jj) 343 343 END IF 344 ! level partially include in Losch_2008 ice shelf boundary layer 344 ! level partially include in Losch_2008 ice shelf boundary layer 345 345 IF ( jk == misfkb_par(ji,jj) ) THEN 346 346 ztc_f = ztc_f - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) ) & … … 371 371 ! 372 372 END_3D 373 ! 373 ! 374 374 END DO 375 375 ! 376 376 IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) ) THEN 377 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 377 IF( l_trdtra .AND. cdtype == 'TRA' ) THEN 378 378 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) ) 379 379 CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) ) -
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 -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90
r12492 r12581 41 41 USE iom ! xIOs server 42 42 USE domqe 43 USE traatfLF ! time filtering (tra_atf_lf routine) 44 USE dynatfLF ! time filtering (dyn_atf_lf routine) 43 45 44 46 IMPLICIT NONE … … 286 288 !! 287 289 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 288 CALL tra_atf ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays289 CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors290 290 CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 291 CALL dom_qe_r3c ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh 292 CALL tra_atf_lf ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays 293 CALL dyn_atf_lf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 291 294 ! 292 295 ! Swap time levels
Note: See TracChangeset
for help on using the changeset viewer.