- Timestamp:
- 2015-11-30T11:47:24+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/step.F90
r5947 r5948 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 & CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF)168 & CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 161 169 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 162 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 163 IF( ln_traldf_grif ) THEN ! before slope for Griffies operator 164 CALL ldf_slp_grif( kstp ) 165 ELSE 166 CALL ldf_slp( kstp, rhd, rn2b ) ! before slope for Madec operator 170 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 171 172 IF( ln_traldf_triad ) THEN 173 CALL ldf_slp_triad( kstp ) ! before slope for triad operator 174 ELSE 175 CALL ldf_slp ( kstp, rhd, rn2b ) ! before slope for standard operator 167 176 ENDIF 168 177 ENDIF 169 #if defined key_traldf_c2d 170 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient 171 #endif 172 #if defined key_traldf_c3d && defined key_traldf_smag 173 CALL ldf_tra_smag( kstp ) ! eddy induced velocity coefficient 174 # endif 175 #if defined key_dynldf_c3d && defined key_dynldf_smag 176 CALL ldf_dyn_smag( kstp ) ! eddy induced velocity coefficient 177 # endif 178 179 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 180 ! Ocean dynamics : hdiv, rot, ssh, e3, wn 181 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 182 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_cur) 178 ! ! eddy diffusivity coeff. and/or eiv coeff. 179 IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp ) 180 181 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 182 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 183 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 184 185 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 183 186 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 184 187 CALL wzv ( kstp ) ! now cross-level velocity 185 188 186 IF( lk_dynspg_ts ) THEN 187 ! In case the time splitting case, update almost all momentum trends here: 188 ! Note that the computation of vertical velocity above, hence "after" sea level 189 ! is necessary to compute momentum advection for the rhs of barotropic loop: 190 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 191 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 192 IF( ln_zps .AND. .NOT. ln_isfcav) & 193 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 194 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 195 IF( ln_zps .AND. ln_isfcav) & 196 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 189 !!gm : why also here ???? 190 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 191 !!gm 192 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 193 194 !!jc: fs simplification 195 !!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) 196 !! but ensures reproductible results 197 !! with previous versions using split-explicit free surface 198 IF( ln_zps .AND. .NOT. ln_isfcav) & ! Partial steps: bottom before horizontal gradient 199 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! of t, s, rd at the last ocean level 200 & rhd, gru , grv ) 201 IF( ln_zps .AND. ln_isfcav) & ! Partial steps: top & bottom before horizontal gradient 202 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & 197 203 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 198 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 199 200 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 201 va(:,:,:) = 0.e0 202 IF( ln_asmiau .AND. & 203 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 204 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! subtract Neptune velocities (simplified) 205 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 206 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 207 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 208 CALL dyn_ldf ( kstp ) ! lateral mixing 209 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! add Neptune velocities (simplified) 210 #if defined key_agrif 211 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 212 #endif 213 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 214 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 215 216 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated at this stage) 217 va_sv(:,:,:) = va(:,:,:) 218 219 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 220 IF( lk_vvl ) 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 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 225 ! diagnostics and outputs (ua, va, tsa used as workspace) 226 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 227 <<<<<<< .working 204 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 205 !!jc: fs simplification 206 207 ua(:,:,:) = 0._wp ! set dynamics trends to zero 208 va(:,:,:) = 0._wp 209 210 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) & 211 CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 212 IF( lk_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends 213 #if defined key_agrif 214 IF(.NOT. Agrif_Root()) & 215 & CALL Agrif_Sponge_dyn ! momentum sponge 216 #endif 217 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 218 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 219 CALL dyn_ldf ( kstp ) ! lateral mixing 220 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 221 CALL dyn_spg ( kstp ) ! surface pressure gradient 222 223 ! With split-explicit free surface, since now transports have been updated and ssha as well 224 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 225 CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case) 226 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 227 CALL wzv ( kstp ) ! now cross-level velocity 228 ENDIF 229 230 CALL dyn_bfr ( kstp ) ! bottom friction 231 CALL dyn_zdf ( kstp ) ! vertical diffusion 232 233 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 234 ! diagnostics and outputs 235 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 228 236 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 229 237 IF( nn_diacfl == 1 ) CALL dia_cfl( kstp ) ! Courant number diagnostics 230 238 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 231 239 IF( lk_diafwb ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 232 IF( ln_diaptr ) CALL dia_ptr( kstp ) ! Poleward TRansports diagnostics233 240 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 234 241 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 235 242 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 236 243 CALL dia_wri( kstp ) ! ocean model: outputs 237 ======= 238 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 239 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 240 IF( .NOT. ln_cpl ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 241 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 242 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 243 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 244 CALL dia_wri( kstp ) ! ocean model: outputs 245 >>>>>>> .merge-right.r5518 246 ! 247 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 244 ! 245 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output 248 246 249 247 #if defined key_top … … 251 249 ! Passive Tracer Model 252 250 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 253 CALL trc_stp( kstp ) ! time-stepping 254 #endif 255 256 257 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 258 ! Active tracers (ua, va used as workspace) 259 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 260 tsa(:,:,:,:) = 0.e0 ! set tracer trends to zero 261 262 IF( ln_asmiau .AND. & 263 & ln_trainc ) CALL tra_asm_inc( kstp ) ! apply tracer assimilation increment 264 CALL tra_sbc ( kstp ) ! surface boundary condition 265 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 266 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 267 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 268 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 269 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 270 CALL tra_adv ( kstp ) ! horizontal & vertical advection 271 IF( lk_zdfkpp ) CALL tra_kpp ( kstp ) ! KPP non-local tracer fluxes 272 CALL tra_ldf ( kstp ) ! lateral mixing 273 274 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 275 276 #if defined key_agrif 277 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 278 #endif 279 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 280 281 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 282 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 283 CALL tra_nxt( kstp ) ! tracer fields at next time step 284 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 285 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 286 IF( ln_zps .AND. .NOT. ln_isfcav) & 287 & CALL zps_hde ( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps: before horizontal gradient 288 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 289 IF( ln_zps .AND. ln_isfcav) & 290 & CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps for top cell (ISF) 291 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 292 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 293 ELSE ! centered hpg (eos then time stepping) 294 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 295 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 296 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 297 IF( ln_zps .AND. .NOT. ln_isfcav) & 298 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 299 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 300 IF( ln_zps .AND. ln_isfcav) & 301 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 302 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 303 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 304 ENDIF 305 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 306 CALL tra_nxt( kstp ) ! tracer fields at next time step 307 ENDIF 308 309 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 310 ! Dynamics (tsa used as workspace) 311 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 312 IF( lk_dynspg_ts ) THEN 313 ! revert to previously computed momentum tendencies 314 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 315 ua(:,:,:) = ua_sv(:,:,:) 316 va(:,:,:) = va_sv(:,:,:) 317 ! Revert now divergence and rotational to previously computed ones 318 !(needed because of the time swap in div_cur, at the beginning of each time step) 319 hdivn(:,:,:) = hdivb(:,:,:) 320 rotn(:,:,:) = rotb(:,:,:) 321 322 CALL dyn_bfr( kstp ) ! bottom friction 323 CALL dyn_zdf( kstp ) ! vertical diffusion 324 ELSE 325 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 326 va(:,:,:) = 0.e0 327 328 IF( ln_asmiau .AND. & 329 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 330 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 331 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 332 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 333 CALL dyn_adv( kstp ) ! advection (vector or flux form) 334 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 335 CALL dyn_ldf( kstp ) ! lateral mixing 336 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 337 #if defined key_agrif 338 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 339 #endif 340 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 341 CALL dyn_bfr( kstp ) ! bottom friction 342 CALL dyn_zdf( kstp ) ! vertical diffusion 343 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 344 ENDIF 345 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 346 347 CALL ssh_swp( kstp ) ! swap of sea surface height 348 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 349 350 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 351 IF( lk_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 352 353 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 354 ! Control and restarts 355 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 356 CALL stp_ctl( kstp, indic ) 357 IF( indic < 0 ) THEN 358 CALL ctl_stop( 'step: indic < 0' ) 359 CALL dia_wri_state( 'output.abort', kstp ) 360 ENDIF 361 IF( kstp == nit000 ) THEN 251 CALL trc_stp ( kstp ) ! time-stepping 252 #endif 253 254 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 255 ! Active tracers 256 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 257 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 258 259 IF( lk_asminc .AND. ln_asmiau .AND. & 260 & ln_trainc ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment 261 CALL tra_sbc ( kstp ) ! surface boundary condition 262 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 263 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 264 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 265 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 266 IF( lk_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends 267 #if defined key_agrif 268 IF(.NOT. Agrif_Root()) & 269 & CALL Agrif_Sponge_tra ! tracers sponge 270 #endif 271 CALL tra_adv ( kstp ) ! horizontal & vertical advection 272 CALL tra_ldf ( kstp ) ! lateral mixing 273 274 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 275 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 276 !!gm 277 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 278 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 279 280 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 281 ! Set boundary conditions and Swap 282 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 283 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 284 !! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 285 !! If so: 286 !! (i) no need to call agrif update at initialization time 287 !! (ii) no need to update "before" fields 288 !! 289 !! Apart from creating new tra_swp/dyn_swp routines, this however: 290 !! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 291 !! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 292 !! e.g. a shift of the feedback interface inside child domain. 293 !! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 294 !! place. 295 !! 296 !!jc2: dynnxt must be the latest call. fse3t_b are indeed updated in that routine 297 CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap 298 CALL dyn_nxt ( kstp ) ! finalize (bcs) velocities at next time step and swap 299 CALL ssh_swp ( kstp ) ! swap of sea surface height 300 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 301 ! 302 303 !!gm : This does not only concern the dynamics ==>>> add a new title 304 !!gm2: why ouput restart before AGRIF update? 305 !! 306 !!jc: That would be better, but see comment above 307 !! 308 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 309 310 #if defined key_agrif 311 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 312 ! AGRIF 313 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 314 CALL Agrif_Integrate_ChildGrids( stp ) 315 316 IF ( Agrif_NbStepint().EQ.0 ) THEN ! AGRIF Update 317 !!jc in fact update i useless at last time step, but do it for global diagnostics 318 CALL Agrif_Update_Tra() ! Update active tracers 319 CALL Agrif_Update_Dyn() ! Update momentum 320 ENDIF 321 #endif 322 IF( ln_diahsb ) CALL dia_hsb ( kstp ) ! - ML - global conservation diagnostics 323 IF( lk_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 324 325 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 326 ! Control 327 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 328 CALL stp_ctl ( kstp, indic ) 329 IF( indic < 0 ) THEN 330 CALL ctl_stop( 'step: indic < 0' ) 331 CALL dia_wri_state( 'output.abort', kstp ) 332 ENDIF 333 IF( kstp == nit000 ) THEN 362 334 CALL iom_close( numror ) ! close input ocean restart file 363 335 IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce 364 IF( lwm.AND.numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice365 ENDIF366 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file336 IF(lwm.AND.numoni /= -1 ) & 337 & CALL FLUSH ( numoni ) ! flush output namelist ice (if exist) 338 ENDIF 367 339 368 340 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 369 341 ! Coupled mode 370 342 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 371 IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 343 !!gm why lk_oasis and not lk_cpl ???? 344 IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 372 345 ! 373 346 #if defined key_iomput
Note: See TracChangeset
for help on using the changeset viewer.
