- Timestamp:
- 2018-07-26T09:50:51+02:00 (6 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpRK3.F90
r9939 r10001 1 MODULE st ep1 MODULE stpRK3 2 2 !!====================================================================== 3 !! *** MODULE st ep***4 !! Time-stepping: manager of the ocean, tracer and ice time stepping3 !! *** MODULE stpRK3 *** 4 !! RK3 Time-stepping: manager of the ocean, tracer and ice 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 !!---------------------------------------------------------------------- 34 35 !!---------------------------------------------------------------------- 36 !! stp : NEMO system time-stepping 6 !! History : 5.0 ! 2018-07 (G. Madec) original code 7 !!---------------------------------------------------------------------- 8 9 !!---------------------------------------------------------------------- 10 !! stp_RK3 : NEMO system RK3 time-stepping scheme 11 !! stp_RK3_init : initialize the RK3 scheme 37 12 !!---------------------------------------------------------------------- 38 13 USE step_oce ! time stepping definition modules … … 43 18 PRIVATE 44 19 45 PUBLIC stp ! called by nemogcm.F90 46 47 !!---------------------------------------------------------------------- 48 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 20 PUBLIC stp_RK3 ! called by nemogcm.F90 21 22 LOGICAL :: l_1st_stg = .TRUE. ! 1st stage only flag 23 LOGICAL :: l_2nd_stg = .TRUE. ! 2nd stage only flag 24 LOGICAL :: l_3rd_stg = .TRUE. ! 3rd stage only flag 25 26 !!---------------------------------------------------------------------- 27 !! NEMO/OCE 5.0 , NEMO Consortium (2018) 49 28 !! $Id$ 50 29 !! Software governed by the CeCILL licence (./LICENSE) … … 53 32 54 33 #if defined key_agrif 55 RECURSIVE SUBROUTINE stp ( )34 RECURSIVE SUBROUTINE stp_RK3( ) 56 35 INTEGER :: kstp ! ocean time-step index 57 36 #else 58 SUBROUTINE stp ( kstp )37 SUBROUTINE stp_RK3( kstp ) 59 38 INTEGER, INTENT(in) :: kstp ! ocean time-step index 60 39 #endif … … 75 54 !! -8- Outputs and diagnostics 76 55 !!---------------------------------------------------------------------- 77 INTEGER :: ji, jj, jk ! dummy loop indice78 INTEGER :: indic ! error indicator if < 056 INTEGER :: ji, jj, jk, jstg ! dummy loop indice 57 INTEGER :: indic ! error indicator if < 0 79 58 !!gm kcall can be removed, I guess 80 59 INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt) … … 93 72 ! 94 73 IF( ln_timing ) CALL timing_start('stp') 95 ! 74 75 !!!======================!!! 76 !!! First STAGE only !!! 77 !!!======================!!! 78 96 79 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 97 80 ! update I/O and calendar … … 114 97 IF( ln_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 115 98 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 99 IF ( ln_diurnal ) CALL stp_diurnal( kstp ) ! diagnose cool skin 100 ! 116 101 117 102 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 118 103 ! Update stochastic parameters and random T/S fluctuations 119 104 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 120 IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters121 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations105 IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters 106 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 122 107 123 108 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 155 140 IF( l_ldfdyn_time ) CALL ldf_dyn( kstp ) ! eddy viscosity coeff. 156 141 157 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 158 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 159 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 160 161 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 162 IF(.NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 163 CALL wzv ( kstp ) ! now cross-level velocity 164 CALL eos ( tsn, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation 142 143 !!!======================!!! 144 !!! Loop over stages !!! 145 !!!======================!!! 146 147 DO jstg = 1, 3 148 149 SELECT CASE( jstg ) 150 CASE( 1 ) ; rDt = rn_Dt / 3._wp 151 CASE( 2 ) ; rDt = rn_Dt / 2._wp 152 CASE( 3 ) ; rDt = rn_Dt 153 END SELECT 154 155 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 156 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 157 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 158 159 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 160 IF(.NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 161 CALL wzv ( kstp ) ! now cross-level velocity 162 CALL eos ( tsb, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation 165 163 166 !!jc: fs simplification 167 !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636) 168 !! but ensures reproductible results 169 !! with previous versions using split-explicit free surface 170 IF( ln_zps .AND. .NOT. ln_isfcav ) & 171 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 172 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 173 IF( ln_zps .AND. ln_isfcav ) & 174 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 175 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 176 !!jc: fs simplification 177 178 ua(:,:,:) = 0._wp ! set dynamics trends to zero 179 va(:,:,:) = 0._wp 180 181 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) & 182 & CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 183 IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends 184 #if defined key_agrif 185 IF(.NOT. Agrif_Root()) & 186 & CALL Agrif_Sponge_dyn ! momentum sponge 187 #endif 188 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 189 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 190 CALL dyn_ldf ( kstp ) ! lateral mixing 191 IF( ln_zdfosm ) CALL dyn_osm ( kstp ) ! OSMOSIS non-local velocity fluxes 192 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 193 CALL dyn_spg ( kstp ) ! surface pressure gradient 194 195 ! With split-explicit free surface, since now transports have been updated and ssha as well 196 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 197 CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case) 198 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 199 CALL wzv ( kstp ) ! now cross-level velocity 200 ENDIF 201 202 CALL dyn_zdf ( kstp ) ! vertical diffusion 203 204 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 205 ! cool skin 206 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 207 IF ( ln_diurnal ) CALL stp_diurnal( kstp ) 208 209 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 210 ! diagnostics and outputs 211 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 212 IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats 213 IF( ln_diacfl ) CALL dia_cfl ( kstp ) ! Courant number diagnostics 214 IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 215 IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports 216 CALL dia_ar5 ( kstp ) ! ar5 diag 217 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 218 CALL dia_wri ( kstp ) ! ocean model: outputs 219 ! 220 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output 221 164 !!jc: fs simplification 165 !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636) 166 !! but ensures reproductible results 167 !! with previous versions using split-explicit free surface 168 IF( ln_zps .AND. .NOT. ln_isfcav ) & 169 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 170 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 171 IF( ln_zps .AND. ln_isfcav ) & 172 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 173 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 174 !!jc: fs simplification 175 176 ua (:,:,:) = 0._wp ! set the RHS of dyn Eq. to zero 177 va (:,:,:) = 0._wp 178 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 179 ! 180 ! ! ================ ! 181 IF ( jstg <= 3 ) THEN ! stages 1 & 2 : ADV, COR, HPG and SPG trends only 182 ! ! ================ ! 183 ! 184 ! !== dynamics ==! 185 CALL dyn_adv( kstp ) ! advection (vector or flux form) 186 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 187 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 188 CALL dyn_spg( kstp ) ! surface pressure gradient 189 ! With split-explicit free surface, since now transports have been updated and ssha as well 190 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 191 CALL div_hor( kstp ) ! Horizontal divergence (2nd call in time-split case) 192 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 193 CALL wzv ( kstp ) ! now cross-level velocity 194 ENDIF 195 !!gm to be added here : time stepping ==>>> un & vn 196 197 198 ! !== tracers ==! 222 199 #if defined key_top 223 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 224 ! Passive Tracer Model 225 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 200 CALL trc_adv( kstp ) ! horizontal & vertical advection 201 #endif 202 CALL tra_adv( kstp ) ! horizontal & vertical advection 203 204 !!gm to be added here : time stepping ==>>> tsn 205 206 ! 207 ! ! ================ ! 208 ELSE ! stage 3 : add all dynamical trends 209 ! ! ================ ! 210 ! 211 ! !== dynamics ==! 212 CALL dyn_adv( kstp ) ! advection (vector or flux form) 213 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 214 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 215 ! 216 IF(l_dynasm) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 217 IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends 218 #if defined key_agrif 219 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 220 #endif 221 IF( ln_zdfosm ) CALL dyn_osm( kstp ) ! OSMOSIS non-local velocity fluxes 222 CALL dyn_ldf( kstp ) ! lateral mixing 223 CALL dyn_spg( kstp ) ! surface pressure gradient 224 ! With split-explicit free surface, since now transports have been updated and ssha as well 225 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 226 CALL div_hor( kstp ) ! Horizontal divergence (2nd call in time-split case) 227 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 228 CALL wzv ( kstp ) ! now cross-level velocity 229 ENDIF 230 CALL dyn_zdf( kstp ) ! vertical diffusion & time-stepping 231 232 ! !== tracers ==! 233 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 234 #if defined key_top 226 235 CALL trc_stp ( kstp ) ! time-stepping 227 236 #endif 228 237 229 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>230 ! Active tracers231 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<232 238 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 233 239 234 IF( lk_asminc .AND. ln_asmiau .AND. &235 & ln_trainc) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment240 CALL tra_adv ( kstp ) ! horizontal & vertical advection 241 IF( l_traasm ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment 236 242 CALL tra_sbc ( kstp ) ! surface boundary condition 237 243 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr … … 241 247 IF( ln_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends 242 248 #if defined key_agrif 243 IF(.NOT. Agrif_Root()) & 244 & CALL Agrif_Sponge_tra ! tracers sponge 245 #endif 246 CALL tra_adv ( kstp ) ! horizontal & vertical advection 249 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 250 #endif 247 251 IF( ln_zdfosm ) CALL tra_osm ( kstp ) ! OSMOSIS non-local tracer fluxes 248 252 IF( lrst_oce .AND. ln_zdfosm ) & … … 255 259 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 256 260 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 261 262 263 ENDIF 264 265 266 267 257 268 258 269 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 278 289 IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 279 290 ! 291 292 293 294 295 !!!==========================!!! 296 !!! end Loop over stages !!! 297 !!!==========================!!! 298 299 END DO 300 301 280 302 IF( ln_diahsb ) CALL dia_hsb ( kstp ) ! - ML - global conservation diagnostics 281 303 … … 296 318 IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components 297 319 #endif 298 IF( ln_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 320 321 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 322 ! diagnostics and outputs 323 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 324 IF( ln_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 325 IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats 326 IF( ln_diacfl ) CALL dia_cfl ( kstp ) ! Courant number diagnostics 327 IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 328 IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports 329 CALL dia_ar5 ( kstp ) ! ar5 diag 330 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 331 CALL dia_wri ( kstp ) ! ocean model: outputs 299 332 300 333 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 323 356 #endif 324 357 ! 325 IF( l_1st_euler ) THEN326 rDt = 2._wp * rn_Dt ! recover Leap-frog time-step327 r1_Dt = 1._wp / rDt328 l_1st_euler = .FALSE.329 ENDIF330 !331 358 IF( ln_timing ) CALL timing_stop('stp') 332 359 ! 333 END SUBROUTINE stp 360 END SUBROUTINE stp_RK3 361 362 363 SUBROUTINE stp_RK3_init 364 !!---------------------------------------------------------------------- 365 !! *** ROUTINE stp_RK3_init *** 366 !! 367 !! ** Purpose : RK3 time stepping initialization 368 !! 369 !! ** Method : 370 !!---------------------------------------------------------------------- 371 372 373 END SUBROUTINE stp_RK3_init 334 374 335 375 !!====================================================================== 336 END MODULE st ep376 END MODULE stpRK3
Note: See TracChangeset
for help on using the changeset viewer.
