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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/LDF/ldftra.F90 @ 10946

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

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