Changeset 5930 for trunk/NEMOGCM/NEMO/OPA_SRC/step.F90
- Timestamp:
- 2015-11-26T17:07:10+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/step.F90
r5836 r5930 28 28 !! 3.7 ! 2014-10 (G. Madec) LDF simplication 29 29 !! - ! 2014-12 (G. Madec) remove KPP scheme 30 !! - ! 2015-11 (J. Chanut) free surface simplification 30 31 !!---------------------------------------------------------------------- 31 32 … … 179 180 180 181 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 181 ! Ocean dynamics : hdiv, ssh, e3, wn 182 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 182 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 183 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 184 183 185 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 184 186 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 185 187 CALL wzv ( kstp ) ! now cross-level velocity 186 188 187 IF( lk_dynspg_ts ) THEN188 ! In case the time splitting case, update almost all momentum trends here:189 ! Note that the computation of vertical velocity above, hence "after" sea level190 ! is necessary to compute momentum advection for the rhs of barotropic loop:191 189 !!gm : why also here ???? 192 IF(ln_sto_eos )CALL sto_pts( tsn ) ! Random T/S fluctuations190 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 193 191 !!gm 194 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 195 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 196 198 IF( ln_zps .AND. .NOT. ln_isfcav) & ! Partial steps: bottom before horizontal gradient 197 199 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! of t, s, rd at the last ocean level … … 201 203 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 202 204 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 203 204 ua(:,:,:) = 0._wp ! set dynamics trends to zero 205 va(:,:,:) = 0._wp 206 IF( lk_asminc .AND. ln_asmiau .AND. & 207 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 208 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 209 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 210 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 211 CALL dyn_ldf ( kstp ) ! lateral mixing 212 #if defined key_agrif 213 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 214 #endif 215 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 216 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 217 218 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated at this stage) 219 va_sv(:,:,:) = va(:,:,:) 220 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 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 227 ! diagnostics and outputs (ua, va, tsa used as workspace) 228 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 229 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 230 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 231 IF(.NOT.ln_cpl ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 232 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 233 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 234 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 235 CALL dia_wri( kstp ) ! ocean model: outputs 236 ! 237 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 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 238 245 239 246 #if defined key_top … … 241 248 ! Passive Tracer Model 242 249 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 243 CALL trc_stp ( kstp )! time-stepping244 #endif 245 246 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 247 ! Active tracers (ua, va used as workspace)248 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 249 tsa(:,:,:,:) = 0._wp! set tracer trends to zero250 CALL trc_stp ( kstp ) ! time-stepping 251 #endif 252 253 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 254 ! Active tracers 255 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 256 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 250 257 251 258 IF( lk_asminc .AND. ln_asmiau .AND. & 252 & ln_trainc ) CALL tra_asm_inc( kstp ) ! apply tracer assimilation increment 253 CALL tra_sbc ( kstp ) ! surface boundary condition 254 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 255 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 256 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 257 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 258 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 259 CALL tra_adv ( kstp ) ! horizontal & vertical advection 260 CALL tra_ldf ( kstp ) ! lateral mixing 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 261 272 262 273 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 263 IF( ln_diaptr ) CALL dia_ptr! Poleward adv/ldf TRansports diagnostics274 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 264 275 !!gm 265 266 #if defined key_agrif 267 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 268 #endif 269 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 270 271 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 272 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 273 CALL tra_nxt( kstp ) ! tracer fields at next time step 274 !!gm : why again a call to sto_pts ??? 275 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 276 !!gm 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, gtui, gtvi, & ! Partial steps for top/bottom cells 283 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 284 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 285 ELSE ! centered hpg (eos then time stepping) 286 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 287 !!gm : why again a call to sto_pts ??? 288 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 289 !!gm 290 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 291 IF( ln_zps .AND. .NOT. ln_isfcav) & 292 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: bottom before horizontal gradient 293 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 294 IF( ln_zps .AND. ln_isfcav) & 295 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top/bottom cells 296 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 297 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 298 ENDIF 299 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 300 CALL tra_nxt( kstp ) ! tracer fields at next time step 301 ENDIF 302 303 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 304 ! Dynamics (tsa used as workspace) 305 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 306 IF( lk_dynspg_ts ) THEN 307 ! revert to previously computed momentum tendencies 308 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 309 ua(:,:,:) = ua_sv(:,:,:) 310 va(:,:,:) = va_sv(:,:,:) 311 312 CALL dyn_bfr( kstp ) ! bottom friction 313 CALL dyn_zdf( kstp ) ! vertical diffusion 314 ELSE 315 ua(:,:,:) = 0._wp ! set dynamics trends to zero 316 va(:,:,:) = 0._wp 317 318 IF( lk_asminc .AND. ln_asmiau .AND. & 319 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 320 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 321 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 322 CALL dyn_adv( kstp ) ! advection (vector or flux form) 323 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 324 CALL dyn_ldf( kstp ) ! lateral mixing 325 #if defined key_agrif 326 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 327 #endif 328 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 329 CALL dyn_bfr( kstp ) ! bottom friction 330 CALL dyn_zdf( kstp ) ! vertical diffusion 331 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 332 ENDIF 333 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 334 335 CALL ssh_swp( kstp ) ! swap of sea surface height 336 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 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 337 300 ! 338 301 339 302 !!gm : This does not only concern the dynamics ==>>> add a new title 340 303 !!gm2: why ouput restart before AGRIF update? 341 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 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 342 308 343 309 #if defined key_agrif … … 345 311 ! AGRIF 346 312 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 347 CALL Agrif_Integrate_ChildGrids( stp ) 348 349 IF ( Agrif_NbStepint().EQ.0 ) THEN 350 CALL Agrif_Update_Tra() ! Update active tracers 351 CALL Agrif_Update_Dyn() ! Update momentum 352 ENDIF 353 #endif 354 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 355 IF( lk_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( lk_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 356 323 357 324 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 358 325 ! Control 359 326 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 360 CALL stp_ctl( kstp, indic )361 IF( indic < 0 362 363 364 ENDIF 365 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 366 333 CALL iom_close( numror ) ! close input ocean restart file 367 334 IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce … … 374 341 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 375 342 !!gm why lk_oasis and not lk_cpl ???? 376 IF( lk_oasis 343 IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 377 344 ! 378 345 #if defined key_iomput
Note: See TracChangeset
for help on using the changeset viewer.