- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/step.F90
r4624 r6225 2 2 !!====================================================================== 3 3 !! *** MODULE step *** 4 !! Time-stepping 4 !! Time-stepping : manager of the ocean, tracer and ice time stepping 5 5 !!====================================================================== 6 6 !! History : OPA ! 1991-03 (G. Madec) Original code … … 24 24 !! - ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 25 25 !! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal 26 !! ! 2012-07 (J. Simeon, G. Madec. C. Ethe) Online coarsening of outputs 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.7 ! 2014-10 (G. Madec) LDF simplication 29 !! - ! 2014-12 (G. Madec) remove KPP scheme 30 !! - ! 2015-11 (J. Chanut) free surface simplification 27 31 !!---------------------------------------------------------------------- 28 32 … … 31 35 !!---------------------------------------------------------------------- 32 36 USE step_oce ! time stepping definition modules 37 ! 38 USE iom ! xIOs server 33 39 34 40 IMPLICIT NONE 35 41 PRIVATE 36 42 37 PUBLIC stp ! called by opa.F90 38 39 !! * Substitutions 40 # include "domzgr_substitute.h90" 41 # include "zdfddm_substitute.h90" 42 !!---------------------------------------------------------------------- 43 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 43 PUBLIC stp ! called by nemogcm.F90 44 45 !!---------------------------------------------------------------------- 46 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 44 47 !! $Id$ 45 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 48 51 49 52 #if defined key_agrif 50 SUBROUTINE stp( )53 RECURSIVE SUBROUTINE stp( ) 51 54 INTEGER :: kstp ! ocean time-step index 52 55 #else … … 67 70 !! -5- Compute the momentum trends 68 71 !! -6- Update the horizontal velocity 69 !! -7- Compute the diagnostics variables (rd,N2, div,cur,w)72 !! -7- Compute the diagnostics variables (rd,N2, hdiv,w) 70 73 !! -8- Outputs and diagnostics 71 74 !!---------------------------------------------------------------------- … … 74 77 INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt) 75 78 !! --------------------------------------------------------------------- 76 77 79 #if defined key_agrif 78 80 kstp = nit000 + Agrif_Nb_Step() 79 ! IF ( Agrif_Root() .and. lwp) Write(*,*) '---' 80 ! IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp 81 IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE. 81 IF( lk_agrif_debug ) THEN 82 IF( Agrif_Root() .and. lwp) WRITE(*,*) '---' 83 IF(lwp) WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint() 84 ENDIF 85 IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE. 82 86 # if defined key_iomput 83 IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( "nemo")87 IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context ) 84 88 # endif 85 89 #endif 86 indic = 0 ! reset to no error condition 87 IF( kstp == nit000 ) THEN 88 CALL iom_init( "nemo" ) ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 89 IF( ln_crs ) CALL iom_init( "nemo_crs" ) ! initialize context for coarse grid 90 ENDIF 91 90 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 91 ! update I/O and calendar 92 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 93 indic = 0 ! reset to no error condition 94 95 IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 96 CALL iom_init( cxios_context ) ! for model grid (including passible AGRIF zoom) 97 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid 98 ENDIF 92 99 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) 93 CALL iom_setkt( kstp - nit000 + 1, "nemo" ) ! say to iom that we are at time step kstp 94 IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, "nemo_crs" ) ! say to iom that we are at time step kstp 95 96 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 97 ! Update data, open boundaries, surface boundary condition (including sea-ice) 98 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 99 IF( lk_tide ) CALL sbc_tide( kstp ) 100 CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp 101 IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp 102 103 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 104 ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice) 105 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 106 IF( lk_tide ) CALL sbc_tide( kstp ) ! update tide potential 107 IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 100 108 IF( lk_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 101 102 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 103 ! clem: moved here for bdy ice purpose 109 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 110 111 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 112 ! Update stochastic parameters and random T/S fluctuations 113 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 114 CALL sto_par( kstp ) ! Stochastic parameters 104 115 105 116 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 106 117 ! Ocean physics update (ua, va, tsa used as workspace) 107 118 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 108 CALL bn2( tsb, rn2b ) ! before Brunt-Vaisala frequency 109 CALL bn2( tsn, rn2 ) ! now Brunt-Vaisala frequency 119 ! THERMODYNAMICS 120 CALL eos_rab( tsb, rab_b ) ! before local thermal/haline expension ratio at T-points 121 CALL eos_rab( tsn, rab_n ) ! now local thermal/haline expension ratio at T-points 122 CALL bn2 ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 123 CALL bn2 ( tsn, rab_n, rn2 ) ! now Brunt-Vaisala frequency 124 110 125 ! 111 126 ! VERTICAL PHYSICS … … 115 130 IF( lk_zdftke ) CALL zdf_tke( kstp ) ! TKE closure scheme for Kz 116 131 IF( lk_zdfgls ) CALL zdf_gls( kstp ) ! GLS closure scheme for Kz 117 IF( lk_zdfkpp ) CALL zdf_kpp( kstp ) ! KPP closure scheme for Kz118 132 IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 119 avt (:,:,:) = rn_avt0 * tmask(:,:,:) 120 avmu(:,:,:) = rn_avm0 * umask(:,:,:) 121 avmv(:,:,:) = rn_avm0 * vmask(:,:,:) 122 ENDIF 133 avt (:,:,:) = rn_avt0 * wmask (:,:,:) 134 avmu(:,:,:) = rn_avm0 * wumask(:,:,:) 135 avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 136 ENDIF 137 123 138 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 124 DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2. e0* rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk) ; END DO139 DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk) ; END DO 125 140 ENDIF 126 141 IF( ln_zdfevd ) CALL zdf_evd( kstp ) ! enhanced vertical eddy diffusivity … … 128 143 IF( lk_zdftmx ) CALL zdf_tmx( kstp ) ! tidal vertical mixing 129 144 130 IF( lk_zdfddm .AND. .NOT. lk_zdfkpp ) & 131 & CALL zdf_ddm( kstp ) ! double diffusive mixing 145 IF( lk_zdfddm ) CALL zdf_ddm( kstp ) ! double diffusive mixing 132 146 133 147 CALL zdf_mxl( kstp ) ! mixed layer depth … … 139 153 ! LATERAL PHYSICS 140 154 ! 141 IF( lk_ldfslp ) THEN ! slope of lateral mixing 142 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 143 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 144 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 145 IF( ln_traldf_grif ) THEN ! before slope for Griffies operator 146 CALL ldf_slp_grif( kstp ) 147 ELSE 148 CALL ldf_slp( kstp, rhd, rn2b ) ! before slope for Madec operator 155 IF( l_ldfslp ) THEN ! slope of lateral mixing 156 !!gm : why this here ???? 157 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 158 !!gm 159 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 160 161 IF( ln_zps .AND. .NOT. ln_isfcav) & 162 & CALL zps_hde ( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 163 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 164 165 IF( ln_zps .AND. ln_isfcav) & 166 & CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 167 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 168 IF( ln_traldf_triad ) THEN 169 CALL ldf_slp_triad( kstp ) ! before slope for triad operator 170 ELSE 171 CALL ldf_slp ( kstp, rhd, rn2b ) ! before slope for standard operator 149 172 ENDIF 150 173 ENDIF 151 #if defined key_traldf_c2d 152 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient 153 #endif 154 #if defined key_traldf_c3d && key_traldf_smag 155 CALL ldf_tra_smag( kstp ) ! eddy induced velocity coefficient 156 # endif 157 #if defined key_dynldf_c3d && key_dynldf_smag 158 CALL ldf_dyn_smag( kstp ) ! eddy induced velocity coefficient 159 # endif 160 161 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 162 ! Ocean dynamics : hdiv, rot, ssh, e3, wn 163 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 164 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_cur) 165 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 166 CALL wzv ( kstp ) ! now cross-level velocity 167 168 IF( lk_dynspg_ts ) THEN 169 ! In case the time splitting case, update almost all momentum trends here: 170 ! Note that the computation of vertical velocity above, hence "after" sea level 171 ! is necessary to compute momentum advection for the rhs of barotropic loop: 172 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 173 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 174 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 175 176 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 177 va(:,:,:) = 0.e0 178 IF( ln_asmiau .AND. & 179 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 180 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! subtract Neptune velocities (simplified) 181 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 182 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 183 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 184 CALL dyn_ldf ( kstp ) ! lateral mixing 185 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! add Neptune velocities (simplified) 186 #if defined key_agrif 187 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 188 #endif 189 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 190 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 191 192 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated at this stage) 193 va_sv(:,:,:) = va(:,:,:) 194 195 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 196 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 197 CALL wzv ( kstp ) ! now cross-level velocity 198 ENDIF 199 174 ! ! eddy diffusivity coeff. and/or eiv coeff. 175 IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp ) 176 177 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 178 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 179 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 180 181 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 182 IF(.NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 183 CALL wzv ( kstp ) ! now cross-level velocity 184 !!gm : why also here ???? 185 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 186 !!gm 187 CALL eos ( tsn, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation 188 189 !!jc: fs simplification 190 !!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) 191 !! but ensures reproductible results 192 !! with previous versions using split-explicit free surface 193 IF( ln_zps .AND. .NOT. ln_isfcav ) & 194 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 195 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 196 IF( ln_zps .AND. ln_isfcav ) & 197 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 198 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 199 !!jc: fs simplification 200 201 ua(:,:,:) = 0._wp ! set dynamics trends to zero 202 va(:,:,:) = 0._wp 203 204 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) & 205 CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 206 IF( lk_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends 207 #if defined key_agrif 208 IF(.NOT. Agrif_Root()) & 209 & CALL Agrif_Sponge_dyn ! momentum sponge 210 #endif 211 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 212 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 213 CALL dyn_ldf ( kstp ) ! lateral mixing 214 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 215 CALL dyn_spg ( kstp ) ! surface pressure gradient 216 217 ! With split-explicit free surface, since now transports have been updated and ssha as well 218 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 219 CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case) 220 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 221 CALL wzv ( kstp ) ! now cross-level velocity 222 ENDIF 223 224 CALL dyn_bfr ( kstp ) ! bottom friction 225 CALL dyn_zdf ( kstp ) ! vertical diffusion 226 227 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 228 ! cool skin 229 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 230 IF ( ln_diurnal ) CALL stp_diurnal( kstp ) 231 200 232 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 201 233 ! diagnostics and outputs (ua, va, tsa used as workspace) 202 234 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 203 235 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 236 IF( nn_diacfl == 1 ) CALL dia_cfl( kstp ) ! Courant number diagnostics 204 237 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 205 IF( lk_diafwb ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 206 IF( ln_diaptr ) CALL dia_ptr( kstp ) ! Poleward TRansports diagnostics 238 IF(.NOT.ln_cpl ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 207 239 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 208 240 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag … … 210 242 CALL dia_wri( kstp ) ! ocean model: outputs 211 243 ! 212 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 213 214 244 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output 245 215 246 #if defined key_top 216 247 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 217 248 ! Passive Tracer Model 218 249 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 219 CALL trc_stp( kstp ) ! time-stepping 220 #endif 221 222 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 223 ! Active tracers (ua, va used as workspace) 224 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 225 tsa(:,:,:,:) = 0.e0 ! set tracer trends to zero 226 227 IF( ln_asmiau .AND. & 228 & ln_trainc ) CALL tra_asm_inc( kstp ) ! apply tracer assimilation increment 229 CALL tra_sbc ( kstp ) ! surface boundary condition 230 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 231 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 232 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 233 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 234 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 235 CALL tra_adv ( kstp ) ! horizontal & vertical advection 236 IF( lk_zdfkpp ) CALL tra_kpp ( kstp ) ! KPP non-local tracer fluxes 237 CALL tra_ldf ( kstp ) ! lateral mixing 238 #if defined key_agrif 239 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 240 #endif 241 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 242 243 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 244 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 245 CALL tra_nxt( kstp ) ! tracer fields at next time step 246 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 247 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv, & ! zps: time filtered hor. derivative 248 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 249 250 ELSE ! centered hpg (eos then time stepping) 251 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 252 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 253 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 254 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 255 ENDIF 256 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 257 CALL tra_nxt( kstp ) ! tracer fields at next time step 258 ENDIF 259 260 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 261 ! Dynamics (tsa used as workspace) 262 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 263 IF( lk_dynspg_ts ) THEN 264 ! revert to previously computed momentum tendencies 265 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 266 ua(:,:,:) = ua_sv(:,:,:) 267 va(:,:,:) = va_sv(:,:,:) 268 ! Revert now divergence and rotational to previously computed ones 269 !(needed because of the time swap in div_cur, at the beginning of each time step) 270 hdivn(:,:,:) = hdivb(:,:,:) 271 rotn(:,:,:) = rotb(:,:,:) 272 273 CALL dyn_bfr( kstp ) ! bottom friction 274 CALL dyn_zdf( kstp ) ! vertical diffusion 275 ELSE 276 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 277 va(:,:,:) = 0.e0 278 279 IF( ln_asmiau .AND. & 280 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 281 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 282 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 283 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 284 CALL dyn_adv( kstp ) ! advection (vector or flux form) 285 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 286 CALL dyn_ldf( kstp ) ! lateral mixing 287 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 288 #if defined key_agrif 289 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 290 #endif 291 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 292 CALL dyn_bfr( kstp ) ! bottom friction 293 CALL dyn_zdf( kstp ) ! vertical diffusion 294 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 295 ENDIF 296 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 297 298 CALL ssh_swp( kstp ) ! swap of sea surface height 299 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 300 301 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 302 IF( lk_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 303 304 IF( lrst_oce .AND. ln_diahsb ) CALL dia_hsb_rst( kstp, 'WRITE' ) 305 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 306 ! Control and restarts 307 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 308 CALL stp_ctl( kstp, indic ) 309 IF( indic < 0 ) THEN 310 CALL ctl_stop( 'step: indic < 0' ) 311 CALL dia_wri_state( 'output.abort', kstp ) 312 ENDIF 313 IF( kstp == nit000 ) THEN 250 CALL trc_stp ( kstp ) ! time-stepping 251 #endif 252 253 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 254 ! Active tracers 255 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 256 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 257 258 IF( lk_asminc .AND. ln_asmiau .AND. & 259 & ln_trainc ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment 260 CALL tra_sbc ( kstp ) ! surface boundary condition 261 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 262 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 263 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 264 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 265 IF( lk_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends 266 #if defined key_agrif 267 IF(.NOT. Agrif_Root()) & 268 & CALL Agrif_Sponge_tra ! tracers sponge 269 #endif 270 CALL tra_adv ( kstp ) ! horizontal & vertical advection 271 CALL tra_ldf ( kstp ) ! lateral mixing 272 273 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 274 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 275 !!gm 276 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 277 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 278 279 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 280 ! Set boundary conditions and Swap 281 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 282 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 283 !! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 284 !! If so: 285 !! (i) no need to call agrif update at initialization time 286 !! (ii) no need to update "before" fields 287 !! 288 !! Apart from creating new tra_swp/dyn_swp routines, this however: 289 !! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 290 !! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 291 !! e.g. a shift of the feedback interface inside child domain. 292 !! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 293 !! place. 294 !! 295 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine 296 CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap 297 CALL dyn_nxt ( kstp ) ! finalize (bcs) velocities at next time step and swap 298 CALL ssh_swp ( kstp ) ! swap of sea surface height 299 IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 300 ! 301 302 !!gm : This does not only concern the dynamics ==>>> add a new title 303 !!gm2: why ouput restart before AGRIF update? 304 !! 305 !!jc: That would be better, but see comment above 306 !! 307 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 308 309 #if defined key_agrif 310 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 311 ! AGRIF 312 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 313 CALL Agrif_Integrate_ChildGrids( stp ) 314 315 IF( Agrif_NbStepint() == 0 ) THEN ! AGRIF Update 316 !!jc in fact update is useless at last time step, but do it for global diagnostics 317 CALL Agrif_Update_Tra() ! Update active tracers 318 CALL Agrif_Update_Dyn() ! Update momentum 319 ENDIF 320 #endif 321 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 322 IF( ln_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 323 324 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 325 ! Control 326 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 327 CALL stp_ctl ( kstp, indic ) 328 IF( indic < 0 ) THEN 329 CALL ctl_stop( 'step: indic < 0' ) 330 CALL dia_wri_state( 'output.abort', kstp ) 331 ENDIF 332 IF( kstp == nit000 ) THEN 314 333 CALL iom_close( numror ) ! close input ocean restart file 315 334 IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce 316 IF(lwm) CALL FLUSH ( numoni ) ! flush output namelist ice 317 ENDIF 318 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 319 320 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 321 ! Trends (ua, va, tsa used as workspace) 322 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 323 IF( nstop == 0 ) THEN 324 IF( lk_trddyn ) CALL trd_dwr( kstp ) ! trends: dynamics 325 IF( lk_trdtra ) CALL trd_twr( kstp ) ! trends: active tracers 326 IF( lk_trdmld ) CALL trd_mld( kstp ) ! trends: Mixed-layer 327 IF( lk_trdvor ) CALL trd_vor( kstp ) ! trends: vorticity budget 335 IF(lwm.AND.numoni /= -1 ) & 336 & CALL FLUSH ( numoni ) ! flush output namelist ice (if exist) 328 337 ENDIF 329 338 … … 331 340 ! Coupled mode 332 341 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 333 IF( lk_cpl ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 342 !!gm why lk_oasis and not lk_cpl ???? 343 IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 334 344 ! 335 345 #if defined key_iomput 336 346 IF( kstp == nitend .OR. indic < 0 ) THEN 337 CALL iom_context_finalize( "nemo") ! needed for XIOS+AGRIF338 IF( ln_crs ) CALL iom_context_finalize( "nemo_crs" ) !347 CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 348 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 339 349 ENDIF 340 350 #endif … … 343 353 ! 344 354 END SUBROUTINE stp 345 346 !!====================================================================== 355 347 356 END MODULE step
Note: See TracChangeset
for help on using the changeset viewer.
