- Timestamp:
- 2015-12-02T11:52:05+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/step.F90
r5682 r5974 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) … … 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 ( 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 86 IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE. 87 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. 88 88 # if defined key_iomput 89 89 IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context ) 90 90 # endif 91 91 #endif 92 indic = 0 ! reset to no error condition 93 IF( kstp == nit000 ) THEN 94 ! must be done after nemo_init for AGRIF+XIOS+OASIS 95 CALL iom_init( cxios_context ) ! iom_put initialization 96 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! initialize context for coarse grid 97 ENDIF 98 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 99 101 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) 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 data, open boundaries, surface boundary condition (including sea-ice) 105 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 106 IF( lk_tide ) CALL sbc_tide( kstp ) 107 IF( lk_bdy ) THEN 108 IF( ln_apr_dyn) CALL sbc_apr( kstp ) ! bdy_dta needs ssh_ib 109 CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 110 ENDIF 111 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 112 ! 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 113 113 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 114 114 ! Update stochastic parameters and random T/S fluctuations 115 115 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 116 CALL sto_par( kstp ) ! Stochastic parameters116 CALL sto_par( kstp ) ! Stochastic parameters 117 117 118 118 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 124 124 CALL bn2 ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 125 125 CALL bn2 ( tsn, rab_n, rn2 ) ! now Brunt-Vaisala frequency 126 126 127 ! 127 128 ! VERTICAL PHYSICS … … 131 132 IF( lk_zdftke ) CALL zdf_tke( kstp ) ! TKE closure scheme for Kz 132 133 IF( lk_zdfgls ) CALL zdf_gls( kstp ) ! GLS closure scheme for Kz 133 IF( lk_zdfkpp ) CALL zdf_kpp( kstp ) ! KPP closure scheme for Kz134 134 IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 135 135 avt (:,:,:) = rn_avt0 * wmask (:,:,:) … … 137 137 avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 138 138 ENDIF 139 139 140 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 140 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 141 142 ENDIF 142 143 IF( ln_zdfevd ) CALL zdf_evd( kstp ) ! enhanced vertical eddy diffusivity … … 144 145 IF( lk_zdftmx ) CALL zdf_tmx( kstp ) ! tidal vertical mixing 145 146 146 IF( lk_zdfddm .AND. .NOT. lk_zdfkpp ) & 147 & CALL zdf_ddm( kstp ) ! double diffusive mixing 147 IF( lk_zdfddm ) CALL zdf_ddm( kstp ) ! double diffusive mixing 148 148 149 149 CALL zdf_mxl( kstp ) ! mixed layer depth … … 155 155 ! LATERAL PHYSICS 156 156 ! 157 IF( lk_ldfslp ) THEN ! slope of lateral mixing 157 IF( l_ldfslp ) THEN ! slope of lateral mixing 158 !!gm : why this here ???? 158 159 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 160 !!gm 159 161 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 162 160 163 IF( ln_zps .AND. .NOT. ln_isfcav) & 161 164 & CALL zps_hde ( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 162 165 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 166 163 167 IF( ln_zps .AND. ln_isfcav) & 164 & 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) 165 169 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 166 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the first ocean level 167 IF( ln_traldf_grif ) THEN ! before slope for Griffies operator 168 CALL ldf_slp_grif( kstp ) 169 ELSE 170 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 171 176 ENDIF 172 177 ENDIF 173 #if defined key_traldf_c2d 174 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient 175 #endif 176 #if defined key_traldf_c3d && defined key_traldf_smag 177 CALL ldf_tra_smag( kstp ) ! eddy induced velocity coefficient 178 # endif 179 #if defined key_dynldf_c3d && defined key_dynldf_smag 180 CALL ldf_dyn_smag( kstp ) ! eddy induced velocity coefficient 181 # endif 182 183 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 184 ! Ocean dynamics : hdiv, rot, ssh, e3, wn 185 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 186 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) 187 186 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 188 187 CALL wzv ( kstp ) ! now cross-level velocity 189 188 190 IF( lk_dynspg_ts ) THEN 191 ! In case the time splitting case, update almost all momentum trends here: 192 ! Note that the computation of vertical velocity above, hence "after" sea level 193 ! is necessary to compute momentum advection for the rhs of barotropic loop: 194 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 195 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 196 IF( ln_zps .AND. .NOT. ln_isfcav) & 197 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 198 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 199 IF( ln_zps .AND. ln_isfcav) & 200 & 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, & 201 203 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 202 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 203 204 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 205 va(:,:,:) = 0.e0 206 IF( lk_asminc .AND. ln_asmiau .AND. & 207 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 208 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! subtract Neptune velocities (simplified) 209 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 210 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 211 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 212 CALL dyn_ldf ( kstp ) ! lateral mixing 213 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! add Neptune velocities (simplified) 214 #if defined key_agrif 215 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 216 #endif 217 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 218 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 219 220 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated at this stage) 221 va_sv(:,:,:) = va(:,:,:) 222 223 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 224 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 225 CALL wzv ( kstp ) ! now cross-level velocity 226 ENDIF 227 228 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 229 ! diagnostics and outputs (ua, va, tsa used as workspace) 230 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 231 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 232 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 233 IF( .NOT. ln_cpl ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 234 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 235 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 236 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 237 CALL dia_wri( kstp ) ! ocean model: outputs 238 ! 239 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 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 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 236 IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats 237 IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 238 IF(.NOT.ln_cpl ) CALL dia_fwb ( kstp ) ! Fresh water budget diagnostics 239 IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports 240 IF( lk_diaar5 ) CALL dia_ar5 ( kstp ) ! ar5 diag 241 IF( lk_diaharm ) CALL dia_harm ( kstp ) ! Tidal harmonic analysis 242 CALL dia_wri ( kstp ) ! ocean model: outputs 243 ! 244 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output 240 245 241 246 #if defined key_top … … 243 248 ! Passive Tracer Model 244 249 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 245 CALL trc_stp( kstp ) ! time-stepping 246 #endif 247 248 249 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 250 ! Active tracers (ua, va used as workspace) 251 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 252 tsa(:,:,:,:) = 0.e0 ! set tracer trends to zero 250 CALL trc_stp ( kstp ) ! time-stepping 251 #endif 252 253 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 254 ! Active tracers 255 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 256 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 253 257 254 258 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 CALL tra_adv ( kstp ) ! horizontal & vertical advection 263 IF( lk_zdfkpp ) CALL tra_kpp ( kstp ) ! KPP non-local tracer fluxes 264 CALL tra_ldf ( kstp ) ! lateral mixing 265 266 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 267 268 #if defined key_agrif 269 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 270 #endif 271 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 272 273 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 274 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 275 CALL tra_nxt( kstp ) ! tracer fields at next time step 276 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 277 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 278 IF( ln_zps .AND. .NOT. ln_isfcav) & 279 & CALL zps_hde ( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps: before horizontal gradient 280 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 281 IF( ln_zps .AND. ln_isfcav) & 282 & CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps for top cell (ISF) 283 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 284 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 285 ELSE ! centered hpg (eos then time stepping) 286 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 287 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 288 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 289 IF( ln_zps .AND. .NOT. ln_isfcav) & 290 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: before horizontal gradient 291 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 292 IF( ln_zps .AND. ln_isfcav) & 293 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps for top cell (ISF) 294 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 295 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 296 ENDIF 297 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 298 CALL tra_nxt( kstp ) ! tracer fields at next time step 299 ENDIF 300 301 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 302 ! Dynamics (tsa used as workspace) 303 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 304 IF( lk_dynspg_ts ) THEN 305 ! revert to previously computed momentum tendencies 306 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 307 ua(:,:,:) = ua_sv(:,:,:) 308 va(:,:,:) = va_sv(:,:,:) 309 ! Revert now divergence and rotational to previously computed ones 310 !(needed because of the time swap in div_cur, at the beginning of each time step) 311 hdivn(:,:,:) = hdivb(:,:,:) 312 rotn(:,:,:) = rotb(:,:,:) 313 314 CALL dyn_bfr( kstp ) ! bottom friction 315 CALL dyn_zdf( kstp ) ! vertical diffusion 316 ELSE 317 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 318 va(:,:,:) = 0.e0 319 320 IF( lk_asminc .AND. ln_asmiau .AND. & 321 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 322 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 323 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 324 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 325 CALL dyn_adv( kstp ) ! advection (vector or flux form) 326 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 327 CALL dyn_ldf( kstp ) ! lateral mixing 328 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 329 #if defined key_agrif 330 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 331 #endif 332 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 333 CALL dyn_bfr( kstp ) ! bottom friction 334 CALL dyn_zdf( kstp ) ! vertical diffusion 335 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 336 ENDIF 337 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 338 339 CALL ssh_swp( kstp ) ! swap of sea surface height 340 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 341 ! 342 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 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. fse3t_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( lk_vvl ) 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 343 308 344 309 #if defined key_agrif … … 346 311 ! AGRIF 347 312 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 348 CALL Agrif_Integrate_ChildGrids( stp ) 349 350 IF ( Agrif_NbStepint().EQ.0 ) THEN 351 CALL Agrif_Update_Tra() ! Update active tracers 352 CALL Agrif_Update_Dyn() ! Update momentum 353 ENDIF 354 #endif 355 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 356 IF( ln_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 313 CALL Agrif_Integrate_ChildGrids( stp ) 314 315 IF ( Agrif_NbStepint().EQ.0 ) THEN ! AGRIF Update 316 !!jc in fact update i 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) 357 323 358 324 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 359 325 ! Control 360 326 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 361 CALL stp_ctl( kstp, indic )362 IF( indic < 0 363 364 365 ENDIF 366 IF( kstp == nit000 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 367 333 CALL iom_close( numror ) ! close input ocean restart file 368 334 IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce 369 IF( lwm.AND.numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice 335 IF(lwm.AND.numoni /= -1 ) & 336 & CALL FLUSH ( numoni ) ! flush output namelist ice (if exist) 370 337 ENDIF 371 338 … … 373 340 ! Coupled mode 374 341 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 375 IF( lk_oasis ) 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 376 344 ! 377 345 #if defined key_iomput … … 383 351 ! 384 352 IF( nn_timing == 1 .AND. kstp == nit000 ) CALL timing_reset 385 !386 353 ! 387 354 END SUBROUTINE stp
Note: See TracChangeset
for help on using the changeset viewer.