Changeset 1492 for trunk/NEMO/OPA_SRC/step.F90
- Timestamp:
- 2009-07-17T11:48:07+02:00 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/step.F90
r1482 r1492 5 5 !!====================================================================== 6 6 !! History : OPA ! 1991-03 (G. Madec) Original code 7 !! 8 !! 9 !! 10 !! 7 !! - ! 1991-11 (G. Madec) 8 !! - ! 1992-06 (M. Imbard) add a first output record 9 !! - ! 1996-04 (G. Madec) introduction of dynspg 10 !! - ! 1996-04 (M.A. Foujols) introduction of passive tracer 11 11 !! 8.0 ! 1997-06 (G. Madec) new architecture of call 12 12 !! 8.2 ! 1997-06 (G. Madec, M. Imbard, G. Roullet) free surface … … 15 15 !! NEMO 1.0 ! 2002-06 (G. Madec) free form, suppress macro-tasking 16 16 !! - ! 2004-08 (C. Talandier) New trends organization 17 !! 18 !! 19 !! 20 !! 17 !! - ! 2005-01 (C. Ethe) Add the KPP closure scheme 18 !! - ! 2005-11 (G. Madec) Reorganisation of tra and dyn calls 19 !! - ! 2006-01 (L. Debreu, C. Mazauric) Agrif implementation 20 !! - ! 2006-07 (S. Masson) restart using iom 21 21 !! 3.2 ! 2009-02 (G. Madec, R. Benshila) reintroduicing z*-coordinate 22 !! - ! 2009-06 (S. Masson, G. Madec) TKE2 restart compatible with key_cpl 22 23 !!---------------------------------------------------------------------- 23 24 … … 123 124 PRIVATE 124 125 125 PUBLIC stp! called by opa.F90126 PUBLIC stp ! called by opa.F90 126 127 127 128 !! * Substitutions … … 129 130 # include "zdfddm_substitute.h90" 130 131 !!---------------------------------------------------------------------- 131 !! OPA 9.0 , LOCEAN-IPSL (2005)132 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 132 133 !! $Id$ 133 134 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 179 180 ! Update data, open boundaries, surface boundary condition (including sea-ice) 180 181 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 181 182 182 IF( lk_dtatem ) CALL dta_tem( kstp ) ! update 3D temperature data 183 183 IF( lk_dtasal ) CALL dta_sal( kstp ) ! update 3D salinity data 184 185 184 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 186 187 185 IF( lk_obc ) CALL obc_dta( kstp ) ! update dynamic and tracer data at open boundaries 188 186 IF( lk_obc ) CALL obc_rad( kstp ) ! compute phase velocities at open boundaries 189 190 187 IF( lk_bdy ) CALL bdy_dta( kstp ) ! update dynamic and tracer data at unstructured open boundary 191 188 192 189 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 193 ! Ocean dynamics : ssh, wn, hdiv, rot ! 194 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 195 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity 196 IF( n_cla == 1 ) CALL div_cla( kstp ) ! Cross Land Advection (Update Hor. divergence) 197 CALL ssh_wzv( kstp ) ! after ssh & vertical velocity 198 199 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 200 ! Ocean physics update 201 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 202 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 203 !----------------------------------------------------------------------- 204 205 CALL bn2( tb, sb, rn2b ) ! before Brunt-Vaisala frequency 206 CALL bn2( tn, sn, rn2 ) ! now Brunt-Vaisala frequency 207 208 !----------------------------------------------------------------------- 209 ! VERTICAL PHYSICS 210 !----------------------------------------------------------------------- 211 ! ! Vertical eddy viscosity and diffusivity coefficients 212 IF( lk_zdfric ) CALL zdf_ric( kstp ) ! Richardson number dependent Kz 213 214 IF( lk_zdftke ) CALL zdf_tke ( kstp ) ! TKE closure scheme for Kz 215 IF( lk_zdftke2) CALL zdf_tke2( kstp ) ! TKE2 closure scheme for Kz 216 217 IF( lk_zdfkpp ) CALL zdf_kpp( kstp ) ! KPP closure scheme for Kz 218 219 IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 190 ! Ocean dynamics : ssh, wn, hdiv, rot ! 191 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 192 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity 193 IF( n_cla == 1 ) CALL div_cla( kstp ) ! Cross Land Advection (Update Hor. divergence) 194 CALL ssh_wzv( kstp ) ! after ssh & vertical velocity 195 196 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 197 ! Ocean physics update (ua, va, ta, sa used as workspace) 198 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 199 CALL bn2( tb, sb, rn2b ) ! before Brunt-Vaisala frequency 200 CALL bn2( tn, sn, rn2 ) ! now Brunt-Vaisala frequency 201 ! 202 ! VERTICAL PHYSICS 203 ! ! Vertical eddy viscosity and diffusivity coefficients 204 IF( lk_zdfric ) CALL zdf_ric( kstp ) ! Richardson number dependent Kz 205 IF( lk_zdftke ) CALL zdf_tke ( kstp ) ! TKE closure scheme for Kz 206 IF( lk_zdftke2 ) CALL zdf_tke2( kstp ) ! TKE2 closure scheme for Kz 207 IF( lk_zdfkpp ) CALL zdf_kpp( kstp ) ! KPP closure scheme for Kz 208 IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 220 209 avt (:,:,:) = avt0 * tmask(:,:,:) 221 210 avmu(:,:,:) = avm0 * umask(:,:,:) 222 211 avmv(:,:,:) = avm0 * vmask(:,:,:) 223 212 ENDIF 224 225 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 213 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 226 214 DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk) ; END DO 227 215 ENDIF 228 229 IF( ln_zdfevd ) CALL zdf_evd( kstp ) ! enhanced vertical eddy diffusivity 230 231 IF( lk_zdftmx ) CALL zdf_tmx( kstp ) ! tidal vertical mixing 216 IF( ln_zdfevd ) CALL zdf_evd( kstp ) ! enhanced vertical eddy diffusivity 217 218 IF( lk_zdftmx ) CALL zdf_tmx( kstp ) ! tidal vertical mixing 232 219 233 220 IF( lk_zdfddm .AND. .NOT. lk_zdfkpp ) & 234 & CALL zdf_ddm( kstp ) ! double diffusive mixing 235 236 237 CALL zdf_bfr( kstp ) ! bottom friction 238 239 CALL zdf_mxl( kstp ) ! mixed layer depth 240 241 IF( lrst_oce .AND. lk_zdftke2 ) & ! write tke2 information in the restart file 242 & CALL tke2_rst( kstp, 'WRITE' ) 243 244 !----------------------------------------------------------------------- 245 ! LATERAL PHYSICS 246 !----------------------------------------------------------------------- 247 IF( lk_ldfslp ) THEN 248 CALL eos( tb, sb, rhd, rhop ) ! before in situ density 249 IF( ln_zps ) CALL zps_hde( kstp, tb, sb, rhd, & ! Partial steps: before horizontal gradient 250 & gtu, gsu, gru, & ! of t, s, rd at the last ocean level 251 & gtv, gsv, grv ) 252 CALL ldf_slp( kstp, rhd, rn2b ) ! before slope of the lateral mixing 221 & CALL zdf_ddm( kstp ) ! double diffusive mixing 222 223 CALL zdf_bfr( kstp ) ! bottom friction 224 225 CALL zdf_mxl( kstp ) ! mixed layer depth 226 227 ! write tke2 information in the restart file 228 IF( lrst_oce .AND. lk_zdftke2 ) CALL tke2_rst( kstp, 'WRITE' ) 229 ! 230 ! LATERAL PHYSICS 231 ! 232 IF( lk_ldfslp ) THEN ! slope of lateral mixing 233 CALL eos( tb, sb, rhd, rhop ) ! before in situ density 234 IF( ln_zps ) CALL zps_hde( kstp, tb, sb, rhd, & ! Partial steps: before horizontal gradient 235 & gtu, gsu, gru, & ! of t, s, rd at the last ocean level 236 & gtv, gsv, grv ) 237 CALL ldf_slp( kstp, rhd, rn2b ) ! before slope of the lateral mixing 253 238 ENDIF 254 239 #if defined key_traldf_c2d 255 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) 240 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient 256 241 # endif 257 242 258 243 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 259 ! diagnostics and outputs 260 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 261 CALL dia_wri( kstp, indic ) 262 IF( lk_floats ) CALL flo_stp( kstp ) 263 IF( lk_diaspr ) CALL dia_spr( kstp ) 264 IF( lk_diahth ) CALL dia_hth( kstp ) 265 IF( lk_diagap ) CALL dia_gap( kstp ) 266 IF( lk_diahdy ) CALL dia_hdy( kstp ) 267 IF( lk_diafwb ) CALL dia_fwb( kstp ) 268 IF( ln_diaptr ) CALL dia_ptr( kstp ) 244 ! diagnostics and outputs (ua, va, ta, sa used as workspace) 245 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 246 CALL dia_wri( kstp, indic ) ! ocean model: outputs 247 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 248 IF( lk_diaspr ) CALL dia_spr( kstp ) ! Surface pressure diagnostics 249 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 250 IF( lk_diagap ) CALL dia_gap( kstp ) ! basin averaged diagnostics 251 IF( lk_diahdy ) CALL dia_hdy( kstp ) ! dynamical heigh diagnostics 252 IF( lk_diafwb ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 253 IF( ln_diaptr ) CALL dia_ptr( kstp ) ! Poleward TRansports diagnostics 269 254 270 255 #if defined key_top … … 272 257 ! Passive Tracer Model 273 258 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 274 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 275 !----------------------------------------------------------------------- 276 CALL trc_stp( kstp ) ! time-stepping 277 #endif 278 279 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 280 ! Active tracers 281 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 282 ! N.B. ua, va arrays are used as workspace in this section 283 !----------------------------------------------------------------------- 259 CALL trc_stp( kstp ) ! time-stepping 260 #endif 261 262 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 263 ! Active tracers (ua, va used as workspace) 264 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 284 265 ta(:,:,:) = 0.e0 ! set tracer trends to zero 285 266 sa(:,:,:) = 0.e0 … … 300 281 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 301 282 302 IF( .NOT. ln_dynhpg_imp ) THEN ! centered hpg (default case) 303 CALL eos( tn, sn, rhd, rhop ) ! now (swap=before) in situ density for dynhpg module 304 IF( ln_zps ) CALL zps_hde( kstp, tn, sn, rhd, & ! Partial steps: now horizontal gradient 305 & gtu, gsu, gru, & ! of t, s, rd at the bottom ocean level 306 & gtv, gsv, grv ) 283 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 284 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 285 CALL tra_nxt ( kstp ) ! tracer fields at next time step 286 CALL eos( ta, sa, rhd, rhop ) ! Time-filtered in situ density for hpg computation 287 IF( ln_zps ) CALL zps_hde( kstp, ta, sa, rhd, & ! Partial steps: time filtered hor. derivative 288 & gtu, gsu, gru, & ! of t, s, rd at the bottom ocean level 289 & gtv, gsv, grv ) 290 291 ELSE ! centered hpg (eos then time stepping) 292 CALL eos( tn, sn, rhd, rhop ) ! now in situ density for hpg computation 293 IF( ln_zps ) CALL zps_hde( kstp, tn, sn, rhd, & ! Partial steps: now horizontal derivative 294 & gtu, gsu, gru, & ! of t, s, rd at the bottom ocean level 295 & gtv, gsv, grv ) 296 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 297 CALL tra_nxt ( kstp ) ! tracer fields at next time step 307 298 ENDIF 308 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 309 CALL tra_nxt ( kstp ) ! tracer fields at next time step 310 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg 311 CALL eos( ta, sa, rhd, rhop ) ! Time-filtered in situ density used in dynhpg module 312 IF( ln_zps ) CALL zps_hde( kstp, ta, sa, rhd, & ! Partial steps: time filtered hor. gradient 313 & gtu, gsu, gru, & ! of t, s, rd at the bottom ocean level 314 & gtv, gsv, grv ) 315 ENDIF 316 317 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 318 ! Dynamics 319 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 320 ! N.B. ta, sa arrays are used as workspace in this section 321 !----------------------------------------------------------------------- 322 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 299 300 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 301 ! Dynamics (ta, sa used as workspace) 302 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 303 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 323 304 va(:,:,:) = 0.e0 324 305 325 CALL dyn_adv( kstp ) 326 CALL dyn_vor( kstp ) 327 CALL dyn_ldf( kstp ) 328 #if defined key_agrif 329 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn 330 #endif 331 CALL dyn_hpg( kstp ) 332 CALL dyn_zdf( kstp ) 306 CALL dyn_adv( kstp ) ! advection (vector or flux form) 307 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 308 CALL dyn_ldf( kstp ) ! lateral mixing 309 #if defined key_agrif 310 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 311 #endif 312 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 313 CALL dyn_zdf( kstp ) ! vertical diffusion 333 314 IF( lk_dynspg_rl ) THEN 334 IF( lk_obc ) CALL obc_spg( kstp ) 315 IF( lk_obc ) CALL obc_spg( kstp ) ! surface pressure gradient at open boundaries 335 316 ENDIF 336 317 indic=0 337 CALL dyn_spg( kstp, indic ) 338 CALL dyn_nxt( kstp ) 339 340 CALL ssh_nxt( kstp ) 318 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 319 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 320 321 CALL ssh_nxt( kstp ) ! sea surface height at next time step 341 322 342 323 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 343 324 ! Control and restarts 344 325 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 345 CALL stp_ctl( kstp, indic ) 346 IF( indic < 0 ) THEN 347 CALL ctl_stop( 'step: indic < 0' ) 348 CALL dia_wri_state( 'output.abort', kstp ) 349 ENDIF 350 IF( kstp == nit000 ) CALL iom_close( numror ) ! close input ocean restart file 351 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 352 IF( lk_obc ) CALL obc_rst_write( kstp ) ! write open boundary restart file 353 354 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 355 ! Trends 356 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 357 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 358 !----------------------------------------------------------------------- 359 360 IF( nstop == 0 ) THEN ! Diagnostics: 361 IF( lk_trddyn ) CALL trd_dwr( kstp ) ! trends: dynamics 362 IF( lk_trdtra ) CALL trd_twr( kstp ) ! trends: active tracers 363 IF( lk_trdmld ) CALL trd_mld( kstp ) ! trends: Mixed-layer 364 IF( lk_trdvor ) CALL trd_vor( kstp ) ! trends: vorticity budget 326 CALL stp_ctl( kstp, indic ) 327 IF( indic < 0 ) THEN 328 CALL ctl_stop( 'step: indic < 0' ) 329 CALL dia_wri_state( 'output.abort', kstp ) 330 ENDIF 331 IF( kstp == nit000 ) CALL iom_close( numror ) ! close input ocean restart file 332 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 333 IF( lk_obc ) CALL obc_rst_write( kstp ) ! write open boundary restart file 334 335 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 336 ! Trends (ua, va, ta, sa used as workspace) 337 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 338 IF( nstop == 0 ) THEN 339 IF( lk_trddyn ) CALL trd_dwr( kstp ) ! trends: dynamics 340 IF( lk_trdtra ) CALL trd_twr( kstp ) ! trends: active tracers 341 IF( lk_trdmld ) CALL trd_mld( kstp ) ! trends: Mixed-layer 342 IF( lk_trdvor ) CALL trd_vor( kstp ) ! trends: vorticity budget 365 343 ENDIF 366 344 … … 368 346 ! Coupled mode 369 347 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 370 371 IF( lk_cpl ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 348 IF( lk_cpl ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 372 349 ! 373 350 !
Note: See TracChangeset
for help on using the changeset viewer.