- 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/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
Note: See TracChangeset
for help on using the changeset viewer.
