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.
tramle.F90 in NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/TRA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/TRA/tramle.F90 @ 12245

Last change on this file since 12245 was 12245, checked in by cguiavarch, 4 years ago

Merge changes from /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser:12178-12226
File size: 21.3 KB
Line 
1MODULE tramle
2   !!======================================================================
3   !!                    ***  MODULE  tramle  ***
4   !! Ocean tracers: Mixed Layer Eddy induced transport
5   !!======================================================================
6   !! History :  3.3  !  2010-08  (G. Madec)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_mle_trp   : update the effective transport with the Mixed Layer Eddy induced transport
11   !!   tra_mle_init  : initialisation of the Mixed Layer Eddy induced transport computation
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers variables
14   USE dom_oce        ! ocean space and time domain variables
15   USE phycst         ! physical constant
16   USE zdfmxl         ! mixed layer depth
17   !
18   USE in_out_manager ! I/O manager
19   USE iom            ! IOM library
20   USE lib_mpp        ! MPP library
21   USE lbclnk         ! lateral boundary condition / mpp link
22
23   USE zdfosm, ONLY  : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   tra_mle_trp        ! routine called in traadv.F90
29   PUBLIC   tra_mle_init   ! routine called in traadv.F90
30
31   !                                    !!* namelist namtra_mle *
32   LOGICAL, PUBLIC ::   ln_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
33   INTEGER         ::      nn_mle           ! MLE type: =0 standard Fox-Kemper ; =1 new formulation
34   INTEGER         ::      nn_mld_uv        ! space interpolation of MLD at u- & v-pts (0=min,1=averaged,2=max)
35   INTEGER         ::      nn_conv          ! =1 no MLE in case of convection ; =0 always MLE
36   REAL(wp)        ::      rn_ce            ! MLE coefficient
37   !                                        ! parameters used in nn_mle = 0 case
38   REAL(wp)        ::      rn_lf               ! typical scale of mixed layer front
39   REAL(wp)        ::      rn_time             ! time scale for mixing momentum across the mixed layer
40   !                                        ! parameters used in nn_mle = 1 case
41   REAL(wp)        ::      rn_lat              ! reference latitude for a 5 km scale of ML front
42   REAL(wp)        ::      rn_rho_c_mle        ! Density criterion for definition of MLD used by FK
43
44   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation
45   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
46   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case
47
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rfu, rfv   ! modified Coriolis parameter (f+tau) at u- & v-pts
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_ft      ! inverse of the modified Coriolis parameter at t-pts
50
51   !! * Substitutions
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
55   !! $Id$
56   !! Software governed by the CeCILL license (see ./LICENSE)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60  SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype )
61    !!----------------------------------------------------------------------
62    !!                  ***  ROUTINE tra_mle_trp  ***
63    !!
64    !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport
65    !!
66    !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced
67    !!              transport are computed as follows :
68    !!                zu_mle = dk[ zpsi_uw ]
69    !!                zv_mle = dk[ zpsi_vw ]
70    !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ]
71    !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc)
72    !!              and added to the input velocity :
73    !!                p.n = p.n + z._mle
74    !!
75    !! ** Action  : - (pun,pvn,pwn) increased by the mle transport
76    !!                CAUTION, the transport is not updated at the last line/raw
77    !!                         this may be a problem for some advection schemes
78    !!
79    !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
80    !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
81    !!----------------------------------------------------------------------
82    INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
83    INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index
84    CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
85    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
86    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
87    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
88    !
89    INTEGER  ::   ji, jj, jk          ! dummy loop indices
90    INTEGER  ::   ii, ij, ik, ikmax   ! local integers
91    REAL(wp) ::   zcuw, zmuw, zc      ! local scalar
92    REAL(wp) ::   zcvw, zmvw          !   -      -
93    INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
94    REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
95    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
96    !!----------------------------------------------------------------------
97    !
98    !
99    IF(ln_osm_mle) THEN
100       ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
101       !
102       !
103       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
104       CASE ( 0 )                                               != min of the 2 neighbour MLDs
105          DO jj = 1, jpjm1
106             DO ji = 1, fs_jpim1   ! vector opt.
107                zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) )
108                zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) )
109             END DO
110          END DO
111       CASE ( 1 )                                               != average of the 2 neighbour MLDs
112          DO jj = 1, jpjm1
113             DO ji = 1, fs_jpim1   ! vector opt.
114                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
115                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
116             END DO
117          END DO
118       CASE ( 2 )                                               != max of the 2 neighbour MLDs
119          DO jj = 1, jpjm1
120             DO ji = 1, fs_jpim1   ! vector opt.
121                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
122                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
123             END DO
124          END DO
125       END SELECT
126       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
127          DO jj = 1, jpjm1
128             DO ji = 1, fs_jpim1   ! vector opt.
129                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
130                     &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   &
131                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
132                !
133                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            &
134                     &           * dbdy_mle(ji,jj)  * MIN( 111.e3_wp , e2v(ji,jj) )   &
135                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
136             END DO
137          END DO
138          !
139       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
140          DO jj = 1, jpjm1
141             DO ji = 1, fs_jpim1   ! vector opt.
142                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               &
143                     &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )
144                !
145                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj)               &
146                     &                  * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) )
147             END DO
148          END DO
149       ENDIF
150
151    ELSE !do not use osn_mle
152       !                                      !==  MLD used for MLE  ==!
153       !                                                ! convert density into buoyancy
154       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
155       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
156          DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m)
157             DO jj = 1, jpj
158                DO ji = 1, jpi                          ! index of the w-level at the ML based
159                   IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
160                END DO
161             END DO
162          END DO
163       ENDIF
164
165       !
166       !
167       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
168       zbm (:,:) = 0._wp
169       zn2 (:,:) = 0._wp
170       DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
171          DO jj = 1, jpj
172             DO ji = 1, jpi
173                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
174                zmld(ji,jj) = zmld(ji,jj) + zc
175                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
176                zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
177             END DO
178          END DO
179       END DO
180
181       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
182       CASE ( 0 )                                               != min of the 2 neighbour MLDs
183          DO jj = 1, jpjm1
184             DO ji = 1, fs_jpim1   ! vector opt.
185                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
186                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
187             END DO
188          END DO
189       CASE ( 1 )                                               != average of the 2 neighbour MLDs
190          DO jj = 1, jpjm1
191             DO ji = 1, fs_jpim1   ! vector opt.
192                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
193                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
194             END DO
195          END DO
196       CASE ( 2 )                                               != max of the 2 neighbour MLDs
197          DO jj = 1, jpjm1
198             DO ji = 1, fs_jpim1   ! vector opt.
199                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
200                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
201             END DO
202          END DO
203       END SELECT
204       !                                                ! convert density into buoyancy
205       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
206       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
207       !
208       !
209       !                                      !==  Magnitude of the MLE stream function  ==!
210       !
211       !                 di[bm]  Ds
212       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
213       !                  e1u   Lf fu                      and the e2u for the "transport"
214       !                                                      (not *e3u as divided by e3u at the end)
215       !
216       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
217          DO jj = 1, jpjm1
218             DO ji = 1, fs_jpim1   ! vector opt.
219                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
220                     &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
221                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
222                !
223                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            &
224                     &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   &
225                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
226             END DO
227          END DO
228          !
229       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
230          DO jj = 1, jpjm1
231             DO ji = 1, fs_jpim1   ! vector opt.
232                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
233                     &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
234                !
235                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               &
236                     &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
237             END DO
238          END DO
239       ENDIF
240       !
241       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
242          DO jj = 1, jpjm1
243             DO ji = 1, fs_jpim1   ! vector opt.
244                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
245                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
246             END DO
247          END DO
248       ENDIF
249       !
250    ENDIF  ! end of lm_osm_mle loop
251    !                                      !==  structure function value at uw- and vw-points  ==!
252    DO jj = 1, jpjm1
253       DO ji = 1, fs_jpim1   ! vector opt.
254          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
255          zhv(ji,jj) = 1._wp / zhv(ji,jj)
256       END DO
257    END DO
258    !
259    zpsi_uw(:,:,:) = 0._wp
260    zpsi_vw(:,:,:) = 0._wp
261    !
262    DO jk = 2, ikmax                                ! start from 2 : surface value = 0
263       DO jj = 1, jpjm1
264          DO ji = 1, fs_jpim1   ! vector opt.
265             zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj)
266             zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj)
267             zcuw = zcuw * zcuw
268             zcvw = zcvw * zcvw
269             zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
270             zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
271             !
272             zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
273             zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
274          END DO
275       END DO
276    END DO
277    !
278    !                                      !==  transport increased by the MLE induced transport ==!
279    DO jk = 1, ikmax
280       DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
281          DO ji = 1, fs_jpim1   ! vector opt.
282             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
283             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
284          END DO
285       END DO
286       DO jj = 2, jpjm1
287          DO ji = fs_2, fs_jpim1   ! vector opt.
288             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
289                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
290          END DO
291       END DO
292    END DO
293
294    IF( cdtype == 'TRA') THEN              !==  outputs  ==!
295       !
296       IF (ln_osm_mle) THEN
297          zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
298       ELSE
299          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
300       END IF
301       CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
302       !
303       ! divide by cross distance to give streamfunction with dimensions m^2/s
304       DO jk = 1, ikmax+1
305          zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
306          zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
307       END DO
308       CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
309       CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
310    ENDIF
311    !
312  END SUBROUTINE tra_mle_trp
313
314
315   SUBROUTINE tra_mle_init
316      !!---------------------------------------------------------------------
317      !!                  ***  ROUTINE tra_mle_init  ***
318      !!
319      !! ** Purpose :   Control the consistency between namelist options for
320      !!              tracer advection schemes and set nadv
321      !!----------------------------------------------------------------------
322      INTEGER  ::   ji, jj, jk   ! dummy loop indices
323      INTEGER  ::   ierr
324      INTEGER ::    ios                 ! Local integer output status for namelist read
325      REAL(wp) ::   z1_t2, zfu, zfv                                !    -         -
326      !
327      NAMELIST/namtra_mle/ ln_mle , nn_mle, rn_ce, rn_lf, rn_time, rn_lat, nn_mld_uv, nn_conv, rn_rho_c_mle
328      !!----------------------------------------------------------------------
329
330      REWIND( numnam_ref )              ! Namelist namtra_mle in reference namelist : Tracer advection scheme
331      READ  ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901)
332901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist', lwp )
333
334      REWIND( numnam_cfg )              ! Namelist namtra_mle in configuration namelist : Tracer advection scheme
335      READ  ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 )
336902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist', lwp )
337      IF(lwm) WRITE ( numond, namtra_mle )
338
339      IF(lwp) THEN                     ! Namelist print
340         WRITE(numout,*)
341         WRITE(numout,*) 'tra_mle_init : mixed layer eddy (MLE) advection acting on tracers'
342         WRITE(numout,*) '~~~~~~~~~~~~~'
343         WRITE(numout,*) '   Namelist namtra_mle : mixed layer eddy advection applied on tracers'
344         WRITE(numout,*) '      use mixed layer eddy (MLE, i.e. Fox-Kemper param) (T/F)      ln_mle       = ', ln_mle
345         WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_mle    = ', nn_mle
346         WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_ce     = ', rn_ce
347         WRITE(numout,*) '         scale of ML front (ML radius of deformation) (rn_mle=0)      rn_lf     = ', rn_lf, 'm'
348         WRITE(numout,*) '         maximum time scale of MLE                    (rn_mle=0)      rn_time   = ', rn_time, 's'
349         WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (rn_mle=1)      rn_lat    = ', rn_lat, 'deg'
350         WRITE(numout,*) '         space interp. of MLD at u-(v-)pts (0=min,1=averaged,2=max)   nn_mld_uv = ', nn_mld_uv
351         WRITE(numout,*) '         =1 no MLE in case of convection ; =0 always MLE              nn_conv   = ', nn_conv
352         WRITE(numout,*) '         Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle
353      ENDIF
354
355
356      IF( ln_osm_mle .AND. ln_mle ) THEN
357         WRITE(numout,*) 'WARNING: You are running with both OSM-FK and default FK'
358         ! CALL ctl_stop('STOP in zdf_osm_init: Cannot run with both OSM-FK and default FK')
359     END IF
360
361     IF(lwp) THEN
362         WRITE(numout,*)
363         IF( ln_mle ) THEN
364            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to tracer advection'
365            IF( nn_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
366            IF( nn_mle == 1 )   WRITE(numout,*) '              New formulation'
367         ELSE
368            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy parametrisation NOT used'
369         ENDIF
370      ENDIF
371      !
372      IF( ln_mle ) THEN                ! MLE initialisation
373         !
374         rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria
375         IF(lwp) WRITE(numout,*)
376         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
377         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rho_c, 'kg/m3'
378         !
379         IF( nn_mle == 0 ) THEN           ! MLE array allocation & initialisation
380            ALLOCATE( rfu(jpi,jpj) , rfv(jpi,jpj) , STAT= ierr )
381            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
382            z1_t2 = 1._wp / ( rn_time * rn_time )
383            DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
384               DO ji = fs_2, jpi   ! vector opt.
385                  zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp
386                  zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp
387                  rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 )
388                  rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
389               END DO
390            END DO
391            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. )
392            !
393         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation
394            rc_f = rn_ce / (  5.e3_wp * 2._wp * omega * SIN( rad * rn_lat )  )
395            !
396         ENDIF
397         !
398         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_mle case)
399         ALLOCATE( r1_ft(jpi,jpj) , STAT= ierr )
400         IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate r1_ft array' )
401         !
402         z1_t2 = 1._wp / ( rn_time * rn_time )
403         r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
404         !
405      ENDIF
406      !
407   END SUBROUTINE tra_mle_init
408
409   !!==============================================================================
410END MODULE tramle
Note: See TracBrowser for help on using the repository browser.