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.
ldftra.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 48.5 KB
Line 
1MODULE ldftra
2   !!======================================================================
3   !!                       ***  MODULE  ldftra  ***
4   !! Ocean physics:  lateral diffusivity coefficients
5   !!=====================================================================
6   !! History :       ! 1997-07  (G. Madec)  from inimix.F split in 2 routines
7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
8   !!            2.0  ! 2005-11  (G. Madec) 
9   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  restructuration/simplification of aht/aeiv specification,
10   !!                 !                                  add velocity dependent coefficient and optional read in file
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ldf_tra_init : initialization, namelist read, and parameters control
15   !!   ldf_tra      : update lateral eddy diffusivity coefficients at each time step
16   !!   ldf_eiv_init : initialization of the eiv coeff. from namelist choices
17   !!   ldf_eiv      : time evolution of the eiv coefficients (function of the growth rate of baroclinic instability)
18   !!   ldf_eiv_trp  : add to the input ocean transport the contribution of the EIV parametrization
19   !!   ldf_eiv_dia  : diagnose the eddy induced velocity from the eiv streamfunction
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers
22   USE dom_oce         ! ocean space and time domain
23   USE phycst          ! physical constants
24   USE ldfslp          ! lateral diffusion: slope of iso-neutral surfaces
25   USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases
26   USE diaptr
27   !
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O module for ehanced bottom friction file
30   USE lib_mpp         ! distribued memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! work arrays
33   USE timing          ! timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   ldf_tra_init   ! called by nemogcm.F90
39   PUBLIC   ldf_tra        ! called by step.F90
40   PUBLIC   ldf_eiv_init   ! called by nemogcm.F90
41   PUBLIC   ldf_eiv        ! called by step.F90
42   PUBLIC   ldf_eiv_trp    ! called by traadv.F90
43   PUBLIC   ldf_eiv_dia    ! called by traldf_iso and traldf_iso_triad.F90
44   
45   !                                   !!* Namelist namtra_ldf : lateral mixing on tracers *
46   !                                    != Operator type =!
47   LOGICAL , PUBLIC ::   ln_traldf_lap       !: laplacian operator
48   LOGICAL , PUBLIC ::   ln_traldf_blp       !: bilaplacian operator
49   !                                    != Direction of action =!
50   LOGICAL , PUBLIC ::   ln_traldf_lev       !: iso-level direction
51   LOGICAL , PUBLIC ::   ln_traldf_hor       !: horizontal (geopotential) direction
52!  LOGICAL , PUBLIC ::   ln_traldf_iso       !: iso-neutral direction                    (see ldfslp)
53!  LOGICAL , PUBLIC ::   ln_traldf_triad     !: griffies triad scheme                    (see ldfslp)
54   LOGICAL , PUBLIC ::   ln_traldf_msc       !: Method of Stabilizing Correction
55!  LOGICAL , PUBLIC ::   ln_triad_iso        !: pure horizontal mixing in ML             (see ldfslp)
56!  LOGICAL , PUBLIC ::   ln_botmix_triad     !: mixing on bottom                         (see ldfslp)
57!  REAL(wp), PUBLIC ::   rn_sw_triad         !: =1/0 switching triad / all 4 triads used (see ldfslp)
58!  REAL(wp), PUBLIC ::   rn_slpmax           !: slope limit                              (see ldfslp)
59   !                                    !=  Coefficients =!
60   INTEGER , PUBLIC ::   nn_aht_ijk_t        !: choice of time & space variations of the lateral eddy diffusivity coef.
61   REAL(wp), PUBLIC ::   rn_aht_0            !:   laplacian lateral eddy diffusivity [m2/s]
62   REAL(wp), PUBLIC ::   rn_bht_0            !: bilaplacian lateral eddy diffusivity [m4/s]
63
64   !                                   !!* Namelist namtra_ldfeiv : eddy induced velocity param. *
65   !                                    != Use/diagnose eiv =!
66   LOGICAL , PUBLIC ::   ln_ldfeiv           !: eddy induced velocity flag
67   LOGICAL , PUBLIC ::   ln_ldfeiv_dia       !: diagnose & output eiv streamfunction and velocity (IOM)
68   !                                    != Coefficients =!
69   INTEGER , PUBLIC ::   nn_aei_ijk_t        !: choice of time/space variation of the eiv coeff.
70   REAL(wp), PUBLIC ::   rn_aeiv_0           !: eddy induced velocity coefficient [m2/s]
71   
72   LOGICAL , PUBLIC ::   l_ldftra_time = .FALSE.   !: flag for time variation of the lateral eddy diffusivity coef.
73   LOGICAL , PUBLIC ::   l_ldfeiv_time = .FALSE.   ! flag for time variation of the eiv coef.
74
75   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahtu, ahtv   !: eddy diffusivity coef. at U- and V-points   [m2/s]
76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aeiu, aeiv   !: eddy induced velocity coeff.                [m2/s]
77
78   REAL(wp) ::   r1_4  = 0.25_wp          ! =1/4
79   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! =1/12
80
81   !! * Substitutions
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
85   !! $Id$
86   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
87   !!----------------------------------------------------------------------
88CONTAINS
89
90   SUBROUTINE ldf_tra_init
91      !!----------------------------------------------------------------------
92      !!                  ***  ROUTINE ldf_tra_init  ***
93      !!
94      !! ** Purpose :   initializations of the tracer lateral mixing coeff.
95      !!
96      !! ** Method  : * the eddy diffusivity coef. specification depends on:
97      !!
98      !!    ln_traldf_lap = T     laplacian operator
99      !!    ln_traldf_blp = T   bilaplacian operator
100      !!
101      !!    nn_aht_ijk_t  =  0 => = constant
102      !!                  !
103      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
104      !!                  !
105      !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file
106      !!                  = 20    = F(i,j)   = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
107      !!                  = 21    = F(i,j,t) = F(growth rate of baroclinic instability)
108      !!                  !
109      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity.nc' file
110      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
111      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
112      !!                                                          or |u|e^3/12 bilaplacian operator )
113      !!              * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init
114      !!           
115      !! ** action  : ahtu, ahtv initialized once for all or l_ldftra_time set to true
116      !!              aeiu, aeiv initialized once for all or l_ldfeiv_time set to true
117      !!----------------------------------------------------------------------
118      INTEGER  ::   jk, jj, ji        ! dummy loop indices
119      INTEGER  ::   ierr, inum, ios   ! local integer
120      REAL(wp) ::   zah0              ! local scalar
121      !
122      NAMELIST/namtra_ldf/ ln_traldf_lap, ln_traldf_blp  ,                   &   ! type of operator
123         &                 ln_traldf_lev, ln_traldf_hor  , ln_traldf_triad,  &   ! acting direction of the operator
124         &                 ln_traldf_iso, ln_traldf_msc  ,  rn_slpmax     ,  &   ! option for iso-neutral operator
125         &                 ln_triad_iso , ln_botmix_triad, rn_sw_triad    ,  &   ! option for triad operator
126         &                 rn_aht_0     , rn_bht_0       , nn_aht_ijk_t          ! lateral eddy coefficient
127      !!----------------------------------------------------------------------
128      !
129      !  Choice of lateral tracer physics
130      ! =================================
131      !
132      REWIND( numnam_ref )              ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers
133      READ  ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901)
134901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist', lwp )
135      !
136      REWIND( numnam_cfg )              ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers
137      READ  ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 )
138902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist', lwp )
139      IF(lwm) WRITE ( numond, namtra_ldf )
140      !
141      IF(lwp) THEN                      ! control print
142         WRITE(numout,*)
143         WRITE(numout,*) 'ldf_tra_init : lateral tracer physics'
144         WRITE(numout,*) '~~~~~~~~~~~~ '
145         WRITE(numout,*) '   Namelist namtra_ldf : lateral mixing parameters (type, direction, coefficients)'
146         !
147         WRITE(numout,*) '      type :'
148         WRITE(numout,*) '         laplacian operator                      ln_traldf_lap   = ', ln_traldf_lap
149         WRITE(numout,*) '         bilaplacian operator                    ln_traldf_blp   = ', ln_traldf_blp
150         !
151         WRITE(numout,*) '      direction of action :'
152         WRITE(numout,*) '         iso-level                               ln_traldf_lev   = ', ln_traldf_lev
153         WRITE(numout,*) '         horizontal (geopotential)               ln_traldf_hor   = ', ln_traldf_hor
154         WRITE(numout,*) '         iso-neutral Madec operator              ln_traldf_iso   = ', ln_traldf_iso
155         WRITE(numout,*) '         iso-neutral triad operator              ln_traldf_triad = ', ln_traldf_triad
156         WRITE(numout,*) '            iso-neutral (Method of Stab. Corr.)  ln_traldf_msc   = ', ln_traldf_msc
157         WRITE(numout,*) '            maximum isoppycnal slope             rn_slpmax       = ', rn_slpmax
158         WRITE(numout,*) '            pure lateral mixing in ML            ln_triad_iso    = ', ln_triad_iso
159         WRITE(numout,*) '            switching triad or not               rn_sw_triad     = ', rn_sw_triad
160         WRITE(numout,*) '            lateral mixing on bottom             ln_botmix_triad = ', ln_botmix_triad
161         !
162         WRITE(numout,*) '      coefficients :'
163         WRITE(numout,*) '         lateral eddy diffusivity   (lap case)   rn_aht_0        = ', rn_aht_0
164         WRITE(numout,*) '         lateral eddy diffusivity (bilap case)   rn_bht_0        = ', rn_bht_0
165         WRITE(numout,*) '         type of time-space variation            nn_aht_ijk_t    = ', nn_aht_ijk_t
166      ENDIF
167      !
168      !                                ! Parameter control
169      !
170      IF( .NOT.ln_traldf_lap .AND. .NOT.ln_traldf_blp ) THEN
171         IF(lwp) WRITE(numout,*) '   No diffusive operator selected. ahtu and ahtv are not allocated'
172         l_ldftra_time = .FALSE.
173         RETURN
174      ENDIF
175      !
176      IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN     ! iso-neutral bilaplacian need MSC
177         IF( .NOT.ln_traldf_msc )   CALL ctl_stop( 'tra_ldf_init: iso-neutral bilaplacian requires ln_traldf_msc=.true.' )
178      ENDIF
179      !
180      !  Space/time variation of eddy coefficients
181      ! ===========================================
182      !                                               ! allocate the aht arrays
183      ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr )
184      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays')
185      !
186!$OMP PARALLEL DO schedule(static) private(jj, ji)
187      DO jj = 1, jpj
188         DO ji = 1, jpi
189            ahtu(ji,jj,jpk) = 0._wp                           ! last level always 0 
190            ahtv(ji,jj,jpk) = 0._wp
191         END DO
192      END DO
193      !
194      !                                               ! value of eddy mixing coef.
195      IF    ( ln_traldf_lap ) THEN   ;   zah0 =      rn_aht_0        !   laplacian operator
196      ELSEIF( ln_traldf_blp ) THEN   ;   zah0 = ABS( rn_bht_0 )      ! bilaplacian operator
197      ENDIF
198      !
199      l_ldftra_time = .FALSE.                         ! no time variation except in case defined below
200      !
201      IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN     ! only if a lateral diffusion operator is used
202         !
203         SELECT CASE(  nn_aht_ijk_t  )                   ! Specification of space time variations of ehtu, ahtv
204         !
205         CASE(   0  )      !==  constant  ==!
206            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = constant = ', rn_aht_0
207!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
208            DO jk = 1, jpk
209               DO jj = 1, jpj
210                  DO ji = 1, jpi
211                     ahtu(ji,jj,jk) = zah0 * umask(ji,jj,jk)
212                     ahtv(ji,jj,jk) = zah0 * vmask(ji,jj,jk)
213                  END DO
214               END DO
215            END DO
216            !
217         CASE(  10  )      !==  fixed profile  ==!
218            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( depth )'
219!$OMP PARALLEL DO schedule(static) private(jj, ji)
220            DO jj = 1, jpj
221               DO ji = 1, jpi
222                  ahtu(ji,jj,1) = zah0 * umask(ji,jj,1)                      ! constant surface value
223                  ahtv(ji,jj,1) = zah0 * vmask(ji,jj,1)
224               END DO
225            END DO
226            CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )
227            !
228         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
229            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file'
230            CALL iom_open( 'eddy_diffusivity_2D.nc', inum )
231            CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) )
232            CALL iom_get ( inum, jpdom_data, 'ahtv_2D', ahtv(:,:,1) )
233            CALL iom_close( inum )
234!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
235            DO jk = 2, jpkm1
236               DO jj = 1, jpj
237                  DO ji = 1, jpi
238                     ahtu(ji,jj,jk) = ahtu(ji,jj,1) * umask(ji,jj,jk)
239                     ahtv(ji,jj,jk) = ahtv(ji,jj,1) * vmask(ji,jj,jk)
240                  END DO
241               END DO
242            END DO
243            !
244         CASE(  20  )      !== fixed horizontal shape  ==!
245            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)'
246            IF( ln_traldf_lap )   CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
247            IF( ln_traldf_blp )   CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
248            !
249         CASE(  21  )      !==  time varying 2D field  ==!
250            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( latitude, longitude, time )'
251            IF(lwp) WRITE(numout,*) '                              = F( growth rate of baroclinic instability )'
252            IF(lwp) WRITE(numout,*) '                              min value = 0.1 * rn_aht_0'
253            IF(lwp) WRITE(numout,*) '                              max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)'
254            IF(lwp) WRITE(numout,*) '                              increased to rn_aht_0 within 20N-20S'
255            !
256            l_ldftra_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
257            !
258            IF( ln_traldf_blp ) THEN
259               CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' )
260            ENDIF
261            !
262         CASE( -30  )      !== fixed 3D shape read in file  ==!
263            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F(i,j,k) read in eddy_diffusivity.nc file'
264            CALL iom_open( 'eddy_diffusivity_3D.nc', inum )
265            CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu )
266            CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv )
267            CALL iom_close( inum )
268!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
269            DO jk = 1, jpkm1
270               DO jj = 1, jpj
271                  DO ji = 1, jpi
272                     ahtu(ji,jj,jk) = ahtu(ji,jj,jk) * umask(ji,jj,jk)
273                     ahtv(ji,jj,jk) = ahtv(ji,jj,jk) * vmask(ji,jj,jk)
274                  END DO
275               END DO
276            END DO
277            !
278         CASE(  30  )      !==  fixed 3D shape  ==!
279            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( latitude, longitude, depth )'
280            IF( ln_traldf_lap )   CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
281            IF( ln_traldf_blp )   CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv )    ! surface value proportional to scale factor
282            !                                                    ! reduction with depth
283            CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )
284            !
285         CASE(  31  )      !==  time varying 3D field  ==!
286            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( latitude, longitude, depth , time )'
287            IF(lwp) WRITE(numout,*) '                                proportional to the velocity : |u|e/12 or |u|e^3/12'
288            !
289            l_ldftra_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
290            !
291         CASE DEFAULT
292            CALL ctl_stop('ldf_tra_init: wrong choice for nn_aht_ijk_t, the type of space-time variation of aht')
293         END SELECT
294         !
295         IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN
296!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
297            DO jk = 1, jpk
298               DO jj = 1, jpj
299                  DO ji = 1, jpi
300                     ahtu(ji,jj,jk) = SQRT( ahtu(ji,jj,jk) )
301                     ahtv(ji,jj,jk) = SQRT( ahtv(ji,jj,jk) )
302                  END DO
303               END DO
304            END DO
305         ENDIF
306         !
307      ENDIF
308      !
309   END SUBROUTINE ldf_tra_init
310
311
312   SUBROUTINE ldf_tra( kt )
313      !!----------------------------------------------------------------------
314      !!                  ***  ROUTINE ldf_tra  ***
315      !!
316      !! ** Purpose :   update at kt the tracer lateral mixing coeff. (aht and aeiv)
317      !!
318      !! ** Method  :   time varying eddy diffusivity coefficients:
319      !!
320      !!    nn_aei_ijk_t = 21    aeiu, aeiv = F(i,j,  t) = F(growth rate of baroclinic instability)
321      !!                                                   with a reduction to 0 in vicinity of the Equator
322      !!    nn_aht_ijk_t = 21    ahtu, ahtv = F(i,j,  t) = F(growth rate of baroclinic instability)
323      !!
324      !!                 = 31    ahtu, ahtv = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
325      !!                                                                     or |u|e^3/12 bilaplacian operator )
326      !!
327      !! ** action  :   ahtu, ahtv   update at each time step   
328      !!                aeiu, aeiv      -       -     -    -   (if ln_ldfeiv=T)
329      !!----------------------------------------------------------------------
330      INTEGER, INTENT(in) ::   kt   ! time step
331      !
332      INTEGER  ::   ji, jj, jk   ! dummy loop indices
333      REAL(wp) ::   zaht, zahf, zaht_min, z1_f20       ! local scalar
334      !!----------------------------------------------------------------------
335      !
336      IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN       ! eddy induced velocity coefficients
337         !                                ! =F(growth rate of baroclinic instability)
338         !                                ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S
339         CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv )
340      ENDIF
341      !
342      SELECT CASE(  nn_aht_ijk_t  )       ! Eddy diffusivity coefficients
343      !
344      CASE(  21  )       !==  time varying 2D field  ==!   = F( growth rate of baroclinic instability )
345         !                                             !   min value rn_aht_0 / 10
346         !                                             !   max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)
347         !                                             !   increase to rn_aht_0 within 20N-20S
348         IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN   ! use the already computed aei.
349!$OMP PARALLEL DO schedule(static) private(jj,ji)
350            DO jj = 1, jpj
351               DO ji = 1, jpi
352                  ahtu(ji,jj,1) = aeiu(ji,jj,1)
353                  ahtv(ji,jj,1) = aeiv(ji,jj,1)
354               END DO
355            END DO
356         ELSE                                            ! compute aht.
357            CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv )
358         ENDIF
359         !
360         z1_f20   = 1._wp / (  2._wp * omega * SIN( rad * 20._wp )  )      ! 1 / ff(20 degrees)   
361         zaht_min = 0.2_wp * rn_aht_0                                      ! minimum value for aht
362!$OMP PARALLEL
363!$OMP DO schedule(static) private(jj,ji,zaht,zahf)
364         DO jj = 1, jpj
365            DO ji = 1, jpi
366               !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg)
367               !!     ==>>>   The Coriolis value is identical for t- & u_points, and for v- and f-points
368               zaht = ( 1._wp -  MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )
369               zahf = ( 1._wp -  MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )
370               ahtu(ji,jj,1) = (  MAX( zaht_min, ahtu(ji,jj,1) ) + zaht  ) * umask(ji,jj,1)     ! min value zaht_min
371               ahtv(ji,jj,1) = (  MAX( zaht_min, ahtv(ji,jj,1) ) + zahf  ) * vmask(ji,jj,1)     ! increase within 20S-20N
372            END DO
373         END DO
374!$OMP DO schedule(static) private(jk,jj,ji)
375         DO jk = 2, jpkm1                             ! deeper value = surface value
376            DO jj = 1, jpj
377               DO ji = 1, jpi
378                  ahtu(ji,jj,jk) = ahtu(ji,jj,1) * umask(ji,jj,jk)
379                  ahtv(ji,jj,jk) = ahtv(ji,jj,1) * vmask(ji,jj,jk)
380               END DO
381            END DO
382         END DO
383!$OMP END PARALLEL
384         !
385      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
386         IF( ln_traldf_lap     ) THEN          !   laplacian operator |u| e /12
387!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
388            DO jk = 1, jpkm1
389               DO jj = 1, jpj
390                  DO ji = 1, jpi
391                     ahtu(ji,jj,jk) = ABS( ub(ji,jj,jk) ) * e1u(ji,jj) * r1_12
392                     ahtv(ji,jj,jk) = ABS( vb(ji,jj,jk) ) * e2v(ji,jj) * r1_12
393                  END DO
394               END DO
395            END DO
396         ELSEIF( ln_traldf_blp ) THEN      ! bilaplacian operator      sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e
397!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
398            DO jk = 1, jpkm1
399               DO jj = 1, jpj
400                  DO ji = 1, jpi
401                     ahtu(ji,jj,jk) = SQRT(  ABS( ub(ji,jj,jk) ) * e1u(ji,jj) * r1_12  ) * e1u(ji,jj)
402                     ahtv(ji,jj,jk) = SQRT(  ABS( vb(ji,jj,jk) ) * e2v(ji,jj) * r1_12  ) * e2v(ji,jj)
403                  END DO
404               END DO
405            END DO
406         ENDIF
407         !
408      END SELECT
409      !
410      CALL iom_put( "ahtu_2d", ahtu(:,:,1) )   ! surface u-eddy diffusivity coeff.
411      CALL iom_put( "ahtv_2d", ahtv(:,:,1) )   ! surface v-eddy diffusivity coeff.
412      CALL iom_put( "ahtu_3d", ahtu(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
413      CALL iom_put( "ahtv_3d", ahtv(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
414      !
415!!gm  : THE IF below is to be checked (comes from Seb)
416      IF( ln_ldfeiv ) THEN
417        CALL iom_put( "aeiu_2d", aeiu(:,:,1) )   ! surface u-EIV coeff.
418        CALL iom_put( "aeiv_2d", aeiv(:,:,1) )   ! surface v-EIV coeff.
419        CALL iom_put( "aeiu_3d", aeiu(:,:,:) )   ! 3D      u-EIV coeff.
420        CALL iom_put( "aeiv_3d", aeiv(:,:,:) )   ! 3D      v-EIV coeff.
421      ENDIF
422      !
423   END SUBROUTINE ldf_tra
424
425
426   SUBROUTINE ldf_eiv_init
427      !!----------------------------------------------------------------------
428      !!                  ***  ROUTINE ldf_eiv_init  ***
429      !!
430      !! ** Purpose :   initialization of the eiv coeff. from namelist choices.
431      !!
432      !! ** Method :
433      !!
434      !! ** Action :   aeiu , aeiv   : EIV coeff. at u- & v-points
435      !!               l_ldfeiv_time : =T if EIV coefficients vary with time
436      !!----------------------------------------------------------------------
437      INTEGER  ::   jk, jj, ji        ! dummy loop indices
438      INTEGER  ::   ierr, inum, ios   ! local integer
439      !
440      NAMELIST/namtra_ldfeiv/ ln_ldfeiv   , ln_ldfeiv_dia,   &    ! eddy induced velocity (eiv)
441         &                    nn_aei_ijk_t, rn_aeiv_0             ! eiv  coefficient
442      !!----------------------------------------------------------------------
443      !
444      REWIND( numnam_ref )              ! Namelist namtra_ldfeiv in reference namelist : eddy induced velocity param.
445      READ  ( numnam_ref, namtra_ldfeiv, IOSTAT = ios, ERR = 901)
446901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldfeiv in reference namelist', lwp )
447      !
448      REWIND( numnam_cfg )              ! Namelist namtra_ldfeiv in configuration namelist : eddy induced velocity param.
449      READ  ( numnam_cfg, namtra_ldfeiv, IOSTAT = ios, ERR = 902 )
450902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldfeiv in configuration namelist', lwp )
451      IF(lwm)  WRITE ( numond, namtra_ldfeiv )
452
453      IF(lwp) THEN                      ! control print
454         WRITE(numout,*)
455         WRITE(numout,*) 'ldf_eiv_init : eddy induced velocity parametrization'
456         WRITE(numout,*) '~~~~~~~~~~~~ '
457         WRITE(numout,*) '   Namelist namtra_ldfeiv : '
458         WRITE(numout,*) '      Eddy Induced Velocity (eiv) param.      ln_ldfeiv     = ', ln_ldfeiv
459         WRITE(numout,*) '      eiv streamfunction & velocity diag.     ln_ldfeiv_dia = ', ln_ldfeiv_dia
460         WRITE(numout,*) '      eddy induced velocity coef.             rn_aeiv_0     = ', rn_aeiv_0
461         WRITE(numout,*) '      type of time-space variation            nn_aei_ijk_t  = ', nn_aei_ijk_t
462         WRITE(numout,*)
463      ENDIF
464      !
465      IF( ln_ldfeiv .AND. ln_traldf_blp )   CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' )
466
467      !                                 ! Parameter control
468      l_ldfeiv_time = .FALSE.   
469      !
470      IF( ln_ldfeiv ) THEN                         ! allocate the aei arrays
471         ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr )
472         IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays')
473         !
474         SELECT CASE( nn_aei_ijk_t )               ! Specification of space time variations of eaiu, aeiv
475         !
476         CASE(   0  )      !==  constant  ==!
477            IF(lwp) WRITE(numout,*) '          eddy induced velocity coef. = constant = ', rn_aeiv_0
478!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
479            DO jk = 1, jpk
480               DO jj = 1, jpj
481                  DO ji = 1, jpi
482                     aeiu(ji,jj,jk) = rn_aeiv_0
483                     aeiv(ji,jj,jk) = rn_aeiv_0
484                  END DO
485               END DO
486            END DO
487            !
488         CASE(  10  )      !==  fixed profile  ==!
489            IF(lwp) WRITE(numout,*) '          eddy induced velocity coef. = F( depth )'
490!$OMP PARALLEL DO schedule(static) private(jj, ji)
491            DO jj = 1, jpj
492               DO ji = 1, jpi
493                  aeiu(ji,jj,1) = rn_aeiv_0                                ! constant surface value
494                  aeiv(ji,jj,1) = rn_aeiv_0
495               END DO
496            END DO
497            CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )
498            !
499         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
500            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file'
501            CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum )
502            CALL iom_get  ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) )
503            CALL iom_get  ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) )
504            CALL iom_close( inum )
505!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
506            DO jk = 2, jpk
507               DO jj = 1, jpj
508                  DO ji = 1, jpi
509                     aeiu(ji,jj,jk) = aeiu(ji,jj,1)
510                     aeiv(ji,jj,jk) = aeiv(ji,jj,1)
511                  END DO
512               END DO
513            END DO
514            !
515         CASE(  20  )      !== fixed horizontal shape  ==!
516            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)'
517            CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv )    ! surface value proportional to scale factor
518            !
519         CASE(  21  )       !==  time varying 2D field  ==!
520            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( latitude, longitude, time )'
521            IF(lwp) WRITE(numout,*) '                              = F( growth rate of baroclinic instability )'
522            !
523            l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90
524            !
525         CASE( -30  )      !== fixed 3D shape read in file  ==!
526            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
527            CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum )
528            CALL iom_get  ( inum, jpdom_data, 'aeiu', aeiu )
529            CALL iom_get  ( inum, jpdom_data, 'aeiv', aeiv )
530            CALL iom_close( inum )
531            !
532         CASE(  30  )       !==  fixed 3D shape  ==!
533            IF(lwp) WRITE(numout,*) '          tracer mixing coef. = F( latitude, longitude, depth )'
534            CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv )    ! surface value proportional to scale factor
535            !                                                 ! reduction with depth
536            CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv )
537            !
538         CASE DEFAULT
539            CALL ctl_stop('ldf_tra_init: wrong choice for nn_aei_ijk_t, the type of space-time variation of aei')
540         END SELECT
541         !
542      ELSE
543          IF(lwp) WRITE(numout,*) '   eddy induced velocity param is NOT used neither diagnosed'
544          ln_ldfeiv_dia = .FALSE.
545      ENDIF
546      !                   
547   END SUBROUTINE ldf_eiv_init
548
549
550   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )
551      !!----------------------------------------------------------------------
552      !!                  ***  ROUTINE ldf_eiv  ***
553      !!
554      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
555      !!              growth rate of baroclinic instability.
556      !!
557      !! ** Method  :   coefficient function of the growth rate of baroclinic instability
558      !!
559      !! Reference : Treguier et al. JPO 1997   ; Held and Larichev JAS 1996
560      !!----------------------------------------------------------------------
561      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
562      REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s]
563      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s]
564      !
565      INTEGER  ::   ji, jj, jk    ! dummy loop indices
566      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei   ! local scalars
567      REAL(wp), DIMENSION(:,:), POINTER ::   zn, zah, zhw, zross, zaeiw   ! 2D workspace
568      !!----------------------------------------------------------------------
569      !
570      IF( nn_timing == 1 )   CALL timing_start('ldf_eiv')
571      !
572      CALL wrk_alloc( jpi,jpj,   zn, zah, zhw, zross, zaeiw )
573      !     
574!$OMP PARALLEL DO schedule(static) private(jj,ji)
575      DO jj = 1, jpj
576         DO ji = 1, jpi
577            zn   (ji,jj) = 0._wp      ! Local initialization
578            zhw  (ji,jj) = 5._wp
579            zah  (ji,jj) = 0._wp
580            zross(ji,jj) = 0._wp
581         END DO
582      END DO
583      !                       ! Compute lateral diffusive coefficient at T-point
584      IF( ln_traldf_triad ) THEN
585         DO jk = 1, jpk
586!$OMP PARALLEL DO schedule(static) private(jj,ji,zn2,ze3w)
587            DO jj = 2, jpjm1
588               DO ji = 2, jpim1
589                  ! Take the max of N^2 and zero then take the vertical sum
590                  ! of the square root of the resulting N^2 ( required to compute
591                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
592                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
593                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
594                  ! Compute elements required for the inverse time scale of baroclinic
595                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
596                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
597                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
598                  zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
599                  zhw(ji,jj) = zhw(ji,jj) + ze3w
600               END DO
601            END DO
602         END DO
603      ELSE
604         DO jk = 1, jpk
605!$OMP PARALLEL DO schedule(static) private(jj,ji,zn2,ze3w)
606            DO jj = 2, jpjm1
607               DO ji = 2, jpim1
608                  ! Take the max of N^2 and zero then take the vertical sum
609                  ! of the square root of the resulting N^2 ( required to compute
610                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
611                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
612                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
613                  ! Compute elements required for the inverse time scale of baroclinic
614                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
615                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
616                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
617                  zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
618                     &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
619                  zhw(ji,jj) = zhw(ji,jj) + ze3w
620               END DO
621            END DO
622         END DO
623      END IF
624
625!$OMP PARALLEL
626!$OMP DO schedule(static) private(jj,ji,zfw)
627      DO jj = 2, jpjm1
628         DO ji = fs_2, fs_jpim1   ! vector opt.
629            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
630            ! Rossby radius at w-point taken < 40km and  > 2km
631            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 )
632            ! Compute aeiw by multiplying Ro^2 and T^-1
633            zaeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
634         END DO
635      END DO
636
637      !                                         !==  Bound on eiv coeff.  ==!
638      z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  )
639!$OMP DO schedule(static) private(jj,ji,zzaei)
640      DO jj = 2, jpjm1
641         DO ji = fs_2, fs_jpim1   ! vector opt.
642            zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)       ! tropical decrease
643            zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0
644         END DO
645      END DO
646!$OMP END PARALLEL
647      CALL lbc_lnk( zaeiw(:,:), 'W', 1. )       ! lateral boundary condition
648      !               
649!$OMP PARALLEL DO schedule(static) private(jj,ji)
650      DO jj = 2, jpjm1                          !== aei at u- and v-points  ==!
651         DO ji = fs_2, fs_jpim1   ! vector opt.
652            paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1)
653            paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji  ,jj+1) ) * vmask(ji,jj,1)
654         END DO
655      END DO
656      CALL lbc_lnk( paeiu(:,:,1), 'U', 1. )   ;   CALL lbc_lnk( paeiv(:,:,1), 'V', 1. )      ! lateral boundary condition
657
658!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
659      DO jk = 2, jpkm1                          !==  deeper values equal the surface one  ==!
660         DO jj = 1, jpj
661            DO ji = 1, jpi
662               paeiu(ji,jj,jk) = paeiu(ji,jj,1) * umask(ji,jj,jk)
663               paeiv(ji,jj,jk) = paeiv(ji,jj,1) * vmask(ji,jj,jk)
664            END DO
665         END DO
666      END DO
667      ! 
668      CALL wrk_dealloc( jpi,jpj,   zn, zah, zhw, zross, zaeiw )
669      !
670      IF( nn_timing == 1 )   CALL timing_stop('ldf_eiv')
671      !
672   END SUBROUTINE ldf_eiv
673
674
675   SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype )
676      !!----------------------------------------------------------------------
677      !!                  ***  ROUTINE ldf_eiv_trp  ***
678      !!
679      !! ** Purpose :   add to the input ocean transport the contribution of
680      !!              the eddy induced velocity parametrization.
681      !!
682      !! ** Method  :   The eddy induced transport is computed from a flux stream-
683      !!              function which depends on the slope of iso-neutral surfaces
684      !!              (see ldf_slp). For example, in the i-k plan :
685      !!                   psi_uw = mk(aeiu) e2u mi(wslpi)   [in m3/s]
686      !!                   Utr_eiv = - dk[psi_uw]
687      !!                   Vtr_eiv = + di[psi_uw]
688      !!                ln_ldfeiv_dia = T : output the associated streamfunction,
689      !!                                    velocity and heat transport (call ldf_eiv_dia)
690      !!
691      !! ** Action  : pun, pvn increased by the eiv transport
692      !!----------------------------------------------------------------------
693      INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
694      INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index
695      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
696      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean transport components   [m3/s]
697      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean transport components   [m3/s]
698      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv                [m3/s]
699      !!
700      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
701      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
702      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
703      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpsi_uw, zpsi_vw
704      !!----------------------------------------------------------------------
705      !
706      IF( nn_timing == 1 )   CALL timing_start( 'ldf_eiv_trp')
707      !
708      CALL wrk_alloc( jpi,jpj,jpk,   zpsi_uw, zpsi_vw )
709
710      IF( kt == kit000 )  THEN
711         IF(lwp) WRITE(numout,*)
712         IF(lwp) WRITE(numout,*) 'ldf_eiv_trp : eddy induced advection on ', cdtype,' :'
713         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
714      ENDIF
715
716     
717!$OMP PARALLEL
718!$OMP DO schedule(static) private(jj,ji)
719      DO jj = 1, jpj
720         DO ji = 1, jpi
721            zpsi_uw(ji,jj, 1 ) = 0._wp   ;   zpsi_vw(ji,jj, 1 ) = 0._wp
722            zpsi_uw(ji,jj,jpk) = 0._wp   ;   zpsi_vw(ji,jj,jpk) = 0._wp
723         END DO
724      END DO
725!$OMP END DO NOWAIT
726      !
727!$OMP DO schedule(static) private(jk,jj,ji)
728      DO jk = 2, jpkm1
729         DO jj = 1, jpjm1
730            DO ji = 1, fs_jpim1   ! vector opt.
731               zpsi_uw(ji,jj,jk) = - 0.25_wp * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   &
732                  &                                       * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * umask(ji,jj,jk)
733               zpsi_vw(ji,jj,jk) = - 0.25_wp * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   &
734                  &                                       * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * vmask(ji,jj,jk)
735            END DO
736         END DO
737      END DO
738      !
739!$OMP DO schedule(static) private(jk,jj,ji)
740      DO jk = 1, jpkm1
741         DO jj = 1, jpjm1
742            DO ji = 1, fs_jpim1   ! vector opt.               
743               pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
744               pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
745            END DO
746         END DO
747      END DO
748!$OMP END DO NOWAIT
749!$OMP DO schedule(static) private(jk,jj,ji)
750      DO jk = 1, jpkm1
751         DO jj = 2, jpjm1
752            DO ji = fs_2, fs_jpim1   ! vector opt.
753               pwn(ji,jj,jk) = pwn(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   &
754                  &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) )
755            END DO
756         END DO
757      END DO
758!$OMP END PARALLEL
759      !
760      !                              ! diagnose the eddy induced velocity and associated heat transport
761      IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
762      !
763      CALL wrk_dealloc( jpi,jpj,jpk,   zpsi_uw, zpsi_vw )
764      !
765      IF( nn_timing == 1 )   CALL timing_stop( 'ldf_eiv_trp')
766      !
767    END SUBROUTINE ldf_eiv_trp
768
769
770   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )
771      !!----------------------------------------------------------------------
772      !!                  ***  ROUTINE ldf_eiv_dia  ***
773      !!
774      !! ** Purpose :   diagnose the eddy induced velocity and its associated
775      !!              vertically integrated heat transport.
776      !!
777      !! ** Method :
778      !!
779      !!----------------------------------------------------------------------
780      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   psi_uw, psi_vw   ! streamfunction   [m3/s]
781      !
782      INTEGER  ::   ji, jj, jk    ! dummy loop indices
783      REAL(wp) ::   zztmp   ! local scalar
784      REAL(wp), DIMENSION(:,:)  , POINTER ::   zw2d   ! 2D workspace
785      REAL(wp), DIMENSION(:,:,:), POINTER ::   zw3d   ! 3D workspace
786      !!----------------------------------------------------------------------
787      !
788      IF( nn_timing == 1 )  CALL timing_start( 'ldf_eiv_dia')
789      !
790      !                                                  !==  eiv stream function: output  ==!
791      CALL lbc_lnk( psi_uw, 'U', -1. )                         ! lateral boundary condition
792      CALL lbc_lnk( psi_vw, 'V', -1. )
793      !
794!!gm      CALL iom_put( "psi_eiv_uw", psi_uw )                 ! output
795!!gm      CALL iom_put( "psi_eiv_vw", psi_vw )
796      !
797      !                                                  !==  eiv velocities: calculate and output  ==!
798      CALL wrk_alloc( jpi,jpj,jpk,   zw3d )
799      !
800!$OMP PARALLEL
801!$OMP DO schedule(static) private(jj,ji)
802      DO jj = 1, jpj
803         DO ji = 1, jpi
804            zw3d(ji,jj,jpk) = 0._wp                            ! bottom value always 0
805         END DO
806      END DO
807!$OMP END DO NOWAIT
808      !
809!$OMP DO schedule(static) private(jk,jj,ji)
810      DO jk = 1, jpkm1                                         ! e2u e3u u_eiv = -dk[psi_uw]
811         DO jj = 1, jpj
812            DO ji = 1, jpi
813               zw3d(ji,jj,jk) = ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )
814            END DO
815         END DO
816      END DO
817!$OMP END PARALLEL
818      CALL iom_put( "uoce_eiv", zw3d )
819      !
820!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
821      DO jk = 1, jpkm1                                         ! e1v e3v v_eiv = -dk[psi_vw]
822         DO jj = 1, jpj
823            DO ji = 1, jpi
824               zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )
825            END DO
826         END DO
827      END DO
828      CALL iom_put( "voce_eiv", zw3d )
829      !
830!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
831      DO jk = 1, jpkm1                                         ! e1 e2 w_eiv = dk[psix] + dk[psix]
832         DO jj = 2, jpjm1
833            DO ji = fs_2, fs_jpim1  ! vector opt.
834               zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    &
835                  &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj)
836            END DO
837         END DO
838      END DO
839      CALL lbc_lnk( zw3d, 'T', 1. )      ! lateral boundary condition
840      CALL iom_put( "woce_eiv", zw3d )
841      !
842      !     
843      !
844      CALL wrk_alloc( jpi,jpj,   zw2d )
845      !
846      zztmp = 0.5_wp * rau0 * rcp 
847      IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
848!$OMP PARALLEL
849!$OMP DO schedule(static) private(jj,ji)
850         DO jj = 1, jpj
851            DO ji = 1, jpi
852               zw2d(ji,jj) = 0._wp
853            END DO
854         END DO
855!$OMP DO schedule(static) private(jk,jj,ji)
856         DO jk = 1, jpk
857            DO jj = 1, jpj
858               DO ji = 1, jpi
859                  zw3d(ji,jj,jk) = 0._wp 
860               END DO
861            END DO
862         END DO
863         DO jk = 1, jpkm1
864!$OMP DO schedule(static) private(jj,ji)
865            DO jj = 2, jpjm1
866               DO ji = fs_2, fs_jpim1   ! vector opt.
867                  zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
868                     &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) ) 
869                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
870               END DO
871            END DO
872         END DO
873!$OMP END PARALLEL
874         CALL lbc_lnk( zw2d, 'U', -1. )
875         CALL lbc_lnk( zw3d, 'U', -1. )
876         CALL iom_put( "ueiv_heattr"  , zztmp * zw2d )                  ! heat transport in i-direction
877         CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
878      ENDIF
879!$OMP PARALLEL
880!$OMP DO schedule(static) private(jj,ji)
881      DO jj = 1, jpj
882         DO ji = 1, jpi
883            zw2d(ji,jj) = 0._wp
884         END DO
885      END DO
886!$OMP DO schedule(static) private(jk,jj,ji)
887      DO jk = 1, jpk
888         DO jj = 1, jpj
889            DO ji = 1, jpi
890               zw3d(ji,jj,jk) = 0._wp
891            END DO
892         END DO
893      END DO
894      DO jk = 1, jpkm1
895!$OMP DO schedule(static) private(jj,ji)
896         DO jj = 2, jpjm1
897            DO ji = fs_2, fs_jpim1   ! vector opt.
898               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
899                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) ) 
900               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
901            END DO
902         END DO
903      END DO
904!$OMP END PARALLEL
905      CALL lbc_lnk( zw2d, 'V', -1. )
906      CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in j-direction
907      CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in j-direction
908      !
909      IF( ln_diaptr )  CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
910      !
911      zztmp = 0.5_wp * 0.5
912      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
913!$OMP PARALLEL
914!$OMP DO schedule(static) private(jj,ji)
915         DO jj = 1, jpj
916            DO ji = 1, jpi
917               zw2d(ji,jj) = 0._wp
918            END DO
919         END DO
920!$OMP DO schedule(static) private(jk,jj,ji)
921         DO jk = 1, jpk
922            DO jj = 1, jpj
923               DO ji = 1, jpi
924                  zw3d(ji,jj,jk) = 0._wp 
925               END DO
926            END DO
927         END DO
928         DO jk = 1, jpkm1
929!$OMP DO schedule(static) private(jj,ji)
930            DO jj = 2, jpjm1
931               DO ji = fs_2, fs_jpim1   ! vector opt.
932                  zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   &
933                     &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) ) 
934                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
935               END DO
936            END DO
937         END DO
938         CALL lbc_lnk( zw2d, 'U', -1. )
939         CALL lbc_lnk( zw3d, 'U', -1. )
940         CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
941         CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
942!$OMP END PARALLEL
943      ENDIF
944!$OMP PARALLEL
945!$OMP DO schedule(static) private(jj,ji)
946      DO jj = 1, jpj
947         DO ji = 1, jpi
948            zw2d(ji,jj) = 0._wp
949         END DO
950      END DO
951!$OMP DO schedule(static) private(jk,jj,ji)
952      DO jk = 1, jpk
953         DO jj = 1, jpj
954            DO ji = 1, jpi
955               zw3d(ji,jj,jk) = 0._wp
956            END DO
957         END DO
958      END DO
959      DO jk = 1, jpkm1
960!$OMP DO schedule(static) private(jj,ji)
961         DO jj = 2, jpjm1
962            DO ji = fs_2, fs_jpim1   ! vector opt.
963               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   &
964                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) ) 
965               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
966            END DO
967         END DO
968      END DO
969!$OMP END PARALLEL
970      CALL lbc_lnk( zw2d, 'V', -1. )
971      CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction
972      CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction
973      !
974      IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
975      !
976      CALL wrk_dealloc( jpi,jpj,   zw2d )
977      CALL wrk_dealloc( jpi,jpj,jpk,   zw3d )
978      !
979      IF( nn_timing == 1 )  CALL timing_stop( 'ldf_eiv_dia')     
980      !
981   END SUBROUTINE ldf_eiv_dia
982
983   !!======================================================================
984END MODULE ldftra
Note: See TracBrowser for help on using the repository browser.