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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP – NEMO

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

Last change on this file since 10922 was 10922, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert IOM, LDF, OBS and SBC directories and compatibility changes elsewhere that these changes enforce. Changes pass SETTE and compare with original trunk results. Outstanding issues (currently with work-arounds) in DIU/step_diu.F90 and fld_bdy_interp within SBC/fldread.F90; proper soltions pending
  • Property svn:keywords set to Id
File size: 29.1 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, Kmm )
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 ) ::   Kmm  ! 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  (:,:) = sshn  (:,:)
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   (:,:,:)        + un   (:,:,:)        * e3u_n(:,:,:) 
117         vn_tm   (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_n(:,:,:) 
118         tsn_tm  (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
119         tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
120         rhop_tm (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:) 
121         avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:) 
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  (:,:)         + sshn  (:,:) 
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    (:,:,:)      = un    (:,:,:)
152         vn_temp    (:,:,:)      = vn    (:,:,:)
153         wn_temp    (:,:,:)      = wn    (:,:,:)
154         tsn_temp   (:,:,:,:)    = tsn   (:,:,:,:)
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  (:,:)        = sshn  (:,:)
172         sshb_temp  (:,:)        = sshb  (:,:)
173         ssha_temp  (:,:)        = ssha  (:,:)
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 (:,:,:)      = hdivn (:,:,:)
185         !
186         ! 2. Create averages and reassign variables
187         un_tm    (:,:,:)        = un_tm   (:,:,:)        + un   (:,:,:)        * e3u_n(:,:,:) 
188         vn_tm    (:,:,:)        = vn_tm   (:,:,:)        + vn   (:,:,:)        * e3v_n(:,:,:) 
189         tsn_tm   (:,:,:,jp_tem) = tsn_tm  (:,:,:,jp_tem) + tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
190         tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
191         rhop_tm (:,:,:)         = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:) 
192         avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:) 
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    (:,:)       + sshn  (:,:) 
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         sshn     (:,:)          = sshn_tm    (:,:) * r1_ndttrcp1 
220         sshb     (:,:)          = 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_n(ji,jj,jk)
264                  z1_ne3u = r1_ndttrcp1  / e3u_n(ji,jj,jk)
265                  z1_ne3v = r1_ndttrcp1  / e3v_n(ji,jj,jk)
266                  z1_ne3w = r1_ndttrcp1  / e3w_n(ji,jj,jk)
267                  !
268                  un   (ji,jj,jk)        = un_tm   (ji,jj,jk)        * z1_ne3u
269                  vn   (ji,jj,jk)        = vn_tm   (ji,jj,jk)        * z1_ne3v
270                  tsn  (ji,jj,jk,jp_tem) = tsn_tm  (ji,jj,jk,jp_tem) * z1_ne3t
271                  tsn  (ji,jj,jk,jp_sal) = 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, Kmm )         ! 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
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      !!-------------------------------------------------------------------
305      !
306      IF(lwp) WRITE(numout,*)
307      IF(lwp) WRITE(numout,*) 'trc_sub_ini : initial set up of the passive tracers substepping'
308      IF(lwp) WRITE(numout,*) '~~~~~~~'
309
310      ierr =  trc_sub_alloc    ()
311      CALL mpp_sum( 'trcsub', ierr )
312      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
313
314      un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
315      vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
316      tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
317      tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
318      rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:) 
319!!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
320      avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:) 
321      IF( l_ldfslp ) THEN
322         wslpi_tm(:,:,:)     = wslpi(:,:,:)
323         wslpj_tm(:,:,:)     = wslpj(:,:,:)
324         uslp_tm (:,:,:)     = uslp (:,:,:)
325         vslp_tm (:,:,:)     = vslp (:,:,:)
326      ENDIF
327      sshn_tm  (:,:) = sshn  (:,:) 
328      rnf_tm   (:,:) = rnf   (:,:) 
329      h_rnf_tm (:,:) = h_rnf (:,:) 
330      hmld_tm  (:,:) = hmld  (:,:)
331
332      ! Physics variables that are set after initialization:
333      fr_i_tm  (:,:) = 0._wp
334      emp_tm   (:,:) = 0._wp
335      fmmflx_tm(:,:)  = 0._wp
336      qsr_tm   (:,:) = 0._wp
337      wndm_tm  (:,:) = 0._wp
338      IF( ln_trabbl ) THEN
339         IF( nn_bbl_ldf == 1 ) THEN
340            ahu_bbl_tm(:,:) = 0._wp
341            ahv_bbl_tm(:,:) = 0._wp
342         ENDIF
343         IF( nn_bbl_adv == 1 ) THEN
344            utr_bbl_tm(:,:) = 0._wp
345            vtr_bbl_tm(:,:) = 0._wp
346         ENDIF
347      ENDIF
348      !
349   END SUBROUTINE trc_sub_ini
350
351
352   SUBROUTINE trc_sub_reset( kt )
353      !!-------------------------------------------------------------------
354      !!                     ***  ROUTINE trc_sub_reset  ***
355      !!                     
356      !! ** Purpose : Reset physics variables averaged for substepping
357      !!
358      !! ** Method  :
359      !!              Compute the averages for sub-stepping
360      !!-------------------------------------------------------------------
361      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
362      INTEGER :: jk                 ! dummy loop indices
363      !!-------------------------------------------------------------------
364      !
365      IF( ln_timing )   CALL timing_start('trc_sub_reset')
366      !
367      !   restore physics variables
368      un    (:,:,:)   =  un_temp    (:,:,:)
369      vn    (:,:,:)   =  vn_temp    (:,:,:)
370      wn    (:,:,:)   =  wn_temp    (:,:,:)
371      tsn   (:,:,:,:) =  tsn_temp   (:,:,:,:)
372      rhop  (:,:,:)   =  rhop_temp  (:,:,:)
373      avs   (:,:,:)   =  avs_temp   (:,:,:)
374      IF( l_ldfslp ) THEN
375         wslpi (:,:,:)=  wslpi_temp (:,:,:)
376         wslpj (:,:,:)=  wslpj_temp (:,:,:)
377         uslp  (:,:,:)=  uslp_temp  (:,:,:)
378         vslp  (:,:,:)=  vslp_temp  (:,:,:)
379      ENDIF
380      sshn  (:,:)     =  sshn_temp  (:,:)
381      sshb  (:,:)     =  sshb_temp  (:,:)
382      ssha  (:,:)     =  ssha_temp  (:,:)
383      rnf   (:,:)     =  rnf_temp   (:,:)
384      h_rnf (:,:)     =  h_rnf_temp (:,:)
385      !
386      hmld  (:,:)     =  hmld_temp  (:,:)
387      fr_i  (:,:)     =  fr_i_temp  (:,:)
388      emp   (:,:)     =  emp_temp   (:,:)
389      fmmflx(:,:)     =  fmmflx_temp(:,:)
390      emp_b (:,:)     =  emp_b_temp (:,:)
391      qsr   (:,:)     =  qsr_temp   (:,:)
392      wndm  (:,:)     =  wndm_temp  (:,:)
393      IF( ln_trabbl ) THEN
394         IF( nn_bbl_ldf == 1 ) THEN
395            ahu_bbl(:,:) = ahu_bbl_temp(:,:) 
396            ahv_bbl(:,:) = ahv_bbl_temp(:,:) 
397         ENDIF
398         IF( nn_bbl_adv == 1 ) THEN
399            utr_bbl(:,:) = utr_bbl_temp(:,:) 
400            vtr_bbl(:,:) = vtr_bbl_temp(:,:) 
401         ENDIF
402      ENDIF
403      !
404      hdivn (:,:,:)   =  hdivn_temp (:,:,:)
405      !                                     
406      ! Start new averages
407         un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
408         vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
409         tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
410         tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
411         rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:) 
412         avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:) 
413      IF( l_ldfslp ) THEN
414         uslp_tm (:,:,:)        = uslp (:,:,:)
415         vslp_tm (:,:,:)        = vslp (:,:,:)
416         wslpi_tm(:,:,:)        = wslpi(:,:,:) 
417         wslpj_tm(:,:,:)        = wslpj(:,:,:)
418      ENDIF
419      !
420      sshb_hold  (:,:) = sshn  (:,:)
421      emp_b_hold (:,:) = emp   (:,:)
422      sshn_tm    (:,:) = sshn  (:,:) 
423      rnf_tm     (:,:) = rnf   (:,:) 
424      h_rnf_tm   (:,:) = h_rnf (:,:) 
425      hmld_tm    (:,:) = hmld  (:,:)
426      fr_i_tm    (:,:) = fr_i  (:,:)
427      emp_tm     (:,:) = emp   (:,:)
428      fmmflx_tm  (:,:) = fmmflx(:,:)
429      qsr_tm     (:,:) = qsr   (:,:)
430      wndm_tm    (:,:) = wndm  (:,:)
431      IF( ln_trabbl ) THEN
432         IF( nn_bbl_ldf == 1 ) THEN
433            ahu_bbl_tm(:,:) = ahu_bbl(:,:) 
434            ahv_bbl_tm(:,:) = ahv_bbl(:,:) 
435         ENDIF
436         IF( nn_bbl_adv == 1 ) THEN
437            utr_bbl_tm(:,:) = utr_bbl(:,:) 
438            vtr_bbl_tm(:,:) = vtr_bbl(:,:) 
439         ENDIF
440      ENDIF
441      !
442      !
443      IF( ln_timing )   CALL timing_stop('trc_sub_reset')
444      !
445   END SUBROUTINE trc_sub_reset
446
447
448   SUBROUTINE trc_sub_ssh( kt, Kmm ) 
449      !!----------------------------------------------------------------------
450      !!                ***  ROUTINE trc_sub_ssh  ***
451      !!                   
452      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
453      !!              and update the now vertical coordinate (ln_linssh=F).
454      !!
455      !! ** Method  : - Using the incompressibility hypothesis, the vertical
456      !!      velocity is computed by integrating the horizontal divergence 
457      !!      from the bottom to the surface minus the scale factor evolution.
458      !!        The boundary conditions are w=0 at the bottom (no flux) and.
459      !!
460      !! ** action  :   ssha    : after sea surface height
461      !!                wn      : now vertical velocity
462      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (ln_linssh=F)
463      !!
464      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
465      !!----------------------------------------------------------------------
466      INTEGER, INTENT(in) ::   kt   ! time step
467      INTEGER, INTENT(in) ::   Kmm  ! ocean time-level index
468      !
469      INTEGER  ::   ji, jj, jk   ! dummy loop indices
470      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
471      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
472      !!---------------------------------------------------------------------
473      !
474      IF( ln_timing )   CALL timing_start('trc_sub_ssh')
475      !
476
477      IF( kt == nittrc000 ) THEN
478         !
479         IF(lwp) WRITE(numout,*)
480         IF(lwp) WRITE(numout,*) 'trc_sub_ssh : after sea surface height and now vertical velocity '
481         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
482         !
483         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
484         !
485      ENDIF
486      !
487!!gm BUG here !   hdivn will include the runoff divergence at the wrong timestep !!!!
488      CALL div_hor( kt, Kmm )                         ! Horizontal divergence & Relative vorticity
489      !
490      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
491      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
492
493      !                                           !------------------------------!
494      !                                           !   After Sea Surface Height   !
495      !                                           !------------------------------!
496      zhdiv(:,:) = 0._wp
497      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
498        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
499      END DO
500      !                                                ! Sea surface elevation time stepping
501      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
502      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
503      z1_rau0 = 0.5 / rau0
504      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
505
506      IF( .NOT.ln_dynspg_ts ) THEN
507      ! These lines are not necessary with time splitting since
508      ! boundary condition on sea level is set during ts loop
509#if defined key_agrif
510      CALL agrif_ssh( kt )
511#endif
512         IF( ln_bdy ) THEN
513            ssha(:,:) = ssha(:,:) * bdytmask(:,:)
514            CALL lbc_lnk( 'trcsub', ssha, 'T', 1. ) 
515         ENDIF
516      ENDIF
517      !
518      !                                           !------------------------------!
519      !                                           !     Now Vertical Velocity    !
520      !                                           !------------------------------!
521      z1_2dt = 1.e0 / z2dt
522      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
523         ! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise
524         wn(:,:,jk) = wn(:,:,jk+1) -   e3t_n(:,:,jk) * hdivn(:,:,jk)        &
525            &                      - ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )    &
526            &                         * tmask(:,:,jk) * z1_2dt
527         IF( ln_bdy ) wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
528      END DO
529      !
530      IF( ln_timing )   CALL timing_stop('trc_sub_ssh')
531      !
532   END SUBROUTINE trc_sub_ssh
533
534
535   INTEGER FUNCTION trc_sub_alloc()
536      !!-------------------------------------------------------------------
537      !!                    *** ROUTINE trc_sub_alloc ***
538      !!-------------------------------------------------------------------
539      USE lib_mpp, ONLY: ctl_stop
540      INTEGER ::  ierr(3)
541      !!-------------------------------------------------------------------
542      !
543      ierr(:) = 0
544      !
545      ALLOCATE( un_temp(jpi,jpj,jpk)      ,  vn_temp(jpi,jpj,jpk)   ,     &
546         &      wn_temp(jpi,jpj,jpk)      ,                               &
547         &      rhop_temp(jpi,jpj,jpk)    ,  rhop_tm(jpi,jpj,jpk)   ,     &
548         &      sshn_temp(jpi,jpj)        ,  sshb_temp(jpi,jpj)     ,     &
549         &      ssha_temp(jpi,jpj)        ,                               &
550         &      rnf_temp(jpi,jpj)         ,  h_rnf_temp(jpi,jpj)    ,     &
551         &      tsn_temp(jpi,jpj,jpk,2)   ,  emp_b_temp(jpi,jpj)    ,     &
552         &      emp_temp(jpi,jpj)         ,  fmmflx_temp(jpi,jpj)   ,     &
553         &      hmld_temp(jpi,jpj)        ,  qsr_temp(jpi,jpj)      ,     &
554         &      fr_i_temp(jpi,jpj)        ,  fr_i_tm(jpi,jpj)       ,     &
555         &      wndm_temp(jpi,jpj)        ,  wndm_tm(jpi,jpj)       ,     &
556         &      avs_tm(jpi,jpj,jpk)       ,  avs_temp(jpi,jpj,jpk)  ,     &
557         &      hdivn_temp(jpi,jpj,jpk)   ,  hdivb_temp(jpi,jpj,jpk),     &
558         &      un_tm(jpi,jpj,jpk)        ,  vn_tm(jpi,jpj,jpk)     ,     &
559         &      sshn_tm(jpi,jpj)          ,  sshb_hold(jpi,jpj)     ,     &
560         &      tsn_tm(jpi,jpj,jpk,2)     ,                               &
561         &      emp_tm(jpi,jpj)           ,  fmmflx_tm(jpi,jpj)     ,     &
562         &      emp_b_hold(jpi,jpj)       ,                               &
563         &      hmld_tm(jpi,jpj)          ,  qsr_tm(jpi,jpj)        ,     &
564         &      rnf_tm(jpi,jpj)           ,  h_rnf_tm(jpi,jpj)      , STAT=ierr(1) ) 
565      !
566      IF( l_ldfslp ) THEN
567         ALLOCATE( uslp_temp(jpi,jpj,jpk) ,  wslpi_temp(jpi,jpj,jpk),     &
568            &      vslp_temp(jpi,jpj,jpk) ,  wslpj_temp(jpi,jpj,jpk),     &
569            &      uslp_tm  (jpi,jpj,jpk) ,  wslpi_tm  (jpi,jpj,jpk),     &
570            &      vslp_tm  (jpi,jpj,jpk) ,  wslpj_tm  (jpi,jpj,jpk), STAT=ierr(2) )
571      ENDIF
572      IF( ln_trabbl ) THEN
573         ALLOCATE( ahu_bbl_temp(jpi,jpj)  , utr_bbl_temp(jpi,jpj)   ,     &
574            &      ahv_bbl_temp(jpi,jpj)  , vtr_bbl_temp(jpi,jpj)   ,     &
575            &      ahu_bbl_tm  (jpi,jpj)  , utr_bbl_tm  (jpi,jpj)   ,     &
576            &      ahv_bbl_tm  (jpi,jpj)  , vtr_bbl_tm  (jpi,jpj)   , STAT=ierr(3) ) 
577      ENDIF
578      !
579      trc_sub_alloc = MAXVAL( ierr )
580      !
581      IF( trc_sub_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trc_sub_alloc: failed to allocate arrays' )
582      !
583   END FUNCTION trc_sub_alloc
584
585#else
586   !!----------------------------------------------------------------------
587   !!   Default key                                     NO passive tracers
588   !!----------------------------------------------------------------------
589CONTAINS
590   SUBROUTINE trc_sub_stp( kt )        ! Empty routine
591      WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
592   END SUBROUTINE trc_sub_stp
593   SUBROUTINE trc_sub_ini        ! Empty routine
594      WRITE(*,*) 'trc_sub_ini: You should not have seen this print! error?', kt
595   END SUBROUTINE trc_sub_ini
596#endif
597
598   !!======================================================================
599END MODULE trcsub
Note: See TracBrowser for help on using the repository browser.