New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
step.F90 in branches/2012/dev_r3452_UKMO9_RESTART/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/2012/dev_r3452_UKMO9_RESTART/NEMOGCM/NEMO/OPA_SRC/step.F90 @ 3594

Last change on this file since 3594 was 3594, checked in by rfurner, 11 years ago

code not tested through SETTEE, builds and runs, but has not been thoroughly tested, so will not be included in 2012 merge, however submitted back to keep record of work done for 2013 developments
  • Property svn:keywords set to Id
File size: 15.9 KB
Line 
1MODULE step
2   !!======================================================================
3   !!                       ***  MODULE step  ***
4   !! Time-stepping    : manager of the ocean, tracer and ice time stepping
5   !!======================================================================
6   !! History :  OPA  !  1991-03  (G. Madec)  Original code
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   !!            8.0  !  1997-06  (G. Madec)  new architecture of call
12   !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
13   !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit
14   !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
15   !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
16   !!             -   !  2004-08  (C. Talandier) New trends organization
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   !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
22   !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl
23   !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
24   !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
25   !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal
26   !!----------------------------------------------------------------------
27
28   !!----------------------------------------------------------------------
29   !!   stp             : OPA system time-stepping
30   !!----------------------------------------------------------------------
31   USE step_oce         ! time stepping definition modules
32#if defined key_top
33   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
34#endif
35#if defined key_agrif
36   USE agrif_opa_sponge ! Momemtum and tracers sponges
37#endif
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   stp   ! called by opa.F90
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "zdfddm_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54#if defined key_agrif
55   SUBROUTINE stp( )
56      INTEGER             ::   kstp   ! ocean time-step index
57#else
58   SUBROUTINE stp( kstp )
59      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
60#endif
61      !!----------------------------------------------------------------------
62      !!                     ***  ROUTINE stp  ***
63      !!                     
64      !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
65      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
66      !!              - Tme stepping  of TRC (passive tracer eqs.)
67      !!
68      !! ** Method  : -1- Update forcings and data 
69      !!              -2- Update ocean physics
70      !!              -3- Compute the t and s trends
71      !!              -4- Update t and s
72      !!              -5- Compute the momentum trends
73      !!              -6- Update the horizontal velocity
74      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
75      !!              -8- Outputs and diagnostics
76      !!----------------------------------------------------------------------
77      INTEGER ::   jk       ! dummy loop indice
78      INTEGER ::   indic    ! error indicator if < 0
79      !! ---------------------------------------------------------------------
80
81#if defined key_agrif
82      kstp = nit000 + Agrif_Nb_Step()
83!      IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
84!      IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
85# if defined key_iomput
86      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap
87# endif   
88#endif   
89                             indic = 0                ! reset to no error condition
90      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
91                             CALL iom_setkt( kstp )   ! say to iom that we are at time step kstp
92
93      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
94      ! Update data, open boundaries, surface boundary condition (including sea-ice)
95      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
96                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
97      IF( lk_tide    )   CALL sbc_tide( kstp )
98      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
99      IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
100      IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries
101
102      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
103      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
104      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
105                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
106
107      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
108      ! Ocean physics update                (ua, va, tsa used as workspace)
109      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
110                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
111                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
112      !
113      !  VERTICAL PHYSICS   
114                         CALL zdf_bfr( kstp )         ! bottom friction
115                         
116      !                                               ! Vertical eddy viscosity and diffusivity coefficients
117      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
118      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
119      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
120      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
121      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
122         avt (:,:,:) = rn_avt0 * tmask(:,:,:)
123         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
124         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
125      ENDIF
126      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
127         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
128      ENDIF
129      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
130
131      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
132
133      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
134         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
135         
136                         CALL zdf_mxl( kstp )         ! mixed layer depth
137
138                                                      ! write TKE or GLS information in the restart file
139      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
140      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
141      !
142      !  LATERAL  PHYSICS
143      !
144      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
145                         CALL eos( tsb, rhd )                ! before in situ density
146         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
147            &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
148         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
149                         CALL ldf_slp_grif( kstp )
150         ELSE
151                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
152         ENDIF
153      ENDIF
154#if defined key_traldf_c2d
155      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
156#endif
157
158      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
159      ! diagnostics and outputs             (ua, va, tsa used as workspace)
160      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
161      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
162      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
163      IF( lk_diafwb  )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
164      IF( ln_diaptr  )   CALL dia_ptr( kstp )         ! Poleward TRansports diagnostics
165      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
166      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
167      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
168                         CALL dia_wri( kstp )         ! ocean model: outputs
169
170#if defined key_top
171      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
172      ! Passive Tracer Model
173      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
174                         CALL trc_stp( kstp )         ! time-stepping
175#endif
176
177      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
178      ! Active tracers                              (ua, va used as workspace)
179      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
180                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
181
182      IF(  ln_asmiau .AND. &
183         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
184                             CALL tra_sbc    ( kstp )       ! surface boundary condition
185      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
186      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
187      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
188      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
189                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection
190      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
191                             CALL tra_ldf    ( kstp )       ! lateral mixing
192#if defined key_agrif
193      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
194#endif
195                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
196
197      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
198         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
199                             CALL tra_nxt( kstp )                ! tracer fields at next time step
200                             CALL eos    ( tsa, rhd, rhop )      ! Time-filtered in situ density for hpg computation
201         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative
202            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
203         
204      ELSE                                                  ! centered hpg  (eos then time stepping)
205                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation
206         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
207            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
208         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
209                             CALL tra_nxt( kstp )                ! tracer fields at next time step
210      ENDIF 
211
212      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
213      ! Dynamics                                    (tsa used as workspace)
214      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
215                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero
216                               va(:,:,:) = 0.e0
217
218      IF(  ln_asmiau .AND. &
219         & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
220      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
221                               CALL dyn_adv( kstp )         ! advection (vector or flux form)
222                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis
223                               CALL dyn_ldf( kstp )         ! lateral mixing
224      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
225#if defined key_agrif
226      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
227#endif
228                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
229                               CALL dyn_bfr( kstp )         ! bottom friction   
230                               CALL dyn_zdf( kstp )         ! vertical diffusion
231                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient
232                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step
233
234                               CALL ssh_nxt( kstp )         ! sea surface height at next time step
235
236      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
237      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
238
239      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
240      ! Control and restarts
241      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
242                               CALL stp_ctl( kstp, indic )
243      IF( indic < 0        )   THEN
244                               CALL ctl_stop( 'step: indic < 0' )
245                               CALL dia_wri_state( 'output.abort', kstp )
246      ENDIF
247      IF( kstp == nit000   )      CALL iom_close( numror )     ! close input  ocean restart file
248      IF( lrst_oce         )      CALL rst_write    ( kstp )   ! write output ocean restart file
249      IF( lrst_oce .and. lk_obc ) CALL obc_rst_write( kstp )   ! write open boundary restart file
250
251      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
252      ! Trends                              (ua, va, tsa used as workspace)
253      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
254      IF( nstop == 0 ) THEN                         
255         IF( lk_trddyn     )   CALL trd_dwr( kstp )         ! trends: dynamics
256         IF( lk_trdtra     )   CALL trd_twr( kstp )         ! trends: active tracers
257         IF( lk_trdmld     )   CALL trd_mld( kstp )         ! trends: Mixed-layer
258         IF( lk_trdvor     )   CALL trd_vor( kstp )         ! trends: vorticity budget
259      ENDIF
260
261      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
262      ! Coupled mode
263      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
264      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
265      !
266      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
267      !
268   END SUBROUTINE stp
269
270   !!======================================================================
271END MODULE step
Note: See TracBrowser for help on using the repository browser.