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/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/tramle.F90 @ 14576

Last change on this file since 14576 was 14576, checked in by hadcv, 3 years ago

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