source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90 @ 10963

Last change on this file since 10963 was 10963, checked in by acc, 19 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TOP routines in top-level TOP directory and all knock on effects of these conversions. SETTE tested.

  • Property svn:keywords set to Id
File size: 29.9 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, uu(Kmm), vv(Kmm), ww
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 ssh(Kbb) 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, Kbb, Kmm, Krhs )
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      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! ocean time-level index
95      !
96      INTEGER ::   ji, jj, jk   ! dummy loop indices
97      REAL(wp)::   z1_ne3t, z1_ne3u, z1_ne3v, z1_ne3w   ! local scalars
98      !!-------------------------------------------------------------------
99      !
100      IF( ln_timing )   CALL timing_start('trc_sub_stp')
101      !
102      IF( kt == nit000 ) THEN
103           IF(lwp) WRITE(numout,*)
104           IF(lwp) WRITE(numout,*) 'trc_sub_stp : substepping of the passive tracers'
105           IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
106           !
107           sshb_hold  (:,:) = ssh   (:,:,Kmm)
108           emp_b_hold (:,:) = emp_b (:,:)
109           !
110           r1_ndttrc        = 1._wp / REAL( nn_dttrc    , wp ) 
111           r1_ndttrcp1      = 1._wp / REAL( nn_dttrc + 1, wp )
112      ENDIF 
113
114      IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
115         !
116         un_tm   (:,:,:)        = un_tm   (:,:,:)        + uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
117         vn_tm   (:,:,:)        = vn_tm   (:,:,:)        + vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
118         tsn_tm  (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) 
119         tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) 
120         rhop_tm (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)            * e3t(:,:,:,Kmm) 
121         avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)            * e3w(:,:,:,Kmm) 
122         IF( l_ldfslp ) THEN
123            uslp_tm (:,:,:)      = uslp_tm (:,:,:)        + uslp (:,:,:)
124            vslp_tm (:,:,:)      = vslp_tm (:,:,:)        + vslp (:,:,:)
125            wslpi_tm(:,:,:)      = wslpi_tm(:,:,:)        + wslpi(:,:,:)
126            wslpj_tm(:,:,:)      = wslpj_tm(:,:,:)        + wslpj(:,:,:)
127         ENDIF
128         IF( ln_trabbl ) THEN
129            IF( nn_bbl_ldf == 1 ) THEN
130               ahu_bbl_tm(:,:)     = ahu_bbl_tm(:,:)        + ahu_bbl(:,:) 
131               ahv_bbl_tm(:,:)     = ahv_bbl_tm(:,:)        + ahv_bbl(:,:) 
132            ENDIF
133            IF( nn_bbl_adv == 1 ) THEN
134               utr_bbl_tm(:,:)     = utr_bbl_tm(:,:)        + utr_bbl(:,:) 
135               vtr_bbl_tm(:,:)     = vtr_bbl_tm(:,:)        + vtr_bbl(:,:) 
136            ENDIF
137         ENDIF 
138         !
139         sshn_tm  (:,:)         = sshn_tm  (:,:)         + ssh   (:,:,Kmm) 
140         rnf_tm   (:,:)         = rnf_tm   (:,:)         + rnf   (:,:) 
141         h_rnf_tm (:,:)         = h_rnf_tm (:,:)         + h_rnf (:,:) 
142         hmld_tm  (:,:)         = hmld_tm  (:,:)         + hmld  (:,:)
143         fr_i_tm  (:,:)         = fr_i_tm  (:,:)         + fr_i  (:,:)
144         emp_tm   (:,:)         = emp_tm   (:,:)         + emp   (:,:) 
145         fmmflx_tm(:,:)         = fmmflx_tm(:,:)         + fmmflx(:,:)
146         qsr_tm   (:,:)         = qsr_tm   (:,:)         + qsr   (:,:)
147         wndm_tm  (:,:)         = wndm_tm  (:,:)         + wndm  (:,:)
148         !
149      ELSE                           !  It is time to substep
150         !   1. set temporary arrays to hold physics/dynamical variables
151         un_temp    (:,:,:)      = uu    (:,:,:,Kmm)
152         vn_temp    (:,:,:)      = vv    (:,:,:,Kmm)
153         wn_temp    (:,:,:)      = ww    (:,:,:)
154         tsn_temp   (:,:,:,:)    = ts    (:,:,:,:,Kmm)
155         rhop_temp  (:,:,:)      = rhop  (:,:,:)   
156         avs_temp   (:,:,:)      = avs   (:,:,:)
157         IF( l_ldfslp ) THEN
158            uslp_temp  (:,:,:)   = uslp  (:,:,:)   ;   wslpi_temp (:,:,:)   = wslpi (:,:,:)
159            vslp_temp  (:,:,:)   = vslp  (:,:,:)   ;   wslpj_temp (:,:,:)   = wslpj (:,:,:)
160         ENDIF
161         IF( ln_trabbl ) THEN
162            IF( nn_bbl_ldf == 1 ) THEN
163               ahu_bbl_temp(:,:)   = ahu_bbl(:,:) 
164               ahv_bbl_temp(:,:)   = ahv_bbl(:,:) 
165            ENDIF
166            IF( nn_bbl_adv == 1 ) THEN
167               utr_bbl_temp(:,:)   = utr_bbl(:,:) 
168               vtr_bbl_temp(:,:)   = vtr_bbl(:,:) 
169            ENDIF
170         ENDIF
171         sshn_temp  (:,:)        = ssh   (:,:,Kmm)
172         sshb_temp  (:,:)        = ssh   (:,:,Kbb)
173         ssha_temp  (:,:)        = ssh   (:,:,Krhs)
174         rnf_temp   (:,:)        = rnf   (:,:)
175         h_rnf_temp (:,:)        = h_rnf (:,:)
176         hmld_temp  (:,:)        = hmld  (:,:)
177         fr_i_temp  (:,:)        = fr_i  (:,:)
178         emp_temp   (:,:)        = emp   (:,:)
179         emp_b_temp (:,:)        = emp_b (:,:)
180         fmmflx_temp(:,:)        = fmmflx(:,:)
181         qsr_temp   (:,:)        = qsr   (:,:)
182         wndm_temp  (:,:)        = wndm  (:,:)
183         !                                    !  Variables reset in trc_sub_ssh
184         hdivn_temp (:,:,:)      = hdiv (:,:,:)
185         !
186         ! 2. Create averages and reassign variables
187         un_tm    (:,:,:)        = un_tm   (:,:,:)        + uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
188         vn_tm    (:,:,:)        = vn_tm   (:,:,:)        + vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
189         tsn_tm   (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) 
190         tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) 
191         rhop_tm  (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)            * e3t(:,:,:,Kmm) 
192         avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)            * e3w(:,:,:,Kmm) 
193         IF( l_ldfslp ) THEN
194            uslp_tm  (:,:,:)     = uslp_tm (:,:,:)        + uslp (:,:,:) 
195            vslp_tm  (:,:,:)     = vslp_tm (:,:,:)        + vslp (:,:,:)
196            wslpi_tm (:,:,:)     = wslpi_tm(:,:,:)        + wslpi(:,:,:) 
197            wslpj_tm (:,:,:)     = wslpj_tm(:,:,:)        + wslpj(:,:,:) 
198         ENDIF
199         IF( ln_trabbl ) THEN
200            IF( nn_bbl_ldf == 1 ) THEN
201               ahu_bbl_tm(:,:)     = ahu_bbl_tm(:,:)        + ahu_bbl(:,:) 
202               ahv_bbl_tm(:,:)     = ahv_bbl_tm(:,:)        + ahv_bbl(:,:) 
203            ENDIF
204            IF( nn_bbl_adv == 1 ) THEN
205               utr_bbl_tm(:,:)     = utr_bbl_tm(:,:)        + utr_bbl(:,:) 
206               vtr_bbl_tm(:,:)     = vtr_bbl_tm(:,:)        + vtr_bbl(:,:) 
207            ENDIF
208         ENDIF
209         sshn_tm  (:,:)          = sshn_tm    (:,:)       + ssh   (:,:,Kmm) 
210         rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:) 
211         h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:) 
212         hmld_tm  (:,:)          = hmld_tm    (:,:)       + hmld  (:,:)
213         fr_i_tm  (:,:)          = fr_i_tm    (:,:)       + fr_i  (:,:)
214         emp_tm   (:,:)          = emp_tm     (:,:)       + emp   (:,:) 
215         fmmflx_tm(:,:)          = fmmflx_tm  (:,:)       + fmmflx(:,:)
216         qsr_tm   (:,:)          = qsr_tm     (:,:)       + qsr   (:,:)
217         wndm_tm  (:,:)          = wndm_tm    (:,:)       + wndm  (:,:)
218         !
219         ssh      (:,:,Kmm)      = sshn_tm    (:,:) * r1_ndttrcp1 
220         ssh      (:,:,Kbb)      = sshb_hold  (:,:)
221         rnf      (:,:)          = rnf_tm     (:,:) * r1_ndttrcp1 
222         h_rnf    (:,:)          = h_rnf_tm   (:,:) * r1_ndttrcp1 
223         hmld     (:,:)          = hmld_tm    (:,:) * r1_ndttrcp1 
224         !  variables that are initialized after averages
225         emp_b    (:,:)          = emp_b_hold (:,:)
226         IF( kt == nittrc000 ) THEN
227            wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrc 
228            qsr   (:,:)          = qsr_tm     (:,:) * r1_ndttrc 
229            emp   (:,:)          = emp_tm     (:,:) * r1_ndttrc 
230            fmmflx(:,:)          = fmmflx_tm  (:,:) * r1_ndttrc 
231            fr_i  (:,:)          = fr_i_tm    (:,:) * r1_ndttrc
232            IF( ln_trabbl ) THEN
233               IF( nn_bbl_ldf == 1 ) THEN
234                  ahu_bbl(:,:)      = ahu_bbl_tm (:,:) * r1_ndttrc 
235                  ahv_bbl(:,:)      = ahv_bbl_tm (:,:) * r1_ndttrc 
236               ENDIF
237               IF( nn_bbl_adv == 1 ) THEN
238                  utr_bbl(:,:)      = utr_bbl_tm (:,:) * r1_ndttrc 
239                  vtr_bbl(:,:)      = vtr_bbl_tm (:,:) * r1_ndttrc 
240               ENDIF
241            ENDIF
242         ELSE
243            wndm  (:,:)          = wndm_tm    (:,:) * r1_ndttrcp1 
244            qsr   (:,:)          = qsr_tm     (:,:) * r1_ndttrcp1 
245            emp   (:,:)          = emp_tm     (:,:) * r1_ndttrcp1 
246            fmmflx(:,:)          = fmmflx_tm  (:,:) * r1_ndttrcp1 
247            fr_i  (:,:)          = fr_i_tm    (:,:) * r1_ndttrcp1 
248            IF( ln_trabbl ) THEN
249               IF( nn_bbl_ldf == 1 ) THEN
250                  ahu_bbl(:,:)      = ahu_bbl_tm (:,:) * r1_ndttrcp1 
251                  ahv_bbl(:,:)      = ahv_bbl_tm (:,:) * r1_ndttrcp1 
252               ENDIF
253               IF( nn_bbl_adv == 1 ) THEN
254                  utr_bbl(:,:)      = utr_bbl_tm (:,:) * r1_ndttrcp1 
255                  vtr_bbl(:,:)      = vtr_bbl_tm (:,:) * r1_ndttrcp1 
256               ENDIF
257            ENDIF
258         ENDIF
259         !
260         DO jk = 1, jpk
261            DO jj = 1, jpj
262               DO ji = 1, jpi
263                  z1_ne3t = r1_ndttrcp1  / e3t(ji,jj,jk,Kmm)
264                  z1_ne3u = r1_ndttrcp1  / e3u(ji,jj,jk,Kmm)
265                  z1_ne3v = r1_ndttrcp1  / e3v(ji,jj,jk,Kmm)
266                  z1_ne3w = r1_ndttrcp1  / e3w(ji,jj,jk,Kmm)
267                  !
268                  uu   (ji,jj,jk,Kmm)        = un_tm   (ji,jj,jk)        * z1_ne3u
269                  vv   (ji,jj,jk,Kmm)        = vn_tm   (ji,jj,jk)        * z1_ne3v
270                  ts   (ji,jj,jk,jp_tem,Kmm) = tsn_tm  (ji,jj,jk,jp_tem) * z1_ne3t
271                  ts   (ji,jj,jk,jp_sal,Kmm) = tsn_tm  (ji,jj,jk,jp_sal) * z1_ne3t
272                  rhop (ji,jj,jk)            = rhop_tm (ji,jj,jk)        * z1_ne3t
273!!gm : BUG ==>> for avs I don't understand the division by e3w
274                  avs  (ji,jj,jk)            = avs_tm  (ji,jj,jk)        * z1_ne3w
275               END DO
276            END DO
277         END DO
278         IF( l_ldfslp ) THEN
279            wslpi(:,:,:)        = wslpi_tm(:,:,:) 
280            wslpj(:,:,:)        = wslpj_tm(:,:,:)
281            uslp (:,:,:)        = uslp_tm (:,:,:)
282            vslp (:,:,:)        = vslp_tm (:,:,:)
283         ENDIF
284         !
285         CALL trc_sub_ssh( kt, Kbb, Kmm, Krhs )         ! after ssh & vertical velocity
286         !
287      ENDIF
288      !
289      IF( ln_timing )   CALL timing_stop('trc_sub_stp')
290      !
291   END SUBROUTINE trc_sub_stp
292
293
294   SUBROUTINE trc_sub_ini( Kmm )
295      !!-------------------------------------------------------------------
296      !!                     ***  ROUTINE trc_sub_ini  ***
297      !!                     
298      !! ** Purpose : Initialize variables needed for sub-stepping passive tracers
299      !!
300      !! ** Method  :
301      !!              Compute the averages for sub-stepping
302      !!-------------------------------------------------------------------
303      INTEGER ::   ierr
304      INTEGER, INTENT( in ) ::   Kmm    ! time level index
305      !!-------------------------------------------------------------------
306      !
307      IF(lwp) WRITE(numout,*)
308      IF(lwp) WRITE(numout,*) 'trc_sub_ini : initial set up of the passive tracers substepping'
309      IF(lwp) WRITE(numout,*) '~~~~~~~'
310
311      ierr =  trc_sub_alloc    ()
312      CALL mpp_sum( 'trcsub', ierr )
313      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
314
315      un_tm   (:,:,:)        = uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
316      vn_tm   (:,:,:)        = vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
317      tsn_tm  (:,:,:,jp_tem) = ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) 
318      tsn_tm  (:,:,:,jp_sal) = ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) 
319      rhop_tm (:,:,:)        = rhop (:,:,:)            * e3t(:,:,:,Kmm) 
320!!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
321      avs_tm  (:,:,:)        = avs  (:,:,:)            * e3w(:,:,:,Kmm) 
322      IF( l_ldfslp ) THEN
323         wslpi_tm(:,:,:)     = wslpi(:,:,:)
324         wslpj_tm(:,:,:)     = wslpj(:,:,:)
325         uslp_tm (:,:,:)     = uslp (:,:,:)
326         vslp_tm (:,:,:)     = vslp (:,:,:)
327      ENDIF
328      sshn_tm  (:,:) = ssh   (:,:,Kmm) 
329      rnf_tm   (:,:) = rnf   (:,:) 
330      h_rnf_tm (:,:) = h_rnf (:,:) 
331      hmld_tm  (:,:) = hmld  (:,:)
332
333      ! Physics variables that are set after initialization:
334      fr_i_tm  (:,:) = 0._wp
335      emp_tm   (:,:) = 0._wp
336      fmmflx_tm(:,:) = 0._wp
337      qsr_tm   (:,:) = 0._wp
338      wndm_tm  (:,:) = 0._wp
339      IF( ln_trabbl ) THEN
340         IF( nn_bbl_ldf == 1 ) THEN
341            ahu_bbl_tm(:,:) = 0._wp
342            ahv_bbl_tm(:,:) = 0._wp
343         ENDIF
344         IF( nn_bbl_adv == 1 ) THEN
345            utr_bbl_tm(:,:) = 0._wp
346            vtr_bbl_tm(:,:) = 0._wp
347         ENDIF
348      ENDIF
349      !
350   END SUBROUTINE trc_sub_ini
351
352
353   SUBROUTINE trc_sub_reset( kt, Kbb, Kmm, Krhs )
354      !!-------------------------------------------------------------------
355      !!                     ***  ROUTINE trc_sub_reset  ***
356      !!                     
357      !! ** Purpose : Reset physics variables averaged for substepping
358      !!
359      !! ** Method  :
360      !!              Compute the averages for sub-stepping
361      !!-------------------------------------------------------------------
362      INTEGER, INTENT( in ) ::  kt             ! ocean time-step index
363      INTEGER, INTENT( in ) ::  Kbb, Kmm, Krhs ! time level index
364      INTEGER :: jk                            ! dummy loop indices
365      !!-------------------------------------------------------------------
366      !
367      IF( ln_timing )   CALL timing_start('trc_sub_reset')
368      !
369      !   restore physics variables
370      uu   (:,:,:,Kmm)   =  un_temp  (:,:,:)
371      vv   (:,:,:,Kmm)   =  vn_temp  (:,:,:)
372      ww   (:,:,:)       =  wn_temp  (:,:,:)
373      ts   (:,:,:,:,Kmm) =  tsn_temp (:,:,:,:)
374      rhop (:,:,:)       =  rhop_temp(:,:,:)
375      avs  (:,:,:)       =  avs_temp (:,:,:)
376      IF( l_ldfslp ) THEN
377         wslpi (:,:,:)   =  wslpi_temp (:,:,:)
378         wslpj (:,:,:)   =  wslpj_temp (:,:,:)
379         uslp  (:,:,:)   =  uslp_temp  (:,:,:)
380         vslp  (:,:,:)   =  vslp_temp  (:,:,:)
381      ENDIF
382      ssh   (:,:,Kmm)     =  sshn_temp  (:,:)
383      ssh   (:,:,Kbb)     =  sshb_temp  (:,:)
384      ssh   (:,:,Krhs)    =  ssha_temp  (:,:)
385      rnf   (:,:)         =  rnf_temp   (:,:)
386      h_rnf (:,:)         =  h_rnf_temp (:,:)
387      !
388      hmld  (:,:)     =  hmld_temp  (:,:)
389      fr_i  (:,:)     =  fr_i_temp  (:,:)
390      emp   (:,:)     =  emp_temp   (:,:)
391      fmmflx(:,:)     =  fmmflx_temp(:,:)
392      emp_b (:,:)     =  emp_b_temp (:,:)
393      qsr   (:,:)     =  qsr_temp   (:,:)
394      wndm  (:,:)     =  wndm_temp  (:,:)
395      IF( ln_trabbl ) THEN
396         IF( nn_bbl_ldf == 1 ) THEN
397            ahu_bbl(:,:) = ahu_bbl_temp(:,:) 
398            ahv_bbl(:,:) = ahv_bbl_temp(:,:) 
399         ENDIF
400         IF( nn_bbl_adv == 1 ) THEN
401            utr_bbl(:,:) = utr_bbl_temp(:,:) 
402            vtr_bbl(:,:) = vtr_bbl_temp(:,:) 
403         ENDIF
404      ENDIF
405      !
406      hdiv (:,:,:)   =  hdivn_temp (:,:,:)
407      !                                     
408      ! Start new averages
409         un_tm   (:,:,:)        = uu   (:,:,:,Kmm)        * e3u(:,:,:,Kmm) 
410         vn_tm   (:,:,:)        = vv   (:,:,:,Kmm)        * e3v(:,:,:,Kmm) 
411         tsn_tm  (:,:,:,jp_tem) = ts   (:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) 
412         tsn_tm  (:,:,:,jp_sal) = ts   (:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) 
413         rhop_tm (:,:,:)        = rhop (:,:,:)            * e3t(:,:,:,Kmm) 
414         avs_tm  (:,:,:)        = avs  (:,:,:)            * e3w(:,:,:,Kmm) 
415      IF( l_ldfslp ) THEN
416         uslp_tm (:,:,:)        = uslp (:,:,:)
417         vslp_tm (:,:,:)        = vslp (:,:,:)
418         wslpi_tm(:,:,:)        = wslpi(:,:,:) 
419         wslpj_tm(:,:,:)        = wslpj(:,:,:)
420      ENDIF
421      !
422      sshb_hold  (:,:) = ssh   (:,:,Kmm )
423      emp_b_hold (:,:) = emp   (:,:)
424      sshn_tm    (:,:) = ssh   (:,:,Kmm) 
425      rnf_tm     (:,:) = rnf   (:,:) 
426      h_rnf_tm   (:,:) = h_rnf (:,:) 
427      hmld_tm    (:,:) = hmld  (:,:)
428      fr_i_tm    (:,:) = fr_i  (:,:)
429      emp_tm     (:,:) = emp   (:,:)
430      fmmflx_tm  (:,:) = fmmflx(:,:)
431      qsr_tm     (:,:) = qsr   (:,:)
432      wndm_tm    (:,:) = wndm  (:,:)
433      IF( ln_trabbl ) THEN
434         IF( nn_bbl_ldf == 1 ) THEN
435            ahu_bbl_tm(:,:) = ahu_bbl(:,:) 
436            ahv_bbl_tm(:,:) = ahv_bbl(:,:) 
437         ENDIF
438         IF( nn_bbl_adv == 1 ) THEN
439            utr_bbl_tm(:,:) = utr_bbl(:,:) 
440            vtr_bbl_tm(:,:) = vtr_bbl(:,:) 
441         ENDIF
442      ENDIF
443      !
444      !
445      IF( ln_timing )   CALL timing_stop('trc_sub_reset')
446      !
447   END SUBROUTINE trc_sub_reset
448
449
450   SUBROUTINE trc_sub_ssh( kt, Kbb, Kmm, Krhs ) 
451      !!----------------------------------------------------------------------
452      !!                ***  ROUTINE trc_sub_ssh  ***
453      !!                   
454      !! ** Purpose :   compute the after ssh (ssh(:,:,Krhs)), the now vertical velocity
455      !!              and update the now vertical coordinate (ln_linssh=F).
456      !!
457      !! ** Method  : - Using the incompressibility hypothesis, the vertical
458      !!      velocity is computed by integrating the horizontal divergence 
459      !!      from the bottom to the surface minus the scale factor evolution.
460      !!        The boundary conditions are w=0 at the bottom (no flux) and.
461      !!
462      !! ** action  :   ssh(:,:,Krhs)    : after sea surface height
463      !!                ww      : now vertical velocity
464      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (ln_linssh=F)
465      !!
466      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
467      !!----------------------------------------------------------------------
468      INTEGER, INTENT(in) ::   kt              ! time step
469      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! ocean time-level index
470      !
471      INTEGER  ::   ji, jj, jk   ! dummy loop indices
472      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
473      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
474      !!---------------------------------------------------------------------
475      !
476      IF( ln_timing )   CALL timing_start('trc_sub_ssh')
477      !
478
479      IF( kt == nittrc000 ) THEN
480         !
481         IF(lwp) WRITE(numout,*)
482         IF(lwp) WRITE(numout,*) 'trc_sub_ssh : after sea surface height and now vertical velocity '
483         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
484         !
485         ww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
486         !
487      ENDIF
488      !
489!!gm BUG here !   hdiv will include the runoff divergence at the wrong timestep !!!!
490      CALL div_hor( kt, Kmm )                         ! Horizontal divergence & Relative vorticity
491      !
492      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
493      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
494
495      !                                           !------------------------------!
496      !                                           !   After Sea Surface Height   !
497      !                                           !------------------------------!
498      zhdiv(:,:) = 0._wp
499      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
500        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
501      END DO
502      !                                                ! Sea surface elevation time stepping
503      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
504      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
505      z1_rau0 = 0.5 / rau0
506      ssh(:,:,Krhs) = (  ssh(:,:,Kbb) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
507
508      IF( .NOT.ln_dynspg_ts ) THEN
509      ! These lines are not necessary with time splitting since
510      ! boundary condition on sea level is set during ts loop
511#if defined key_agrif
512      CALL agrif_ssh( kt )
513#endif
514         IF( ln_bdy ) THEN
515            ssh(:,:,Krhs) = ssh(:,:,Krhs) * bdytmask(:,:)
516            CALL lbc_lnk( 'trcsub', ssh(:,:,Krhs), 'T', 1. ) 
517         ENDIF
518      ENDIF
519      !
520      !                                           !------------------------------!
521      !                                           !     Now Vertical Velocity    !
522      !                                           !------------------------------!
523      z1_2dt = 1.e0 / z2dt
524      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
525         ! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise
526         ww(:,:,jk) = ww(:,:,jk+1) -   e3t(:,:,jk,Kmm ) * hdiv(:,:,jk)          &
527            &                      - ( e3t(:,:,jk,Krhs) - e3t (:,:,jk,Kbb) )    &
528            &                         * tmask(:,:,jk) * z1_2dt
529         IF( ln_bdy ) ww(:,:,jk) = ww(:,:,jk) * bdytmask(:,:)
530      END DO
531      !
532      IF( ln_timing )   CALL timing_stop('trc_sub_ssh')
533      !
534   END SUBROUTINE trc_sub_ssh
535
536
537   INTEGER FUNCTION trc_sub_alloc()
538      !!-------------------------------------------------------------------
539      !!                    *** ROUTINE trc_sub_alloc ***
540      !!-------------------------------------------------------------------
541      USE lib_mpp, ONLY: ctl_stop
542      INTEGER ::  ierr(3)
543      !!-------------------------------------------------------------------
544      !
545      ierr(:) = 0
546      !
547      ALLOCATE( un_temp(jpi,jpj,jpk)      ,  vn_temp(jpi,jpj,jpk)   ,     &
548         &      wn_temp(jpi,jpj,jpk)      ,                               &
549         &      rhop_temp(jpi,jpj,jpk)    ,  rhop_tm(jpi,jpj,jpk)   ,     &
550         &      sshn_temp(jpi,jpj)        ,  sshb_temp(jpi,jpj)     ,     &
551         &      ssha_temp(jpi,jpj)        ,                               &
552         &      rnf_temp(jpi,jpj)         ,  h_rnf_temp(jpi,jpj)    ,     &
553         &      tsn_temp(jpi,jpj,jpk,2)   ,  emp_b_temp(jpi,jpj)    ,     &
554         &      emp_temp(jpi,jpj)         ,  fmmflx_temp(jpi,jpj)   ,     &
555         &      hmld_temp(jpi,jpj)        ,  qsr_temp(jpi,jpj)      ,     &
556         &      fr_i_temp(jpi,jpj)        ,  fr_i_tm(jpi,jpj)       ,     &
557         &      wndm_temp(jpi,jpj)        ,  wndm_tm(jpi,jpj)       ,     &
558         &      avs_tm(jpi,jpj,jpk)       ,  avs_temp(jpi,jpj,jpk)  ,     &
559         &      hdivn_temp(jpi,jpj,jpk)   ,  hdivb_temp(jpi,jpj,jpk),     &
560         &      un_tm(jpi,jpj,jpk)        ,  vn_tm(jpi,jpj,jpk)     ,     &
561         &      sshn_tm(jpi,jpj)          ,  sshb_hold(jpi,jpj)     ,     &
562         &      tsn_tm(jpi,jpj,jpk,2)     ,                               &
563         &      emp_tm(jpi,jpj)           ,  fmmflx_tm(jpi,jpj)     ,     &
564         &      emp_b_hold(jpi,jpj)       ,                               &
565         &      hmld_tm(jpi,jpj)          ,  qsr_tm(jpi,jpj)        ,     &
566         &      rnf_tm(jpi,jpj)           ,  h_rnf_tm(jpi,jpj)      , STAT=ierr(1) ) 
567      !
568      IF( l_ldfslp ) THEN
569         ALLOCATE( uslp_temp(jpi,jpj,jpk) ,  wslpi_temp(jpi,jpj,jpk),     &
570            &      vslp_temp(jpi,jpj,jpk) ,  wslpj_temp(jpi,jpj,jpk),     &
571            &      uslp_tm  (jpi,jpj,jpk) ,  wslpi_tm  (jpi,jpj,jpk),     &
572            &      vslp_tm  (jpi,jpj,jpk) ,  wslpj_tm  (jpi,jpj,jpk), STAT=ierr(2) )
573      ENDIF
574      IF( ln_trabbl ) THEN
575         ALLOCATE( ahu_bbl_temp(jpi,jpj)  , utr_bbl_temp(jpi,jpj)   ,     &
576            &      ahv_bbl_temp(jpi,jpj)  , vtr_bbl_temp(jpi,jpj)   ,     &
577            &      ahu_bbl_tm  (jpi,jpj)  , utr_bbl_tm  (jpi,jpj)   ,     &
578            &      ahv_bbl_tm  (jpi,jpj)  , vtr_bbl_tm  (jpi,jpj)   , STAT=ierr(3) ) 
579      ENDIF
580      !
581      trc_sub_alloc = MAXVAL( ierr )
582      !
583      IF( trc_sub_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trc_sub_alloc: failed to allocate arrays' )
584      !
585   END FUNCTION trc_sub_alloc
586
587#else
588   !!----------------------------------------------------------------------
589   !!   Default key                                     NO passive tracers
590   !!----------------------------------------------------------------------
591CONTAINS
592   SUBROUTINE trc_sub_stp( kt, Kbb, Kmm, Krhs )        ! Empty routine
593      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices
594      WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
595   END SUBROUTINE trc_sub_stp
596   SUBROUTINE trc_sub_ini( Kmm )        ! Empty routine
597      INTEGER, INTENT( in ) ::   Kmm    ! time level index
598      WRITE(*,*) 'trc_sub_ini: You should not have seen this print! error?', kt
599   END SUBROUTINE trc_sub_ini
600#endif
601
602   !!======================================================================
603END MODULE trcsub
Note: See TracBrowser for help on using the repository browser.