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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/trcsub.F90 @ 9119

Last change on this file since 9119 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch

  • Property svn:keywords set to Id
File size: 29.3 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_opa_update
27   USE agrif_opa_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 (2017)
79   !! $Id$
80   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
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( nn_timing == 1 )  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( nn_timing == 1 )  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( nn_timing == 1 )  CALL timing_start('trc_sub_ini')
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      IF( lk_mpp    )   CALL mpp_sum( ierr )
313      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
314
315      un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
316      vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
317      tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
318      tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
319      rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:) 
320!!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
321      avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:) 
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  (:,:) = sshn  (:,:) 
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      IF( nn_timing == 1 )  CALL timing_stop('trc_sub_ini')
351      !
352   END SUBROUTINE trc_sub_ini
353
354   SUBROUTINE trc_sub_reset( kt )
355      !!-------------------------------------------------------------------
356      !!                     ***  ROUTINE trc_sub_reset  ***
357      !!                     
358      !! ** Purpose : Reset physics variables averaged for substepping
359      !!
360      !! ** Method  :
361      !!              Compute the averages for sub-stepping
362      !!-------------------------------------------------------------------
363      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
364      INTEGER :: jk                 ! dummy loop indices
365      !!-------------------------------------------------------------------
366      !
367      IF( nn_timing == 1 )  CALL timing_start('trc_sub_reset')
368      !
369      !   restore physics variables
370      un    (:,:,:)   =  un_temp    (:,:,:)
371      vn    (:,:,:)   =  vn_temp    (:,:,:)
372      wn    (:,:,:)   =  wn_temp    (:,:,:)
373      tsn   (:,:,:,:) =  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      sshn  (:,:)     =  sshn_temp  (:,:)
383      sshb  (:,:)     =  sshb_temp  (:,:)
384      ssha  (:,:)     =  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      hdivn (:,:,:)   =  hdivn_temp (:,:,:)
407      !                                     
408      ! Start new averages
409         un_tm   (:,:,:)        = un   (:,:,:)        * e3u_n(:,:,:) 
410         vn_tm   (:,:,:)        = vn   (:,:,:)        * e3v_n(:,:,:) 
411         tsn_tm  (:,:,:,jp_tem) = tsn  (:,:,:,jp_tem) * e3t_n(:,:,:) 
412         tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:) 
413         rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:) 
414         avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:) 
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  (:,:) = sshn  (:,:)
423      emp_b_hold (:,:) = emp   (:,:)
424      sshn_tm    (:,:) = sshn  (:,:) 
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( nn_timing == 1 )  CALL timing_stop('trc_sub_reset')
446      !
447   END SUBROUTINE trc_sub_reset
448
449
450   SUBROUTINE trc_sub_ssh( kt ) 
451      !!----------------------------------------------------------------------
452      !!                ***  ROUTINE trc_sub_ssh  ***
453      !!                   
454      !! ** Purpose :   compute the after ssh (ssha), 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  :   ssha    : after sea surface height
463      !!                wn      : 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      !
470      INTEGER  ::   ji, jj, jk   ! dummy loop indices
471      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
472      REAL(wp), POINTER, DIMENSION(:,:) :: zhdiv
473      !!---------------------------------------------------------------------
474      !
475      IF( nn_timing == 1 )  CALL timing_start('trc_sub_ssh')
476      !
477      ! Allocate temporary workspace
478      CALL wrk_alloc( jpi,jpj,   zhdiv )
479
480      IF( kt == nittrc000 ) THEN
481         !
482         IF(lwp) WRITE(numout,*)
483         IF(lwp) WRITE(numout,*) 'trc_sub_ssh : after sea surface height and now vertical velocity '
484         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
485         !
486         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
487         !
488      ENDIF
489      !
490!!gm BUG here !   hdivn will include the runoff divergence at the wrong timestep !!!!
491      CALL div_hor( kt )                              ! Horizontal divergence & Relative vorticity
492      !
493      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
494      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
495
496      !                                           !------------------------------!
497      !                                           !   After Sea Surface Height   !
498      !                                           !------------------------------!
499      zhdiv(:,:) = 0._wp
500      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
501        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
502      END DO
503      !                                                ! Sea surface elevation time stepping
504      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
505      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
506      z1_rau0 = 0.5 / rau0
507      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
508
509      IF( .NOT.ln_dynspg_ts ) THEN
510      ! These lines are not necessary with time splitting since
511      ! boundary condition on sea level is set during ts loop
512#if defined key_agrif
513      CALL agrif_ssh( kt )
514#endif
515         IF( ln_bdy ) THEN
516            ssha(:,:) = ssha(:,:) * bdytmask(:,:)
517            CALL lbc_lnk( ssha, 'T', 1. ) 
518         ENDIF
519      ENDIF
520      !
521      !                                           !------------------------------!
522      !                                           !     Now Vertical Velocity    !
523      !                                           !------------------------------!
524      z1_2dt = 1.e0 / z2dt
525      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
526         ! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise
527         wn(:,:,jk) = wn(:,:,jk+1) -   e3t_n(:,:,jk) * hdivn(:,:,jk)        &
528            &                      - ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )    &
529            &                         * tmask(:,:,jk) * z1_2dt
530         IF( ln_bdy ) wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
531      END DO
532      !
533      CALL wrk_dealloc( jpi,jpj,   zhdiv )
534      !
535      IF( nn_timing == 1 )  CALL timing_stop('trc_sub_ssh')
536      !
537   END SUBROUTINE trc_sub_ssh
538
539
540   INTEGER FUNCTION trc_sub_alloc()
541      !!-------------------------------------------------------------------
542      !!                    *** ROUTINE trc_sub_alloc ***
543      !!-------------------------------------------------------------------
544      USE lib_mpp, ONLY: ctl_warn
545      INTEGER ::  ierr(3)
546      !!-------------------------------------------------------------------
547      !
548      ierr(:) = 0
549      !
550      ALLOCATE( un_temp(jpi,jpj,jpk)      ,  vn_temp(jpi,jpj,jpk)   ,     &
551         &      wn_temp(jpi,jpj,jpk)      ,                               &
552         &      rhop_temp(jpi,jpj,jpk)    ,  rhop_tm(jpi,jpj,jpk)   ,     &
553         &      sshn_temp(jpi,jpj)        ,  sshb_temp(jpi,jpj)     ,     &
554         &      ssha_temp(jpi,jpj)        ,                               &
555         &      rnf_temp(jpi,jpj)         ,  h_rnf_temp(jpi,jpj)    ,     &
556         &      tsn_temp(jpi,jpj,jpk,2)   ,  emp_b_temp(jpi,jpj)    ,     &
557         &      emp_temp(jpi,jpj)         ,  fmmflx_temp(jpi,jpj)   ,     &
558         &      hmld_temp(jpi,jpj)        ,  qsr_temp(jpi,jpj)      ,     &
559         &      fr_i_temp(jpi,jpj)        ,  fr_i_tm(jpi,jpj)       ,     &
560         &      wndm_temp(jpi,jpj)        ,  wndm_tm(jpi,jpj)       ,     &
561         &      avs_tm(jpi,jpj,jpk)       ,  avs_temp(jpi,jpj,jpk)  ,     &
562         &      hdivn_temp(jpi,jpj,jpk)   ,  hdivb_temp(jpi,jpj,jpk),     &
563         &      un_tm(jpi,jpj,jpk)        ,  vn_tm(jpi,jpj,jpk)     ,     &
564         &      sshn_tm(jpi,jpj)          ,  sshb_hold(jpi,jpj)     ,     &
565         &      tsn_tm(jpi,jpj,jpk,2)     ,                               &
566         &      emp_tm(jpi,jpj)           ,  fmmflx_tm(jpi,jpj)     ,     &
567         &      emp_b_hold(jpi,jpj)       ,                               &
568         &      hmld_tm(jpi,jpj)          ,  qsr_tm(jpi,jpj)        ,     &
569         &      rnf_tm(jpi,jpj)           ,  h_rnf_tm(jpi,jpj)      , STAT=ierr(1) ) 
570      !
571      IF( l_ldfslp ) THEN
572         ALLOCATE( uslp_temp(jpi,jpj,jpk) ,  wslpi_temp(jpi,jpj,jpk),     &
573            &      vslp_temp(jpi,jpj,jpk) ,  wslpj_temp(jpi,jpj,jpk),     &
574            &      uslp_tm  (jpi,jpj,jpk) ,  wslpi_tm  (jpi,jpj,jpk),     &
575            &      vslp_tm  (jpi,jpj,jpk) ,  wslpj_tm  (jpi,jpj,jpk), STAT=ierr(2) )
576      ENDIF
577      IF( ln_trabbl ) THEN
578         ALLOCATE( ahu_bbl_temp(jpi,jpj)  , utr_bbl_temp(jpi,jpj)   ,     &
579            &      ahv_bbl_temp(jpi,jpj)  , vtr_bbl_temp(jpi,jpj)   ,     &
580            &      ahu_bbl_tm  (jpi,jpj)  , utr_bbl_tm  (jpi,jpj)   ,     &
581            &      ahv_bbl_tm  (jpi,jpj)  , vtr_bbl_tm  (jpi,jpj)   , STAT=ierr(3) ) 
582      ENDIF
583      !
584      trc_sub_alloc = MAXVAL( ierr )
585      !
586      IF( trc_sub_alloc /= 0 )   CALL ctl_warn('trc_sub_alloc: failed to allocate arrays')
587      !
588   END FUNCTION trc_sub_alloc
589
590#else
591   !!----------------------------------------------------------------------
592   !!   Default key                                     NO passive tracers
593   !!----------------------------------------------------------------------
594CONTAINS
595   SUBROUTINE trc_sub_stp( kt )        ! Empty routine
596      WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
597   END SUBROUTINE trc_sub_stp
598   SUBROUTINE trc_sub_ini        ! Empty routine
599      WRITE(*,*) 'trc_sub_ini: You should not have seen this print! error?', kt
600   END SUBROUTINE trc_sub_ini
601#endif
602
603   !!======================================================================
604END MODULE trcsub
Note: See TracBrowser for help on using the repository browser.