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 NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/LDF – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/LDF/ldftra.F90 @ 11671

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

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