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/trunk/src/OCE/LDF – NEMO

source: NEMO/trunk/src/OCE/LDF/ldftra.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

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