- Timestamp:
- 2020-03-26T15:59:52+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE
- Files:
-
- 1 added
- 1 copied
-
. (added)
-
step.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/step.F90) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
r12529 r12614 2 2 !!====================================================================== 3 3 !! *** MODULE step *** 4 !! Time-stepping : manager of the ocean, tracer and icetime stepping4 !! Time-stepping : manager of the shallow water equation time stepping 5 5 !!====================================================================== 6 !! History : OPA ! 1991-03 (G. Madec) Original code 7 !! - ! 1991-11 (G. Madec) 8 !! - ! 1992-06 (M. Imbard) add a first output record 9 !! - ! 1996-04 (G. Madec) introduction of dynspg 10 !! - ! 1996-04 (M.A. Foujols) introduction of passive tracer 11 !! 8.0 ! 1997-06 (G. Madec) new architecture of call 12 !! 8.2 ! 1997-06 (G. Madec, M. Imbard, G. Roullet) free surface 13 !! - ! 1999-02 (G. Madec, N. Grima) hpg implicit 14 !! - ! 2000-07 (J-M Molines, M. Imbard) Open Bondary Conditions 15 !! NEMO 1.0 ! 2002-06 (G. Madec) free form, suppress macro-tasking 16 !! - ! 2004-08 (C. Talandier) New trends organization 17 !! - ! 2005-01 (C. Ethe) Add the KPP closure scheme 18 !! - ! 2005-11 (G. Madec) Reorganisation of tra and dyn calls 19 !! - ! 2006-01 (L. Debreu, C. Mazauric) Agrif implementation 20 !! - ! 2006-07 (S. Masson) restart using iom 21 !! 3.2 ! 2009-02 (G. Madec, R. Benshila) reintroduicing z*-coordinate 22 !! - ! 2009-06 (S. Masson, G. Madec) TKE restart compatible with key_cpl 23 !! 3.3 ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 24 !! - ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 25 !! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal 26 !! 3.6 ! 2012-07 (J. Simeon, G. Madec. C. Ethe) Online coarsening of outputs 27 !! 3.6 ! 2014-04 (F. Roquet, G. Madec) New equations of state 28 !! 3.6 ! 2014-10 (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 29 !! 3.7 ! 2014-10 (G. Madec) LDF simplication 30 !! - ! 2014-12 (G. Madec) remove KPP scheme 31 !! - ! 2015-11 (J. Chanut) free surface simplification (remove filtered free surface) 32 !! 4.0 ! 2017-05 (G. Madec) introduction of the vertical physics manager (zdfphy) 33 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 34 !!---------------------------------------------------------------------- 35 36 !!---------------------------------------------------------------------- 37 !! stp : OPA system time-stepping 6 !! History : NEMO ! 2020-03 (A. Nasser, G. Madec) Original code from 4.0.2 7 !!---------------------------------------------------------------------- 8 9 !!---------------------------------------------------------------------- 10 !! stp : Shallow Water time-stepping 38 11 !!---------------------------------------------------------------------- 39 12 USE step_oce ! time stepping definition modules 13 USE phycst ! physical constants 14 USE usrdef_nam 40 15 ! 41 USE iom ! xIOs server 16 USE iom ! xIOs server 42 17 43 18 IMPLICIT NONE … … 45 20 46 21 PUBLIC stp ! called by nemogcm.F90 47 22 48 23 !!---------------------------------------------------------------------- 49 24 !! time level indices 50 25 !!---------------------------------------------------------------------- 51 INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !! used by nemo_init 52 26 INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !! used by nemo_init 27 28 !! * Substitutions 29 # include "do_loop_substitute.h90" 53 30 !!---------------------------------------------------------------------- 54 31 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 68 45 !! *** ROUTINE stp *** 69 46 !! 70 !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.) 71 !! - Time stepping of SI3 (dynamic and thermodynamic eqs.) 72 !! - Time stepping of TRC (passive tracer eqs.) 47 !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.) 73 48 !! 74 !! ** Method : -1- Update forcings and data 75 !! -2- Update ocean physics 76 !! -3- Compute the t and s trends 77 !! -4- Update t and s 78 !! -5- Compute the momentum trends 79 !! -6- Update the horizontal velocity 80 !! -7- Compute the diagnostics variables (rd,N2, hdiv,w) 81 !! -8- Outputs and diagnostics 49 !! ** Method : -1- Update forcings 50 !! -2- Update the ssh at Naa 51 !! -3- Compute the momentum trends (Nrhs) 52 !! -4- Update the horizontal velocity 53 !! -5- Apply Asselin time filter to uu,vv,ssh 54 !! -6- Outputs and diagnostics 82 55 !!---------------------------------------------------------------------- 83 56 INTEGER :: ji, jj, jk ! dummy loop indice … … 85 58 !!gm kcall can be removed, I guess 86 59 INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt) 60 REAL(wp):: z1_2rho0 ! local scalars 61 62 REAL(wp) :: zue3a, zue3n, zue3b ! local scalars 63 REAL(wp) :: zve3a, zve3n, zve3b ! - - 64 REAL(wp) :: ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 87 65 !! --------------------------------------------------------------------- 88 66 #if defined key_agrif … … 105 83 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 106 84 ! 107 IF( l_1st_euler ) THEN 108 ! start or restart with Euler 1st time-step 109 rDt = rn_Dt 85 IF( l_1st_euler ) THEN ! start or restart with Euler 1st time-step 86 rDt = rn_Dt 110 87 r1_Dt = 1._wp / rDt 111 88 ENDIF 89 90 IF ( kstp == nit000 ) ww(:,:,:) = 0._wp ! initialize vertical velocity one for all to zero 91 112 92 ! 113 93 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 132 112 IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 133 113 IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries 134 IF( ln_isf ) CALL isf_stp ( kstp, Nnn ) 135 CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice) 136 137 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 138 ! Update stochastic parameters and random T/S fluctuations 139 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 140 IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters 141 IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn) ) ! Random T/S fluctuations 114 CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition 142 115 143 116 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 144 117 ! Ocean physics update 145 118 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 146 ! THERMODYNAMICS147 CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points148 CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points149 CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency150 CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency151 152 ! VERTICAL PHYSICS153 CALL zdf_phy( kstp, Nbb, Nnn, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)154 155 119 ! LATERAL PHYSICS 156 !157 IF( l_ldfslp ) THEN ! slope of lateral mixing158 CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) ) ! before in situ density159 160 IF( ln_zps .AND. .NOT. ln_isfcav) &161 & CALL zps_hde ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient162 & rhd, gru , grv ) ! of t, s, rd at the last ocean level163 164 IF( ln_zps .AND. ln_isfcav) &165 & CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)166 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level167 IF( ln_traldf_triad ) THEN168 CALL ldf_slp_triad( kstp, Nbb, Nnn ) ! before slope for triad operator169 ELSE170 CALL ldf_slp ( kstp, rhd, rn2b, Nbb, Nnn ) ! before slope for standard operator171 ENDIF172 ENDIF173 120 ! ! eddy diffusivity coeff. 174 IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp, Nbb, Nnn ) ! and/or eiv coeff.175 121 IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff. 176 122 … … 180 126 181 127 CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor) 128 182 129 IF( .NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors 183 CALL wzv ( kstp, Nbb, Nnn, ww, Naa ) ! now cross-level velocity 184 IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning 185 CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) ) ! now in situ density for hpg computation 186 187 130 188 131 uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero 189 132 vv(:,:,:,Nrhs) = 0._wp 190 133 191 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) &192 & CALL dyn_asm_inc ( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! apply dynamics assimilation increment193 134 IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs ) ! bdy damping trends 135 194 136 #if defined key_agrif 195 137 IF(.NOT. Agrif_Root()) & … … 197 139 #endif 198 140 CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS 141 199 142 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 143 200 144 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing 201 IF( ln_zdfosm ) CALL dyn_osm( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes ==> RHS 202 CALL dyn_hpg( kstp, Nnn , uu, vv, Nrhs ) ! horizontal gradient of Hydrostatic pressure 203 CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa ) ! surface pressure gradient 204 205 ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well 206 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 207 CALL div_hor ( kstp, Nbb, Nnn ) ! Horizontal divergence (2nd call in time-split case) 208 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa, kcall=2 ) ! after vertical scale factors (update depth average component) 209 ENDIF 210 CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion 211 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 212 CALL wzv ( kstp, Nbb, Nnn, ww, Naa ) ! now cross-level velocity 213 IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning 214 ENDIF 215 216 217 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 218 ! cool skin 219 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 220 IF ( ln_diurnal ) CALL diurnal_layers( kstp ) 221 222 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 223 ! diagnostics and outputs 224 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 225 IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats 226 IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics 227 CALL dia_hth ( kstp, Nnn ) ! Thermocline depth (20 degres isotherm depth) 228 IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports 229 CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag 230 CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics 231 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 232 IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output 233 IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics 234 IF( lk_diamlr ) CALL dia_mlr ! Update time used in multiple-linear-regression analysis 235 236 #if defined key_top 237 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 238 ! Passive Tracer Model 239 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 240 CALL trc_stp ( kstp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 241 #endif 242 243 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 244 ! Active tracers 245 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 246 ts(:,:,:,:,Nrhs) = 0._wp ! set tracer trends to zero 247 248 IF( lk_asminc .AND. ln_asmiau .AND. & 249 & ln_trainc ) CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs ) ! apply tracer assimilation increment 250 CALL tra_sbc ( kstp, Nnn, ts, Nrhs ) ! surface boundary condition 251 IF( ln_traqsr ) CALL tra_qsr ( kstp, Nnn, ts, Nrhs ) ! penetrative solar radiation qsr 252 IF( ln_isf ) CALL tra_isf ( kstp, Nnn, ts, Nrhs ) ! ice shelf heat flux 253 IF( ln_trabbc ) CALL tra_bbc ( kstp, Nnn, ts, Nrhs ) ! bottom heat flux 254 IF( ln_trabbl ) CALL tra_bbl ( kstp, Nbb, Nnn, ts, Nrhs ) ! advective (and/or diffusive) bottom boundary layer scheme 255 IF( ln_tradmp ) CALL tra_dmp ( kstp, Nbb, Nnn, ts, Nrhs ) ! internal damping trends 256 IF( ln_bdy ) CALL bdy_tra_dmp( kstp, Nbb, ts, Nrhs ) ! bdy damping trends 257 #if defined key_agrif 258 IF(.NOT. Agrif_Root()) & 259 & CALL Agrif_Sponge_tra ! tracers sponge 260 #endif 261 CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS 262 IF( ln_zdfosm ) CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS 263 IF( lrst_oce .AND. ln_zdfosm ) & 264 & CALL osm_rst ( kstp, Nnn, 'WRITE' ) ! write OSMOSIS outputs + ww (so must do here) to restarts 265 CALL tra_ldf ( kstp, Nbb, Nnn, ts, Nrhs ) ! lateral mixing 266 267 CALL tra_zdf ( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing and after tracer fields 268 IF( ln_zdfnpc ) CALL tra_npc ( kstp, Nnn, Nrhs, ts, Naa ) ! update after fields by non-penetrative convection 145 146 !!an - calcul du gradient de pression horizontal (explicit) 147 DO_3D_00_00( 1, jpkm1 ) 148 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 149 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 150 END_3D 151 ! 152 ! add wind stress forcing and layer linear friction to the RHS 153 z1_2rho0 = 0.5_wp * r1_rho0 154 DO_3D_00_00(1,jpkm1) 155 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) & 156 & - rn_rfr * uu(ji,jj,jk,Nbb) 157 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) & 158 & - rn_rfr * vv(ji,jj,jk,Nbb) 159 END_3D 160 !!an 161 162 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 163 ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3 164 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 165 166 !!an futur module dyn_nxt (a la place de dyn_atf) 167 168 IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity 169 IF( l_1st_euler ) THEN ! Euler time stepping (no Asselin filter) 170 DO_3D_00_00(1,jpkm1) 171 uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 172 vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 173 END_3D 174 ELSE ! Leap Frog time stepping + Asselin filter 175 DO_3D_11_11(1,jpkm1) 176 zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 177 zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 178 ! ! Asselin time filter on u,v (Nnn) 179 uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua) 180 vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva) 181 ! 182 ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn) + e3u(ji,jj,jk,Naa) ) 183 ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn) + e3v(ji,jj,jk,Naa) ) 184 ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn) + e3t(ji,jj,jk,Naa) ) 185 ! 186 e3u(ji,jj,jk,Nnn) = ze3u_tf 187 e3v(ji,jj,jk,Nnn) = ze3v_tf 188 e3t(ji,jj,jk,Nnn) = ze3t_tf 189 ! 190 uu(ji,jj,jk,Naa) = zua 191 vv(ji,jj,jk,Naa) = zva 192 END_3D 193 ENDIF 194 ! 195 ELSE ! flux form : applied on thickness weighted velocity 196 IF( l_1st_euler ) THEN ! Euler time stepping (no Asselin filter) 197 DO_3D_00_00(1,jpkm1) 198 zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 199 zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 200 ! ! LF time stepping 201 zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 202 zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 203 ! 204 uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa) 205 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 206 END_3D 207 ELSE ! Leap Frog time stepping + Asselin filter 208 DO_3D_11_11(1,jpkm1) 209 zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn) 210 zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn) 211 zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 212 zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 213 ! ! LF time stepping 214 zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 215 zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 216 ! ! Asselin time filter on e3u/v/t 217 ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn) + e3u(ji,jj,jk,Naa) ) 218 ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn) + e3v(ji,jj,jk,Naa) ) 219 ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn) + e3t(ji,jj,jk,Naa) ) 220 ! ! Asselin time filter on u,v (Nnn) 221 uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_tf 222 vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_tf 223 ! 224 e3u(ji,jj,jk,Nnn) = ze3u_tf 225 e3v(ji,jj,jk,Nnn) = ze3v_tf 226 e3t(ji,jj,jk,Nnn) = ze3t_tf 227 ! 228 uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa) 229 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 230 END_3D 231 ENDIF 232 ENDIF 233 234 CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., & !* local domain boundaries 235 & uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. ) 236 237 !!an 269 238 270 239 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 271 240 ! Set boundary conditions, time filter and swap time levels 272 241 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 273 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 274 !! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 275 !! If so: 276 !! (i) no need to call agrif update at initialization time 277 !! (ii) no need to update "before" fields 278 !! 279 !! Apart from creating new tra_swp/dyn_swp routines, this however: 280 !! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 281 !! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 282 !! e.g. a shift of the feedback interface inside child domain. 283 !! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 284 !! place. 285 !! 286 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 287 CALL tra_atf ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays 288 CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 289 CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 290 ! 242 243 !!an TO BE ADDED : dyn_nxt 244 !! CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 245 !!an TO BE ADDED : a simplifier 246 ! CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 247 248 IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps 249 ! ! filtering "now" field 250 ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) ) 251 ENDIF 252 253 !!an 254 255 291 256 ! Swap time levels 292 257 Nrhs = Nbb … … 295 260 Naa = Nrhs 296 261 ! 297 IF(.NOT.ln_linssh) CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 298 ! 299 IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics 300 301 !!gm : This does not only concern the dynamics ==>>> add a new title 302 !!gm2: why ouput restart before AGRIF update? 303 !! 304 !!jc: That would be better, but see comment above 305 !! 262 CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 263 264 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 265 ! diagnostics and outputs 266 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 267 IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats 268 IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics 269 270 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 271 272 ! 306 273 IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file 307 IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters274 308 275 309 276 #if defined key_agrif … … 324 291 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 325 292 CALL stp_ctl ( kstp, Nbb, Nnn, indic ) 293 326 294 327 295 IF( kstp == nit000 ) THEN ! 1st time step only … … 331 299 ENDIF 332 300 333 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>334 ! Coupled mode335 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<336 !!gm why lk_oasis and not lk_cpl ????337 IF( lk_oasis ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges338 301 ! 339 302 #if defined key_iomput … … 341 304 CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 342 305 IF(lrxios) CALL iom_context_finalize( crxios_context ) 343 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !344 306 ENDIF 345 307 #endif
Note: See TracChangeset
for help on using the changeset viewer.
