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/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/step.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

  1. Initial restructuring.
  2. Use fldread to read open boundary data.
  • Property svn:keywords set to Id
File size: 15.5 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   !!----------------------------------------------------------------------
26
27   !!----------------------------------------------------------------------
28   !!   stp             : OPA system time-stepping
29   !!----------------------------------------------------------------------
30   USE step_oce         ! time stepping definition modules
31#if defined key_top
32   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
33#endif
34#if defined key_agrif
35   USE agrif_opa_sponge ! Momemtum and tracers sponges
36#endif
37   USE asminc           ! assimilation increments    (tra_asm_inc, dyn_asm_inc routines)
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      IF( lk_dtatem  )   CALL dta_tem( kstp )         ! update 3D temperature data
97      IF( lk_dtasal  )   CALL dta_sal( kstp )         ! update 3D salinity data
98                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
99      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
100
101      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
102      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
103      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
104                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
105
106      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
107      ! Ocean physics update                (ua, va, ta, sa used as workspace)
108      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
109                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
110                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
111      !
112      !  VERTICAL PHYSICS   
113                         CALL zdf_bfr( kstp )         ! bottom friction
114                         
115      !                                               ! Vertical eddy viscosity and diffusivity coefficients
116      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
117      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
118      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
119      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
120      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
121         avt (:,:,:) = rn_avt0 * tmask(:,:,:)
122         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
123         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
124      ENDIF
125      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
126         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
127      ENDIF
128      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
129
130      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
131
132      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
133         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
134         
135                         CALL zdf_mxl( kstp )         ! mixed layer depth
136
137                                                      ! write TKE or GLS information in the restart file
138      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
139      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
140      !
141      !  LATERAL  PHYSICS
142      !
143      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
144                         CALL eos( tsb, rhd )                ! before in situ density
145         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
146            &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
147         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
148                         CALL ldf_slp_grif( kstp )
149         ELSE
150                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
151         ENDIF
152      ENDIF
153#if defined key_traldf_c2d
154      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
155#endif
156
157      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
158      ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
159      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
160      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
161      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
162      IF( lk_diafwb  )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
163      IF( ln_diaptr  )   CALL dia_ptr( kstp )         ! Poleward TRansports diagnostics
164      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
165                         CALL dia_wri( kstp )         ! ocean model: outputs
166
167#if defined key_top
168      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
169      ! Passive Tracer Model
170      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
171                         CALL trc_stp( kstp )         ! time-stepping
172#endif
173
174      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
175      ! Active tracers                              (ua, va used as workspace)
176      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
177                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
178
179      IF(  ln_asmiau .AND. &
180         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
181                             CALL tra_sbc    ( kstp )       ! surface boundary condition
182      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
183      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
184      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
185      IF( lk_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
186                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection
187      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
188                             CALL tra_ldf    ( kstp )       ! lateral mixing
189#if defined key_agrif
190                             CALL tra_unswap
191      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
192                             CALL tra_swap
193#endif
194                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
195
196      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
197         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
198                             CALL tra_nxt( kstp )                ! tracer fields at next time step
199                             CALL eos    ( tsa, rhd, rhop )      ! Time-filtered in situ density for hpg computation
200         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative
201            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
202         
203      ELSE                                                  ! centered hpg  (eos then time stepping)
204                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation
205         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
206            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
207         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
208                             CALL tra_nxt( kstp )                ! tracer fields at next time step
209      ENDIF
210                             CALL tra_unswap                ! udate T & S 3D arrays  (to be suppressed)
211
212      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
213      ! Dynamics                                    (ta, sa 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                               CALL dyn_adv( kstp )         ! advection (vector or flux form)
221                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis
222                               CALL dyn_ldf( kstp )         ! lateral mixing
223#if defined key_agrif
224      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
225#endif
226                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
227                               CALL dyn_bfr( kstp )         ! bottom friction   
228                               CALL dyn_zdf( kstp )         ! vertical diffusion
229                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient
230                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step
231
232                               CALL ssh_nxt( kstp )         ! sea surface height at next time step
233
234      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
235      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
236
237      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
238      ! Control and restarts
239      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240                               CALL stp_ctl( kstp, indic )
241      IF( indic < 0        )   THEN
242                               CALL ctl_stop( 'step: indic < 0' )
243                               CALL dia_wri_state( 'output.abort', kstp )
244      ENDIF
245      IF( kstp == nit000   )   CALL iom_close( numror )     ! close input  ocean restart file
246      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
247
248      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
249      ! Trends                              (ua, va, ta, sa used as workspace)
250      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
251      IF( nstop == 0 ) THEN                         
252         IF( lk_trddyn     )   CALL trd_dwr( kstp )         ! trends: dynamics
253         IF( lk_trdtra     )   CALL trd_twr( kstp )         ! trends: active tracers
254         IF( lk_trdmld     )   CALL trd_mld( kstp )         ! trends: Mixed-layer
255         IF( lk_trdvor     )   CALL trd_vor( kstp )         ! trends: vorticity budget
256      ENDIF
257
258      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
259      ! Coupled mode
260      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
261      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
262      !
263      !
264   END SUBROUTINE stp
265
266   !!======================================================================
267END MODULE step
Note: See TracBrowser for help on using the repository browser.