source: branches/2013/dev_r3996_CMCC6_topbc/NEMOGCM/NEMO/OPA_SRC/step.F90 @ 4011

Last change on this file since 4011 was 4011, checked in by vichi, 8 years ago

Make a generic interface for trcdta when using other BGCM

This change introduces a more general trcdta structure that
is not strictly dependent on the number of tracers defined
in PISCES. The loop on the number of tracers is moved outside
trcdta and the tracer info and array is passed as an argument.
This allows to use trcdta as a library subroutine by the BFM and
other models.
NOTE: it must be tested throughly with all the PISCES configurations

This commit also updates the GYRE_BFM configuration and corrects
some minor missing cpp keys and real type definitions

  • Property svn:keywords set to Id
File size: 17.6 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
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   stp   ! called by opa.F90
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "zdfddm_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48#if defined key_agrif
49   SUBROUTINE stp( )
50      INTEGER             ::   kstp   ! ocean time-step index
51#else
52   SUBROUTINE stp( kstp )
53      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
54#endif
55      !!----------------------------------------------------------------------
56      !!                     ***  ROUTINE stp  ***
57      !!
58      !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
59      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
60      !!              - Tme stepping  of TRC (passive tracer eqs.)
61      !!
62      !! ** Method  : -1- Update forcings and data
63      !!              -2- Update ocean physics
64      !!              -3- Compute the t and s trends
65      !!              -4- Update t and s
66      !!              -5- Compute the momentum trends
67      !!              -6- Update the horizontal velocity
68      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
69      !!              -8- Outputs and diagnostics
70      !!----------------------------------------------------------------------
71      INTEGER ::   jk       ! dummy loop indice
72      INTEGER ::   indic    ! error indicator if < 0
73      !! ---------------------------------------------------------------------
74
75#if defined key_agrif
76      kstp = nit000 + Agrif_Nb_Step()
77!      IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
78!      IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
79# if defined key_iomput
80      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap
81# endif
82#endif
83                             indic = 0                ! reset to no error condition
84      IF( kstp == nit000 )   CALL iom_init            ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
85      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
86                             CALL iom_setkt( kstp - nit000 + 1 )   ! say to iom that we are at time step kstp
87
88      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
89      ! Update data, open boundaries, surface boundary condition (including sea-ice)
90      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
91                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
92      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp )
93      IF( lk_tide    )   CALL sbc_tide( kstp )
94      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
95      IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
96      IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries
97
98      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
99      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
100      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
101                         CALL ssh_wzv( kstp )         ! after ssh & vertical velocity
102
103      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
104      ! Ocean physics update                (ua, va, tsa used as workspace)
105      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
106                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
107                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
108      !
109      !  VERTICAL PHYSICS
110                         CALL zdf_bfr( kstp )         ! bottom friction
111
112      !                                               ! Vertical eddy viscosity and diffusivity coefficients
113      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
114      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
115      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
116      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
117      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
118         avt (:,:,:) = rn_avt0 * tmask(:,:,:)
119         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
120         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
121      ENDIF
122      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
123         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
124      ENDIF
125      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
126
127      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
128
129      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
130         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
131
132                         CALL zdf_mxl( kstp )         ! mixed layer depth
133
134                                                      ! write TKE or GLS information in the restart file
135      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
136      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
137      !
138      !  LATERAL  PHYSICS
139      !
140      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
141                         CALL eos( tsb, rhd )                ! before in situ density
142         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
143            &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
144         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
145                         CALL ldf_slp_grif( kstp )
146         ELSE
147                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
148         ENDIF
149      ENDIF
150#if defined key_traldf_c2d
151      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
152#endif
153#if defined key_traldf_c3d && key_traldf_smag
154                          CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
155#  endif
156#if defined key_dynldf_c3d && key_dynldf_smag
157                          CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
158#  endif
159
160
161      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
162      ! diagnostics and outputs             (ua, va, tsa used as workspace)
163      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
164      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
165      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
166      IF( lk_diafwb  )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
167      IF( ln_diaptr  )   CALL dia_ptr( kstp )         ! Poleward TRansports diagnostics
168      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
169      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
170      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
171                         CALL dia_wri( kstp )         ! ocean model: outputs
172
173#if defined key_top
174      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
175      ! Passive Tracer Model
176      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
177                         CALL trc_stp( kstp )         ! time-stepping
178#endif
179
180      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
181      ! Active tracers                              (ua, va used as workspace)
182      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
183                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
184
185!write(numout,*) "MAV kt",kstp
186!write(numout,'(a5,3(1x,f21.18))') "INIn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
187!write(numout,'(a5,3(1x,f21.18))') "INIa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
188      IF(  ln_asmiau .AND. &
189         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
190                             CALL tra_sbc    ( kstp )       ! surface boundary condition
191      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
192      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
193      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
194      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
195      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
196                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection
197!write(numout,'(a5,3(1x,f21.18))') "ADVn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
198!write(numout,'(a5,3(1x,f21.18))') "ADVa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
199      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
200                             CALL tra_ldf    ( kstp )       ! lateral mixing
201!write(numout,'(a5,3(1x,f21.18))') "LDFn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
202!write(numout,'(a5,3(1x,f21.18))') "LDFa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
203#if defined key_agrif
204      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
205#endif
206                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
207!do jk=1,jpk
208!write(numout,'(a5,3(1x,f21.18))') "ZDFn:",tsn(5,10,jk,jp_tem),tsn(5,10,jk,jp_sal),tmask(5,10,jk)
209!write(numout,'(a5,3(1x,f21.18))') "ZDFa:",tsa(5,10,jk,jp_tem),tsa(5,10,jk,jp_sal),ssha(5,10)
210!end do
211
212      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
213         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
214                             CALL tra_nxt( kstp )                ! tracer fields at next time step
215                             CALL eos    ( tsa, rhd, rhop )      ! Time-filtered in situ density for hpg computation
216         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative
217            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
218
219      ELSE                                                  ! centered hpg  (eos then time stepping)
220                             CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation
221         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative
222            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
223!write(numout,'(a5,3(1x,f21.18))') "ZPSn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)
224!write(numout,'(a5,3(1x,f21.18))') "ZPSa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
225         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
226                             CALL tra_nxt( kstp )                ! tracer fields at next time step
227!write(numout,'(a5,3(1x,f21.18))') "NXTn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(25,11)
228!write(numout,'(a5,3(1x,f21.18))') "NXTa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)
229      ENDIF
230
231      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
232      ! Dynamics                                    (tsa used as workspace)
233      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
234                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero
235                               va(:,:,:) = 0.e0
236
237      IF(  ln_asmiau .AND. &
238         & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
239      IF( ln_bkgwri )          CALL asm_bkg_wri( kstp )     ! output background fields
240      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
241      IF( lk_bdy           )   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
242                               CALL dyn_adv( kstp )         ! advection (vector or flux form)
243                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis
244                               CALL dyn_ldf( kstp )         ! lateral mixing
245      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
246#if defined key_agrif
247      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
248#endif
249                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
250                               CALL dyn_bfr( kstp )         ! bottom friction
251                               CALL dyn_zdf( kstp )         ! vertical diffusion
252                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient
253                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step
254
255                               CALL ssh_nxt( kstp )         ! sea surface height at next time step
256
257      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
258      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
259
260      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
261      ! Control and restarts
262      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
263                               CALL stp_ctl( kstp, indic )
264      IF( indic < 0        )   THEN
265                               CALL ctl_stop( 'step: indic < 0' )
266                               CALL dia_wri_state( 'output.abort', kstp )
267      ENDIF
268      IF( kstp == nit000   )   CALL iom_close( numror )     ! close input  ocean restart file
269      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
270      IF( lk_obc           )   CALL obc_rst_write( kstp )   ! write open boundary restart file
271
272      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
273      ! Trends                              (ua, va, tsa used as workspace)
274      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
275      IF( nstop == 0 ) THEN
276         IF( lk_trddyn     )   CALL trd_dwr( kstp )         ! trends: dynamics
277         IF( lk_trdtra     )   CALL trd_twr( kstp )         ! trends: active tracers
278         IF( lk_trdmld     )   CALL trd_mld( kstp )         ! trends: Mixed-layer
279         IF( lk_trdvor     )   CALL trd_vor( kstp )         ! trends: vorticity budget
280      ENDIF
281
282      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
283      ! Coupled mode
284      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
285      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
286      !
287#if defined key_iomput
288      IF( kstp == nitend .OR. indic < 0 )   CALL xios_context_finalize() ! needed for XIOS+AGRIF
289#endif
290      !
291      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
292      !
293   END SUBROUTINE stp
294
295   !!======================================================================
296END MODULE step
Note: See TracBrowser for help on using the repository browser.