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