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.
trcsub.F90 in NEMO/trunk/src/TOP – NEMO

source: NEMO/trunk/src/TOP/trcsub.F90 @ 11612

Last change on this file since 11612 was 10425, checked in by smasson, 5 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 29.0 KB
Line 
1MODULE trcsub
2   !!======================================================================
3   !!                       ***  MODULE trcsubstp  ***
4   !!   TOP :   Averages physics variables for TOP substepping.
5   !!======================================================================
6   !! History :  1.0  !  2011-10  (K. Edwards)  Original
7   !!----------------------------------------------------------------------
8#if defined key_top
9   !!----------------------------------------------------------------------
10   !!   trc_sub       : passive tracer system sub-stepping
11   !!----------------------------------------------------------------------
12   USE oce_trc        ! ocean dynamics and active tracers variables
13   USE trc
14   USE trabbl         ! bottom boundary layer
15   USE zdf_oce
16   USE domvvl
17   USE divhor         ! horizontal divergence
18   USE sbcrnf   , ONLY: h_rnf, nk_rnf    ! River runoff
19   USE bdy_oce  , ONLY: ln_bdy, bdytmask ! BDY
20   !
21   USE prtctl_trc     ! Print control for debbuging
22   USE in_out_manager !
23   USE iom
24   USE lbclnk
25#if defined key_agrif
26   USE agrif_oce_update
27   USE agrif_oce_interp
28#endif
29
30   IMPLICIT NONE
31
32   PUBLIC   trc_sub_stp     ! called by trc_stp
33   PUBLIC   trc_sub_ini     ! called by trc_ini to initialize substepping arrays.
34   PUBLIC   trc_sub_reset   ! called by trc_stp to reset physics variables
35   PUBLIC   trc_sub_ssh     ! called by trc_stp to reset physics variables
36
37   REAL(wp) :: r1_ndttrc     ! = 1 /  nn_dttrc
38   REAL(wp) :: r1_ndttrcp1   ! = 1 / (nn_dttrc+1)
39
40
41   !! averaged and temporary saved variables  (needed when a larger passive tracer time-step is used)
42   !! ----------------------------------------------------------------
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::    un_tm ,   un_temp   !: i-horizontal velocity average     [m/s]
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::    vn_tm ,   vn_temp   !: j-horizontal velocity average     [m/s]
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::              wn_temp   !: hold current values of avt, un, vn, wn
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsn_tm ,  tsn_temp   !: t/s average     [m/s]
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   avs_tm ,  avs_temp   !: vertical diffusivity coeff. at  w-point   [m2/s]
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  rhop_tm , rhop_temp   !:
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  sshn_tm , sshn_temp   !: average ssh for the now step [m]
50
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::    rnf_tm ,    rnf_temp   !: river runoff
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  h_rnf_tm ,  h_rnf_temp   !: depth in metres to the bottom of the relevant grid box
53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hmld_tm ,   hmld_temp   !: mixed layer depth average [m]
54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   fr_i_tm ,   fr_i_temp   !: average ice fraction     [m/s]
55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::    emp_tm ,    emp_temp   !: freshwater budget: volume flux [Kg/m2/s]
56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     :: fmmflx_tm , fmmflx_temp   !: freshwater budget: freezing/melting [Kg/m2/s]
57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     :: emp_b_hold,  emp_b_temp   !: hold emp from the beginning of each sub-stepping[m] 
58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::    qsr_tm ,    qsr_temp   !: solar radiation average [m]
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   wndm_tm ,   wndm_temp   !: 10m wind average [m]
60   !
61   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  sshb_hold   !:hold sshb from the beginning of each sub-stepping[m] 
62   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sshb_temp, ssha_temp
63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  hdivn_temp, rotn_temp
64   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  hdivb_temp, rotb_temp
65   !
66   !                                                    !!- bottom boundary layer param (ln_trabbl=T)
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  ahu_bbl_tm, ahu_bbl_temp  ! BBL diffusive i-coef.
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  ahv_bbl_tm, ahv_bbl_temp  ! BBL diffusive j-coef.
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  utr_bbl_tm, utr_bbl_temp  ! BBL u-advection
70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  vtr_bbl_tm, vtr_bbl_temp  ! BBL v-advection
71
72   !                                                      !!- iso-neutral slopes (if l_ldfslp=T)
73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   uslp_temp, vslp_temp, wslpi_temp, wslpj_temp   !: hold current values
74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   uslp_tm  , vslp_tm  , wslpi_tm  , wslpj_tm     !: time mean
75
76
77   !!----------------------------------------------------------------------
78   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
79   !! $Id$
80   !! Software governed by the CeCILL license (see ./LICENSE)
81   !!----------------------------------------------------------------------
82CONTAINS
83
84   SUBROUTINE trc_sub_stp( kt )
85      !!-------------------------------------------------------------------
86      !!                     ***  ROUTINE trc_stp  ***
87      !!                     
88      !! ** Purpose : Average variables needed for sub-stepping passive tracers
89      !!
90      !! ** Method  : Called every timestep to increment _tm (time mean) variables
91      !!              on TOP steps, calculate averages.
92      !!-------------------------------------------------------------------
93      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
94      !
95      INTEGER ::   ji, jj, jk   ! dummy loop indices
96      REAL(wp)::   z1_ne3t, z1_ne3u, z1_ne3v, z1_ne3w   ! local scalars
97      !!-------------------------------------------------------------------
98      !
99      IF( ln_timing )   CALL timing_start('trc_sub_stp')
100      !
101      IF( kt == nit000 ) THEN
102           IF(lwp) WRITE(numout,*)
103           IF(lwp) WRITE(numout,*) 'trc_sub_stp : substepping of the passive tracers'
104           IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
105           !
106           sshb_hold  (:,:) = sshn  (:,:)
107           emp_b_hold (:,:) = emp_b (:,:)
108           !
109           r1_ndttrc        = 1._wp / REAL( nn_dttrc    , wp ) 
110           r1_ndttrcp1      = 1._wp / REAL( nn_dttrc + 1, wp )
111      ENDIF 
112
113      IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
114         !
115         un_tm   (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * e3u_n(:,:,:) 
116         vn_tm   (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_n(:,:,:) 
117         tsn_tm  (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
118         tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
119         rhop_tm (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:) 
120         avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:) 
121         IF( l_ldfslp ) THEN
122            uslp_tm (:,:,:)      = uslp_tm (:,:,:)        + uslp (:,:,:)
123            vslp_tm (:,:,:)      = vslp_tm (:,:,:)        + vslp (:,:,:)
124            wslpi_tm(:,:,:)      = wslpi_tm(:,:,:)        + wslpi(:,:,:)
125            wslpj_tm(:,:,:)      = wslpj_tm(:,:,:)        + wslpj(:,:,:)
126         ENDIF
127         IF( ln_trabbl ) THEN
128            IF( nn_bbl_ldf == 1 ) THEN
129               ahu_bbl_tm(:,:)     = ahu_bbl_tm(:,:)        + ahu_bbl(:,:) 
130               ahv_bbl_tm(:,:)     = ahv_bbl_tm(:,:)        + ahv_bbl(:,:) 
131            ENDIF
132            IF( nn_bbl_adv == 1 ) THEN
133               utr_bbl_tm(:,:)     = utr_bbl_tm(:,:)        + utr_bbl(:,:) 
134               vtr_bbl_tm(:,:)     = vtr_bbl_tm(:,:)        + vtr_bbl(:,:) 
135            ENDIF
136         ENDIF 
137         !
138         sshn_tm  (:,:)         = sshn_tm  (:,:)         + sshn  (:,:) 
139         rnf_tm   (:,:)         = rnf_tm   (:,:)         + rnf   (:,:) 
140         h_rnf_tm (:,:)         = h_rnf_tm (:,:)         + h_rnf (:,:) 
141         hmld_tm  (:,:)         = hmld_tm  (:,:)         + hmld  (:,:)
142         fr_i_tm  (:,:)         = fr_i_tm  (:,:)         + fr_i  (:,:)
143         emp_tm   (:,:)         = emp_tm   (:,:)         + emp   (:,:) 
144         fmmflx_tm(:,:)         = fmmflx_tm(:,:)         + fmmflx(:,:)
145         qsr_tm   (:,:)         = qsr_tm   (:,:)         + qsr   (:,:)
146         wndm_tm  (:,:)         = wndm_tm  (:,:)         + wndm  (:,:)
147         !
148      ELSE                           !  It is time to substep
149         !   1. set temporary arrays to hold physics/dynamical variables
150         un_temp    (:,:,:)      = un    (:,:,:)
151         vn_temp    (:,:,:)      = vn    (:,:,:)
152         wn_temp    (:,:,:)      = wn    (:,:,:)
153         tsn_temp   (:,:,:,:)    = tsn   (:,:,:,:)
154         rhop_temp  (:,:,:)      = rhop  (:,:,:)   
155         avs_temp   (:,:,:)      = avs   (:,:,:)
156         IF( l_ldfslp ) THEN
157            uslp_temp  (:,:,:)   = uslp  (:,:,:)   ;   wslpi_temp (:,:,:)   = wslpi (:,:,:)
158            vslp_temp  (:,:,:)   = vslp  (:,:,:)   ;   wslpj_temp (:,:,:)   = wslpj (:,:,:)
159         ENDIF
160         IF( ln_trabbl ) THEN
161            IF( nn_bbl_ldf == 1 ) THEN
162               ahu_bbl_temp(:,:)   = ahu_bbl(:,:) 
163               ahv_bbl_temp(:,:)   = ahv_bbl(:,:) 
164            ENDIF
165            IF( nn_bbl_adv == 1 ) THEN
166               utr_bbl_temp(:,:)   = utr_bbl(:,:) 
167               vtr_bbl_temp(:,:)   = vtr_bbl(:,:) 
168            ENDIF
169         ENDIF
170         sshn_temp  (:,:)        = sshn  (:,:)
171         sshb_temp  (:,:)        = sshb  (:,:)
172         ssha_temp  (:,:)        = ssha  (:,:)
173         rnf_temp   (:,:)        = rnf   (:,:)
174         h_rnf_temp (:,:)        = h_rnf (:,:)
175         hmld_temp  (:,:)        = hmld  (:,:)
176         fr_i_temp  (:,:)        = fr_i  (:,:)
177         emp_temp   (:,:)        = emp   (:,:)
178         emp_b_temp (:,:)        = emp_b (:,:)
179         fmmflx_temp(:,:)        = fmmflx(:,:)
180         qsr_temp   (:,:)        = qsr   (:,:)
181         wndm_temp  (:,:)        = wndm  (:,:)
182         !                                    !  Variables reset in trc_sub_ssh
183         hdivn_temp (:,:,:)      = hdivn (:,:,:)
184         !
185         ! 2. Create averages and reassign variables
186         un_tm    (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * e3u_n(:,:,:) 
187         vn_tm    (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_n(:,:,:) 
188         tsn_tm   (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
189         tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
190         rhop_tm (:,:,:)         = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:) 
191         avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:) 
192         IF( l_ldfslp ) THEN
193            uslp_tm  (:,:,:)     = uslp_tm (:,:,:)        + uslp (:,:,:) 
194            vslp_tm  (:,:,:)     = vslp_tm (:,:,:)        + vslp (:,:,:)
195            wslpi_tm (:,:,:)     = wslpi_tm(:,:,:)        + wslpi(:,:,:) 
196            wslpj_tm (:,:,:)     = wslpj_tm(:,:,:)        + wslpj(:,:,:) 
197         ENDIF
198         IF( ln_trabbl ) THEN
199            IF( nn_bbl_ldf == 1 ) THEN
200               ahu_bbl_tm(:,:)     = ahu_bbl_tm(:,:)        + ahu_bbl(:,:) 
201               ahv_bbl_tm(:,:)     = ahv_bbl_tm(:,:)        + ahv_bbl(:,:) 
202            ENDIF
203            IF( nn_bbl_adv == 1 ) THEN
204               utr_bbl_tm(:,:)     = utr_bbl_tm(:,:)        + utr_bbl(:,:) 
205               vtr_bbl_tm(:,:)     = vtr_bbl_tm(:,:)        + vtr_bbl(:,:) 
206            ENDIF
207         ENDIF
208         sshn_tm  (:,:)          = sshn_tm    (:,:)       + sshn  (:,:) 
209         rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:) 
210         h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:) 
211         hmld_tm  (:,:)          = hmld_tm    (:,:)       + hmld  (:,:)
212         fr_i_tm  (:,:)          = fr_i_tm    (:,:)       + fr_i  (:,:)
213         emp_tm   (:,:)          = emp_tm     (:,:)       + emp   (:,:) 
214         fmmflx_tm(:,:)          = fmmflx_tm  (:,:)       + fmmflx(:,:)
215         qsr_tm   (:,:)          = qsr_tm     (:,:)       + qsr   (:,:)
216         wndm_tm  (:,:)          = wndm_tm    (:,:)       + wndm  (:,:)
217         !
218         sshn     (:,:)          = sshn_tm    (:,:) * r1_ndttrcp1 
219         sshb     (:,:)          = sshb_hold  (:,:)
220         rnf      (:,:)          = rnf_tm     (:,:) * r1_ndttrcp1 
221         h_rnf    (:,:)          = h_rnf_tm   (:,:) * r1_ndttrcp1 
222         hmld     (:,:)          = hmld_tm    (:,:) * r1_ndttrcp1 
223         !  variables that are initialized after averages
224         emp_b    (:,:) = emp_b_hold (:,:)
225         IF( kt == nittrc000 ) THEN
226            wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrc 
227            qsr   (:,:)          = qsr_tm     (:,:) * r1_ndttrc 
228            emp   (:,:)          = emp_tm     (:,:) * r1_ndttrc 
229            fmmflx(:,:)          = fmmflx_tm  (:,:) * r1_ndttrc 
230            fr_i  (:,:)          = fr_i_tm    (:,:) * r1_ndttrc
231            IF( ln_trabbl ) THEN
232               IF( nn_bbl_ldf == 1 ) THEN
233                  ahu_bbl(:,:)      = ahu_bbl_tm (:,:) * r1_ndttrc 
234                  ahv_bbl(:,:)      = ahv_bbl_tm (:,:) * r1_ndttrc 
235               ENDIF
236               IF( nn_bbl_adv == 1 ) THEN
237                  utr_bbl(:,:)      = utr_bbl_tm (:,:) * r1_ndttrc 
238                  vtr_bbl(:,:)      = vtr_bbl_tm (:,:) * r1_ndttrc 
239               ENDIF
240            ENDIF
241         ELSE
242            wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrcp1 
243            qsr   (:,:)          = qsr_tm     (:,:) * r1_ndttrcp1 
244            emp   (:,:)          = emp_tm     (:,:) * r1_ndttrcp1 
245            fmmflx(:,:)          = fmmflx_tm  (:,:) * r1_ndttrcp1 
246            fr_i  (:,:)          = fr_i_tm    (:,:) * r1_ndttrcp1 
247            IF( ln_trabbl ) THEN
248               IF( nn_bbl_ldf == 1 ) THEN
249                  ahu_bbl(:,:)      = ahu_bbl_tm (:,:) * r1_ndttrcp1 
250                  ahv_bbl(:,:)      = ahv_bbl_tm (:,:) * r1_ndttrcp1 
251               ENDIF
252               IF( nn_bbl_adv == 1 ) THEN
253                  utr_bbl(:,:)      = utr_bbl_tm (:,:) * r1_ndttrcp1 
254                  vtr_bbl(:,:)      = vtr_bbl_tm (:,:) * r1_ndttrcp1 
255               ENDIF
256            ENDIF
257         ENDIF
258         !
259         DO jk = 1, jpk
260            DO jj = 1, jpj
261               DO ji = 1, jpi
262                  z1_ne3t = r1_ndttrcp1  / e3t_n(ji,jj,jk)
263                  z1_ne3u = r1_ndttrcp1  / e3u_n(ji,jj,jk)
264                  z1_ne3v = r1_ndttrcp1  / e3v_n(ji,jj,jk)
265                  z1_ne3w = r1_ndttrcp1  / e3w_n(ji,jj,jk)
266                  !
267                  un   (ji,jj,jk)        = un_tm   (ji,jj,jk)        * z1_ne3u
268                  vn   (ji,jj,jk)        = vn_tm   (ji,jj,jk)        * z1_ne3v
269                  tsn  (ji,jj,jk,jp_tem) = tsn_tm  (ji,jj,jk,jp_tem) * z1_ne3t
270                  tsn  (ji,jj,jk,jp_sal) = tsn_tm  (ji,jj,jk,jp_sal) * z1_ne3t
271                  rhop (ji,jj,jk)        = rhop_tm (ji,jj,jk)        * z1_ne3t
272!!gm : BUG ==>> for avs I don't understand the division by e3w
273                  avs  (ji,jj,jk)        = avs_tm  (ji,jj,jk)        * z1_ne3w
274               END DO
275            END DO
276         END DO
277         IF( l_ldfslp ) THEN
278            wslpi(:,:,:)        = wslpi_tm(:,:,:) 
279            wslpj(:,:,:)        = wslpj_tm(:,:,:)
280            uslp (:,:,:)        = uslp_tm (:,:,:)
281            vslp (:,:,:)        = vslp_tm (:,:,:)
282         ENDIF
283         !
284         CALL trc_sub_ssh( kt )         ! after ssh & vertical velocity
285         !
286      ENDIF
287      !
288      IF( ln_timing )   CALL timing_stop('trc_sub_stp')
289      !
290   END SUBROUTINE trc_sub_stp
291
292
293   SUBROUTINE trc_sub_ini
294      !!-------------------------------------------------------------------
295      !!                     ***  ROUTINE trc_sub_ini  ***
296      !!                     
297      !! ** Purpose : Initialize variables needed for sub-stepping passive tracers
298      !!
299      !! ** Method  :
300      !!              Compute the averages for sub-stepping
301      !!-------------------------------------------------------------------
302      INTEGER ::   ierr
303      !!-------------------------------------------------------------------
304      !
305      IF(lwp) WRITE(numout,*)
306      IF(lwp) WRITE(numout,*) 'trc_sub_ini : initial set up of the passive tracers substepping'
307      IF(lwp) WRITE(numout,*) '~~~~~~~'
308
309      ierr =  trc_sub_alloc    ()
310      CALL mpp_sum( 'trcsub', ierr )
311      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
312
313      un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
314      vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
315      tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
316      tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
317      rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:) 
318!!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
319      avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:) 
320      IF( l_ldfslp ) THEN
321         wslpi_tm(:,:,:)     = wslpi(:,:,:)
322         wslpj_tm(:,:,:)     = wslpj(:,:,:)
323         uslp_tm (:,:,:)     = uslp (:,:,:)
324         vslp_tm (:,:,:)     = vslp (:,:,:)
325      ENDIF
326      sshn_tm  (:,:) = sshn  (:,:) 
327      rnf_tm   (:,:) = rnf   (:,:) 
328      h_rnf_tm (:,:) = h_rnf (:,:) 
329      hmld_tm  (:,:) = hmld  (:,:)
330
331      ! Physics variables that are set after initialization:
332      fr_i_tm  (:,:) = 0._wp
333      emp_tm   (:,:) = 0._wp
334      fmmflx_tm(:,:)  = 0._wp
335      qsr_tm   (:,:) = 0._wp
336      wndm_tm  (:,:) = 0._wp
337      IF( ln_trabbl ) THEN
338         IF( nn_bbl_ldf == 1 ) THEN
339            ahu_bbl_tm(:,:) = 0._wp
340            ahv_bbl_tm(:,:) = 0._wp
341         ENDIF
342         IF( nn_bbl_adv == 1 ) THEN
343            utr_bbl_tm(:,:) = 0._wp
344            vtr_bbl_tm(:,:) = 0._wp
345         ENDIF
346      ENDIF
347      !
348   END SUBROUTINE trc_sub_ini
349
350
351   SUBROUTINE trc_sub_reset( kt )
352      !!-------------------------------------------------------------------
353      !!                     ***  ROUTINE trc_sub_reset  ***
354      !!                     
355      !! ** Purpose : Reset physics variables averaged for substepping
356      !!
357      !! ** Method  :
358      !!              Compute the averages for sub-stepping
359      !!-------------------------------------------------------------------
360      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
361      INTEGER :: jk                 ! dummy loop indices
362      !!-------------------------------------------------------------------
363      !
364      IF( ln_timing )   CALL timing_start('trc_sub_reset')
365      !
366      !   restore physics variables
367      un    (:,:,:)   =  un_temp    (:,:,:)
368      vn    (:,:,:)   =  vn_temp    (:,:,:)
369      wn    (:,:,:)   =  wn_temp    (:,:,:)
370      tsn   (:,:,:,:) =  tsn_temp   (:,:,:,:)
371      rhop  (:,:,:)   =  rhop_temp  (:,:,:)
372      avs   (:,:,:)   =  avs_temp   (:,:,:)
373      IF( l_ldfslp ) THEN
374         wslpi (:,:,:)=  wslpi_temp (:,:,:)
375         wslpj (:,:,:)=  wslpj_temp (:,:,:)
376         uslp  (:,:,:)=  uslp_temp  (:,:,:)
377         vslp  (:,:,:)=  vslp_temp  (:,:,:)
378      ENDIF
379      sshn  (:,:)     =  sshn_temp  (:,:)
380      sshb  (:,:)     =  sshb_temp  (:,:)
381      ssha  (:,:)     =  ssha_temp  (:,:)
382      rnf   (:,:)     =  rnf_temp   (:,:)
383      h_rnf (:,:)     =  h_rnf_temp (:,:)
384      !
385      hmld  (:,:)     =  hmld_temp  (:,:)
386      fr_i  (:,:)     =  fr_i_temp  (:,:)
387      emp   (:,:)     =  emp_temp   (:,:)
388      fmmflx(:,:)     =  fmmflx_temp(:,:)
389      emp_b (:,:)     =  emp_b_temp (:,:)
390      qsr   (:,:)     =  qsr_temp   (:,:)
391      wndm  (:,:)     =  wndm_temp  (:,:)
392      IF( ln_trabbl ) THEN
393         IF( nn_bbl_ldf == 1 ) THEN
394            ahu_bbl(:,:) = ahu_bbl_temp(:,:) 
395            ahv_bbl(:,:) = ahv_bbl_temp(:,:) 
396         ENDIF
397         IF( nn_bbl_adv == 1 ) THEN
398            utr_bbl(:,:) = utr_bbl_temp(:,:) 
399            vtr_bbl(:,:) = vtr_bbl_temp(:,:) 
400         ENDIF
401      ENDIF
402      !
403      hdivn (:,:,:)   =  hdivn_temp (:,:,:)
404      !                                     
405      ! Start new averages
406         un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
407         vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
408         tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
409         tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
410         rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:) 
411         avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:) 
412      IF( l_ldfslp ) THEN
413         uslp_tm (:,:,:)        = uslp (:,:,:)
414         vslp_tm (:,:,:)        = vslp (:,:,:)
415         wslpi_tm(:,:,:)        = wslpi(:,:,:) 
416         wslpj_tm(:,:,:)        = wslpj(:,:,:)
417      ENDIF
418      !
419      sshb_hold  (:,:) = sshn  (:,:)
420      emp_b_hold (:,:) = emp   (:,:)
421      sshn_tm    (:,:) = sshn  (:,:) 
422      rnf_tm     (:,:) = rnf   (:,:) 
423      h_rnf_tm   (:,:) = h_rnf (:,:) 
424      hmld_tm    (:,:) = hmld  (:,:)
425      fr_i_tm    (:,:) = fr_i  (:,:)
426      emp_tm     (:,:) = emp   (:,:)
427      fmmflx_tm  (:,:) = fmmflx(:,:)
428      qsr_tm     (:,:) = qsr   (:,:)
429      wndm_tm    (:,:) = wndm  (:,:)
430      IF( ln_trabbl ) THEN
431         IF( nn_bbl_ldf == 1 ) THEN
432            ahu_bbl_tm(:,:) = ahu_bbl(:,:) 
433            ahv_bbl_tm(:,:) = ahv_bbl(:,:) 
434         ENDIF
435         IF( nn_bbl_adv == 1 ) THEN
436            utr_bbl_tm(:,:) = utr_bbl(:,:) 
437            vtr_bbl_tm(:,:) = vtr_bbl(:,:) 
438         ENDIF
439      ENDIF
440      !
441      !
442      IF( ln_timing )   CALL timing_stop('trc_sub_reset')
443      !
444   END SUBROUTINE trc_sub_reset
445
446
447   SUBROUTINE trc_sub_ssh( kt ) 
448      !!----------------------------------------------------------------------
449      !!                ***  ROUTINE trc_sub_ssh  ***
450      !!                   
451      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
452      !!              and update the now vertical coordinate (ln_linssh=F).
453      !!
454      !! ** Method  : - Using the incompressibility hypothesis, the vertical
455      !!      velocity is computed by integrating the horizontal divergence 
456      !!      from the bottom to the surface minus the scale factor evolution.
457      !!        The boundary conditions are w=0 at the bottom (no flux) and.
458      !!
459      !! ** action  :   ssha    : after sea surface height
460      !!                wn      : now vertical velocity
461      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (ln_linssh=F)
462      !!
463      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
464      !!----------------------------------------------------------------------
465      INTEGER, INTENT(in) ::   kt   ! time step
466      !
467      INTEGER  ::   ji, jj, jk   ! dummy loop indices
468      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
469      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
470      !!---------------------------------------------------------------------
471      !
472      IF( ln_timing )   CALL timing_start('trc_sub_ssh')
473      !
474
475      IF( kt == nittrc000 ) THEN
476         !
477         IF(lwp) WRITE(numout,*)
478         IF(lwp) WRITE(numout,*) 'trc_sub_ssh : after sea surface height and now vertical velocity '
479         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
480         !
481         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
482         !
483      ENDIF
484      !
485!!gm BUG here !   hdivn will include the runoff divergence at the wrong timestep !!!!
486      CALL div_hor( kt )                              ! Horizontal divergence & Relative vorticity
487      !
488      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
489      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
490
491      !                                           !------------------------------!
492      !                                           !   After Sea Surface Height   !
493      !                                           !------------------------------!
494      zhdiv(:,:) = 0._wp
495      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
496        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
497      END DO
498      !                                                ! Sea surface elevation time stepping
499      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
500      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
501      z1_rau0 = 0.5 / rau0
502      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
503
504      IF( .NOT.ln_dynspg_ts ) THEN
505      ! These lines are not necessary with time splitting since
506      ! boundary condition on sea level is set during ts loop
507#if defined key_agrif
508      CALL agrif_ssh( kt )
509#endif
510         IF( ln_bdy ) THEN
511            ssha(:,:) = ssha(:,:) * bdytmask(:,:)
512            CALL lbc_lnk( 'trcsub', ssha, 'T', 1. ) 
513         ENDIF
514      ENDIF
515      !
516      !                                           !------------------------------!
517      !                                           !     Now Vertical Velocity    !
518      !                                           !------------------------------!
519      z1_2dt = 1.e0 / z2dt
520      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
521         ! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise
522         wn(:,:,jk) = wn(:,:,jk+1) -   e3t_n(:,:,jk) * hdivn(:,:,jk)        &
523            &                      - ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )    &
524            &                         * tmask(:,:,jk) * z1_2dt
525         IF( ln_bdy ) wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
526      END DO
527      !
528      IF( ln_timing )   CALL timing_stop('trc_sub_ssh')
529      !
530   END SUBROUTINE trc_sub_ssh
531
532
533   INTEGER FUNCTION trc_sub_alloc()
534      !!-------------------------------------------------------------------
535      !!                    *** ROUTINE trc_sub_alloc ***
536      !!-------------------------------------------------------------------
537      USE lib_mpp, ONLY: ctl_stop
538      INTEGER ::  ierr(3)
539      !!-------------------------------------------------------------------
540      !
541      ierr(:) = 0
542      !
543      ALLOCATE( un_temp(jpi,jpj,jpk)      ,  vn_temp(jpi,jpj,jpk)   ,     &
544         &      wn_temp(jpi,jpj,jpk)      ,                               &
545         &      rhop_temp(jpi,jpj,jpk)    ,  rhop_tm(jpi,jpj,jpk)   ,     &
546         &      sshn_temp(jpi,jpj)        ,  sshb_temp(jpi,jpj)     ,     &
547         &      ssha_temp(jpi,jpj)        ,                               &
548         &      rnf_temp(jpi,jpj)         ,  h_rnf_temp(jpi,jpj)    ,     &
549         &      tsn_temp(jpi,jpj,jpk,2)   ,  emp_b_temp(jpi,jpj)    ,     &
550         &      emp_temp(jpi,jpj)         ,  fmmflx_temp(jpi,jpj)   ,     &
551         &      hmld_temp(jpi,jpj)        ,  qsr_temp(jpi,jpj)      ,     &
552         &      fr_i_temp(jpi,jpj)        ,  fr_i_tm(jpi,jpj)       ,     &
553         &      wndm_temp(jpi,jpj)        ,  wndm_tm(jpi,jpj)       ,     &
554         &      avs_tm(jpi,jpj,jpk)       ,  avs_temp(jpi,jpj,jpk)  ,     &
555         &      hdivn_temp(jpi,jpj,jpk)   ,  hdivb_temp(jpi,jpj,jpk),     &
556         &      un_tm(jpi,jpj,jpk)        ,  vn_tm(jpi,jpj,jpk)     ,     &
557         &      sshn_tm(jpi,jpj)          ,  sshb_hold(jpi,jpj)     ,     &
558         &      tsn_tm(jpi,jpj,jpk,2)     ,                               &
559         &      emp_tm(jpi,jpj)           ,  fmmflx_tm(jpi,jpj)     ,     &
560         &      emp_b_hold(jpi,jpj)       ,                               &
561         &      hmld_tm(jpi,jpj)          ,  qsr_tm(jpi,jpj)        ,     &
562         &      rnf_tm(jpi,jpj)           ,  h_rnf_tm(jpi,jpj)      , STAT=ierr(1) ) 
563      !
564      IF( l_ldfslp ) THEN
565         ALLOCATE( uslp_temp(jpi,jpj,jpk) ,  wslpi_temp(jpi,jpj,jpk),     &
566            &      vslp_temp(jpi,jpj,jpk) ,  wslpj_temp(jpi,jpj,jpk),     &
567            &      uslp_tm  (jpi,jpj,jpk) ,  wslpi_tm  (jpi,jpj,jpk),     &
568            &      vslp_tm  (jpi,jpj,jpk) ,  wslpj_tm  (jpi,jpj,jpk), STAT=ierr(2) )
569      ENDIF
570      IF( ln_trabbl ) THEN
571         ALLOCATE( ahu_bbl_temp(jpi,jpj)  , utr_bbl_temp(jpi,jpj)   ,     &
572            &      ahv_bbl_temp(jpi,jpj)  , vtr_bbl_temp(jpi,jpj)   ,     &
573            &      ahu_bbl_tm  (jpi,jpj)  , utr_bbl_tm  (jpi,jpj)   ,     &
574            &      ahv_bbl_tm  (jpi,jpj)  , vtr_bbl_tm  (jpi,jpj)   , STAT=ierr(3) ) 
575      ENDIF
576      !
577      trc_sub_alloc = MAXVAL( ierr )
578      !
579      IF( trc_sub_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trc_sub_alloc: failed to allocate arrays' )
580      !
581   END FUNCTION trc_sub_alloc
582
583#else
584   !!----------------------------------------------------------------------
585   !!   Default key                                     NO passive tracers
586   !!----------------------------------------------------------------------
587CONTAINS
588   SUBROUTINE trc_sub_stp( kt )        ! Empty routine
589      WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
590   END SUBROUTINE trc_sub_stp
591   SUBROUTINE trc_sub_ini        ! Empty routine
592      WRITE(*,*) 'trc_sub_ini: You should not have seen this print! error?', kt
593   END SUBROUTINE trc_sub_ini
594#endif
595
596   !!======================================================================
597END MODULE trcsub
Note: See TracBrowser for help on using the repository browser.