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.
traadv_mle.F90 in branches/2013/dev_r3948_NOC_FK/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2013/dev_r3948_NOC_FK/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90 @ 3959

Last change on this file since 3959 was 3959, checked in by agn, 11 years ago

Gurvan's FK + separate rn_rho_c_mle

File size: 21.8 KB
Line 
1MODULE traadv_mle
2   !!======================================================================
3   !!                    ***  MODULE  traadv_mle  ***
4   !! Ocean tracers: Mixed Layer Eddy induced transport
5   !!======================================================================
6   !! History :  3.3  !  2010-08  (G. Madec)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_adv_mle      : update the effective transport with the Mixed Layer Eddy induced transport
11   !!   tra_adv_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   USE lbclnk         ! lateral boundary condition / mpp link
18   USE in_out_manager ! I/O manager
19   USE iom            ! IOM library
20   USE lib_mpp        ! MPP library
21   USE wrk_nemo       ! work arrays
22   USE timing         ! Timing
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   tra_adv_mle        ! routine called in traadv.F90
28   PUBLIC   tra_adv_mle_init   ! routine called in traadv.F90
29
30   !                                               !!* namelist namtra_adv_mle *
31   LOGICAL, PUBLIC ::   ln_mle    = .TRUE.          ! flag to activate the Mixed Layer Eddy (MLE) parameterisation
32   INTEGER         ::   nn_mle    = 0               ! MLE type: =0 standard Fox-Kemper ; =1 new formulation
33   INTEGER         ::   nn_mld_uv = 0               ! space interpolation of MLD at u- & v-pts (0=min,1=averaged,2=max)
34   INTEGER         ::   nn_conv   = 0               ! =1 no MLE in case of convection ; =0 always MLE
35   REAL(wp)        ::   rn_ce   = 0.06_wp           ! MLE coefficient
36   !                                                ! parameters used in nn_mle = 0 case
37   REAL(wp)        ::   rn_lf   = 5.e+3_wp               ! typical scale of mixed layer front
38   REAL(wp)        ::   rn_time = 2._wp * 86400._wp      ! time scale for mixing momentum across the mixed layer
39   !                                                ! parameters used in nn_mle = 1 case
40   REAL(wp)        ::   rn_lat  = 20._wp                 ! reference latitude for a 5 km scale of ML front
41   REAL(wp)        ::   rn_rho_c_mle  = 0.01        ! Density criterion for definition of MLD used by FK
42
43   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation
44   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
45   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case
46
47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rfu, rfv   ! modified Coriolis parameter (f+tau) at u- & v-pts
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_ft      ! inverse of the modified Coriolis parameter at t-pts
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
55   !! $Id:$
56   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE tra_adv_mle( kt, kit000, pu, pv, pw, cdtype )
61      !!----------------------------------------------------------------------
62      !!                  ***  ROUTINE adv_mle  ***
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      !
83      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
84      INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index
85      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
86      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
87      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
88      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
89      !
90      INTEGER  ::   ji, jj, jk   ! dummy loop indices
91      INTEGER  ::   ikmax        ! temporary integer
92      REAL(wp) ::   zcuw, zmuw   ! local scalar
93      REAL(wp) ::   zcvw, zmvw   !   -      -
94      REAL(wp) ::   zc                                     !   -      -
95
96      INTEGER  ::   ii, ij, ik              ! local integers
97      INTEGER, DIMENSION(3) ::   ilocu      !
98      INTEGER, DIMENSION(2) ::   ilocs      !
99      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
100      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpsi_uw, zpsi_vw
101      INTEGER, POINTER, DIMENSION(:,:) :: inml_mle
102      !!----------------------------------------------------------------------
103
104      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mle')
105      CALL wrk_alloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH)
106      CALL wrk_alloc( jpi, jpj, jpk, zpsi_uw, zpsi_vw)
107      CALL wrk_alloc( jpi, jpj, inml_mle)
108!!gm   Caution HERE : * compute the nmln from the 10m density to deal with the diurnal cycle
109!!gm                  * check that nmld is the number of T-level not w-levels...
110
111      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
112         DO jj = 1, jpj
113            DO ji = 1, jpi
114               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
115            END DO
116         END DO
117      END DO
118      ikmax = MAXVAL( inml_mle(:,:) )                  ! max level of the computation
119      !
120      zbm(:,:) = 0._wp   ;   zmld(:,:) = 0._wp; zn2(:,:) = 0._wp   ! Horizontal shape of the MLE
121      !
122
123
124   IF(lwp) WRITE(numout,*) 'mle 0'
125
126           ilocs = MAXLOC( inml_mle(:,:) - 1 )
127           ii = ilocs(1)
128           ij = ilocs(2)
129            ! ii=120
130            ! ij=110
131             !  inml_mle(ii,ij) = Max( 4, inml_mle(ii,ij) )
132            WRITE(numout,*) 'mle zmld max', ii,ij, inml_mle(ii,ij), 'ikmax', ikmax
133
134
135      DO jk = 1, ikmax                             ! MLD and mean buoyancy and N2 over the mixed layer
136   IF(lwp) WRITE(numout,*) 'mle zc',jk,' c ', REAL( MIN( MAX( 0, inml_mle(ii,ij)-jk ) , 1  )  )
137         DO jj = 1, jpj
138            DO ji = 1, jpi
139               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! e3t being 0 outside the ML
140               zmld(ji,jj) = zmld(ji,jj) + zc
141               zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) / rau0
142               zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
143            END DO
144         END DO
145      END DO
146
147      SELECT CASE( nn_mld_uv )                           ! MLD at u- & v-pts
148      CASE ( 0 )                                               != min of the 2 neighbour MLDs
149         DO jj = 1, jpjm1
150            DO ji = 1, fs_jpim1   ! vector opt.
151               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
152               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
153            END DO
154         END DO
155      CASE ( 1 )                                               != average of the 2 neighbour MLDs
156         DO jj = 1, jpjm1
157            DO ji = 1, fs_jpim1   ! vector opt.
158               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
159               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
160            END DO
161         END DO
162      CASE ( 2 )                                               != max of the 2 neighbour MLDs
163         DO jj = 1, jpjm1
164            DO ji = 1, fs_jpim1   ! vector opt.
165               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
166               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
167            END DO
168         END DO
169      END SELECT
170
171      zbm(:,:) = + grav * zbm(:,:) / MAX( fse3t(:,:,1), zmld(:,:) )
172
173
174   write(numout,*) 'n2: max ', MAXVAL(zn2(:,:) , MASK = tmask(:,:,1) ==1. ),    &
175         &            ' min ', MINVAL(zn2(:,:) , MASK = tmask(:,:,1) ==1. )
176   write(numout,*) 'h: max ', MAXVAL(zhu(1:jpim1,1:jpjm1) ), ' min ', MINVAL(zhu(1:jpim1,1:jpjm1) )
177   write(numout,*) 'h: max ', MAXVAL(zhv(1:jpim1,1:jpjm1) ), ' min ', MINVAL(zhv(1:jpim1,1:jpjm1) )
178
179      IF(lwp) write(numout,*) 'mle: max zbm', MAXVAL(-zbm(:,:) , MASK = tmask(:,:,1) ==1. ) / grav,    &
180         &                    ' min '       , MINVAL(-zbm(:,:) , MASK = tmask(:,:,1) ==1. ) / grav
181
182
183      !
184      ! Magnitude of the MLE stream function:
185      !
186      !                 di[bm]  Ds
187      ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
188      !                  e1u   Lf fu                      and the e2u for the "transport"
189      !                                                      (not *e3u as divided by e3u at the end)
190! MLD criteria set in zdfmxl: rho_c = 0.01 kg/m3
191      zpsim_u(:,:) = 0._wp
192      zpsim_v(:,:) = 0._wp
193      !
194
195      IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
196         DO jj = 1, jpjm1
197            DO ji = 1, fs_jpim1   ! vector opt.
198               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
199                  &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp          , e1u(ji,jj)                 )   &
200                  &           / (         e1u(ji,jj)          * MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
201                  !
202               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            &
203                  &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp          , e2v(ji,jj)                 )   &
204                  &           / (         e1u(ji,jj)          * MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
205            END DO
206         END DO
207         !
208      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
209         DO jj = 1, jpjm1
210            DO ji = 1, fs_jpim1   ! vector opt.
211               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj) / e1u(ji,jj)          &
212                  &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
213                  !
214               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj) / e1u(ji,jj)          &
215                  &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
216            END DO
217         END DO
218      ENDIF
219      !
220      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
221         DO jj = 1, jpjm1
222            DO ji = 1, fs_jpim1   ! vector opt.
223               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
224               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
225            END DO
226         END DO
227      ENDIF
228      !
229      !
230      !                                        ! structure function value at uw- and vw-points
231      zhu(:,:) = 1._wp / zhu(:,:)                   ! hu --> 1/hu
232      zhv(:,:) = 1._wp / zhv(:,:)
233      zpsi_uw(:,:,:) = 0._wp
234      zpsi_vw(:,:,:) = 0._wp
235      !
236      DO jk = 2, ikmax                                ! start from 2 : surface value = 0
237         DO jj = 1, jpjm1
238            DO ji = 1, fs_jpim1   ! vector opt.
239               zcuw = 1._wp - ( fsdepw(ji+1,jj,jk) + fsdepw(ji,jj,jk) ) * zhu(ji,jj)
240               zcvw = 1._wp - ( fsdepw(ji,jj+1,jk) + fsdepw(ji,jj,jk) ) * zhv(ji,jj)
241               zcuw = zcuw * zcuw
242               zcvw = zcvw * zcvw
243               zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
244               zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
245               !
246               zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
247               zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
248            END DO
249         END DO
250      END DO
251
252!!!gm start
253   write(numout,*) 'psi max u ', MAXVAL(ABS(zpsim_u(1:jpim1,1:jpjm1)*umask(1:jpim1,1:jpjm1,1) )),    &
254         &                ' v ', MAXVAL(ABS(zpsim_v(1:jpim1,1:jpjm1)*vmask(1:jpim1,1:jpjm1,1) ))
255
256            ilocs = MAXLOC( ABS(zpsim_v(1:jpim1,1:jpjm1)*vmask(1:jpim1,1:jpjm1,1) )) !, inml_mle(1:jpim1,1:jpjm1 )
257            ii = ilocs(1)
258            ij = ilocs(2)
259               WRITE(numout,*)        'hu, hv'      , 1./zhu(ii,ij), 1./zhv(ii,ij)
260               WRITE(numout,*)
261               WRITE(numout,*)        'psim at i,j = ', ii, ij, zpsim_u(ii,ij)
262!           ji=ii
263!           jj=ij
264
265      DO jk = 2, ikmax                                ! start from 2 : surface value = 0
266               zcuw = 1.e0 - ( fsdepw(ji+1,jj,jk) + fsdepw(ji,jj,jk) ) * zhu(ji,jj)
267               zcvw = 1.e0 - ( fsdepw(ji,jj+1,jk) + fsdepw(ji,jj,jk) ) * zhv(ji,jj)
268               zcuw = zcuw * zcuw
269               zcvw = zcvw * zcvw
270               zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
271               zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
272               !
273               WRITE(numout,*)        jk, zmuw, zmvw
274      END DO
275
276   IF(lwp) WRITE(numout,*) 'mle d'
277      !!!gm
278      IF(lwp) write(numout,*) 'max uw', MAXVAL(zpsi_uw,MASK = umask==1.), 'vw', MAXVAL(zpsi_vw,MASK = vmask==1.)
279      IF(lwp) write(numout,*) 'min uw', MinVAL(zpsi_uw,MASK = umask==1.), 'vw', MinVAL(zpsi_vw,MASK = vmask==1.)
280      !!!gm
281
282      IF(lwp) write(numout,*) 'rho i i1', rhop(ii,ij,1), rhop(ii+1,ij,1)
283      IF(lwp) write(numout,*) 'zbm i i1', zbm(ii,ij), zbm(ii+1,ij)
284      IF(lwp) write(numout,*) 'Dbm i i1', zbm(ii+1,ij)-zbm(ii,ij)
285      IF(lwp) write(numout,*) 'zmld i i1', zmld(ii,ij), zmld(ii+1,ij)
286      IF(lwp) write(numout,*) 'Dmld i i1', zmld(ii+1,ij)-zmld(ii,ij)
287      IF(lwp) write(numout,*) 'hu   i   ', zhu(ii,ij)
288      IF(lwp) write(numout,*) 'mld  i i1', inml_mle(ii,ij), inml_mle(ii+1,ij)
289
290      IF(lwp) write(numout,*) 'rho j j1', rhop(ii,ij,1), rhop(ii,ij+1,1)
291      IF(lwp) write(numout,*) 'zbm j j1', zbm(ii,ij), zbm(ii,ij+1)
292      IF(lwp) write(numout,*) 'Dbm j j1', zbm(ii,ij+1)-zbm(ii,ij)
293      IF(lwp) write(numout,*) 'zmld j j1', zmld(ii,ij), zmld(ii,ij+1)
294      IF(lwp) write(numout,*) 'Dmld j j1', zmld(ii,ij+1)-zmld(ii,ij)
295      IF(lwp) write(numout,*) 'hv   i   ', zhv(ii,ij)
296      IF(lwp) write(numout,*) 'mld  i i1', inml_mle(ii,ij), inml_mle(ii,ij+1)
297
298      WRITE(numout,*) 'mle zpsi: i,j,k',ii,ij, 'mask', tmask(ii,ij,1), 'mld', zmld(ii,ij)
299      DO jk = 1, jpk
300         WRITE(numout,*) jk, zpsi_uw(ii,ij,jk), ' u - v', zpsi_vw(ii,ij,jk)
301      END DO
302!!!gm end
303
304
305
306
307
308      DO jk = 1, ikmax
309         DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
310            DO ji = 1, fs_jpim1   ! vector opt.
311               pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
312               pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
313            END DO
314         END DO
315         DO jj = 2, jpjm1
316            DO ji = fs_2, fs_jpim1   ! vector opt.
317               pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
318                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
319            END DO
320         END DO
321      END DO
322
323!!!gm start
324      IF(lwp) THEN
325         write(numout,*) jk, 'velocity min/max'
326         Do jk=1,jpk
327         write(numout,*) jk, 'uw', MAXVAL(ABS((zpsi_uw(:,:,jk)-zpsi_uw(:,:,jk+1))/( e2u(:,:)*fse3u(:,:,jk)))),  &
328            &                'vw', MAXVAL(ABS((zpsi_vw(:,:,jk)-zpsi_vw(:,:,jk+1))/( e1v(:,:)*fse3v(:,:,jk))))
329         end do
330      endif
331!!!gm end
332
333      IF( cdtype == 'TRA') THEN
334         !
335         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
336         !
337         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
338         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
339         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
340      ENDIF
341      CALL wrk_dealloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH)
342      CALL wrk_dealloc( jpi, jpj, jpk, zpsi_uw, zpsi_vw)
343      CALL wrk_dealloc( jpi, jpj, inml_mle)
344
345      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mle')
346      !
347   END SUBROUTINE tra_adv_mle
348
349
350   SUBROUTINE tra_adv_mle_init
351      !!---------------------------------------------------------------------
352      !!                  ***  ROUTINE tra_adv_mle_init  ***
353      !!
354      !! ** Purpose :   Control the consistency between namelist options for
355      !!              tracer advection schemes and set nadv
356      !!----------------------------------------------------------------------
357      INTEGER  ::   ji, jj, jk   ! dummy loop indices
358      INTEGER  ::   ierr
359      REAL(wp) ::   z1_t2, zfu, zfv                                !    -         -
360      !
361      NAMELIST/namtra_adv_mle/ ln_mle , nn_mle, rn_ce, rn_lf, rn_time, rn_lat, nn_mld_uv, nn_conv, rn_rho_c_mle
362      !!----------------------------------------------------------------------
363
364      REWIND ( numnam )                ! Read Namelist namtra_adv_mle : mixed layer eddy advection acting on tracers
365      READ   ( numnam, namtra_adv_mle )
366
367      IF(lwp) THEN                     ! Namelist print
368         WRITE(numout,*)
369         WRITE(numout,*) 'tra_adv_mle_init : mixed layer eddy (MLE) advection acting on tracers'
370         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
371         WRITE(numout,*) '   Namelist namtra_adv_mle : mixed layer eddy advection on tracers'
372         WRITE(numout,*) '      use mixed layer eddy (MLE, i.e. Fox-Kemper param) (T/F)      ln_mle    = ', ln_mle
373         WRITE(numout,*) '      MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_mle    = ', nn_mle
374         WRITE(numout,*) '      magnitude of the MLE (typical value: 0.06 to 0.08)           rn_ce     = ', rn_ce
375         WRITE(numout,*) '      scale of ML front (ML radius of deformation) (rn_mle=0)      rn_lf     = ', rn_lf, 'm'
376         WRITE(numout,*) '      maximum time scale of MLE                    (rn_mle=0)      rn_time   = ', rn_time, 's'
377         WRITE(numout,*) '      reference latitude (degrees) of MLE coef.    (rn_mle=1)      rn_lat    = ', rn_lat, 'deg'
378         WRITE(numout,*) '      space interp. of MLD at u-(v-)pts (0=min,1=averaged,2=max)   nn_mld_uv = ', nn_mld_uv
379         WRITE(numout,*) '      =1 no MLE in case of convection ; =0 always MLE              nn_conv   = ', nn_conv
380         WRITE(numout,*) '      Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle
381      ENDIF
382      !
383      IF(lwp) THEN
384         WRITE(numout,*)
385         IF( ln_mle ) THEN
386            WRITE(numout,*) '   Mixed Layer Eddy induced transport added to tracer advection'
387            IF( nn_mle == 0 )   WRITE(numout,*) '   Fox-Kemper et al 2010 formulation'
388            IF( nn_mle == 1 )   WRITE(numout,*) '   New formulation'
389         ELSE
390            WRITE(numout,*) '   Mixed Layer Eddy parametrisation NOT used'
391         ENDIF
392      ENDIF
393      !
394      IF( ln_mle ) THEN                ! MLE initialisation
395         !
396         rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria
397         IF(lwp) WRITE(numout,*)
398         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
399         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rho_c, 'kg/m3'
400         !
401         IF( nn_mle == 0 ) THEN           ! MLE array allocation & initialisation
402            ALLOCATE( rfu(jpi,jpj) , rfv(jpi,jpj) , STAT= ierr )
403            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
404            z1_t2 = 1._wp / ( rn_time * rn_time )
405            DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
406               DO ji = fs_2, jpi   ! vector opt.
407                  zfu = ( ff(ji,jj) + ff(ji,jj-1) ) * 0.5_wp
408                  zfv = ( ff(ji,jj) + ff(ji-1,jj) ) * 0.5_wp
409                  rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 )
410                  rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
411               END DO
412            END DO
413            CALL lbc_lnk( rfu, 'U', 1. )   ;   CALL lbc_lnk( rfv, 'V', 1. )
414            !
415         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation
416            rc_f = rn_ce / (  5.e3_wp * 2._wp * omega * SIN( rad * rn_lat )  )
417            !
418         ENDIF
419         !
420         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_mle case)
421         ALLOCATE( r1_ft(jpi,jpj) , STAT= ierr )
422         IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate r1_ft array' )
423         !
424         z1_t2 = 1._wp / ( rn_time * rn_time )
425         r1_ft(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )
426         r1_ft(:,:) = 1._wp / SQRT(  r1_ft(:,:) * r1_ft(:,:) + z1_t2 )
427         !
428      ENDIF
429      !
430   END SUBROUTINE tra_adv_mle_init
431
432   !!==============================================================================
433END MODULE traadv_mle
Note: See TracBrowser for help on using the repository browser.