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/UKMO/NEMO_4.0.3_GM_rossby_radius_ramp/src/OCE/LDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.3_GM_rossby_radius_ramp/src/OCE/LDF/ldftra.F90 @ 13883

Last change on this file since 13883 was 13883, checked in by davestorkey, 4 years ago

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