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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
File size: 19.0 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              ! flag to activate the Mixed Layer Eddy (MLE) parameterisation
32   INTEGER         ::   nn_mle              ! MLE type: =0 standard Fox-Kemper ; =1 new formulation
33   INTEGER         ::   nn_mld_uv           ! space interpolation of MLD at u- & v-pts (0=min,1=averaged,2=max)
34   INTEGER         ::   nn_conv             ! =1 no MLE in case of convection ; =0 always MLE
35   REAL(wp)        ::   rn_ce               ! MLE coefficient
36   !                                           ! parameters used in nn_mle = 0 case
37   REAL(wp)        ::   rn_lf                  ! typical scale of mixed layer front
38   REAL(wp)        ::   rn_time             ! time scale for mixing momentum across the mixed layer
39   !                                             ! parameters used in nn_mle = 1 case
40   REAL(wp)        ::   rn_lat                   ! reference latitude for a 5 km scale of ML front
41   REAL(wp)        ::   rn_rho_c_mle         ! 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      !
109      !                                      !==  MLD used for MLE  ==!
110      !                                                ! compute from the 10m density to deal with the diurnal cycle
111      inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
112      DO jk = jpkm1, nlb10, -1                         ! from the bottom to nlb10 (10m)
113         DO jj = 1, jpj
114            DO ji = 1, jpi                             ! index of the w-level at the ML based
115               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
116            END DO
117         END DO
118      END DO
119      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
120      !
121      !
122      zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
123      zbm (:,:) = 0._wp
124      zn2 (:,:) = 0._wp
125      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
126         DO jj = 1, jpj
127            DO ji = 1, jpi
128               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
129               zmld(ji,jj) = zmld(ji,jj) + zc
130               zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
131               zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
132            END DO
133         END DO
134      END DO
135
136      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
137      CASE ( 0 )                                               != min of the 2 neighbour MLDs
138         DO jj = 1, jpjm1
139            DO ji = 1, fs_jpim1   ! vector opt.
140               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
141               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
142            END DO
143         END DO
144      CASE ( 1 )                                               != average of the 2 neighbour MLDs
145         DO jj = 1, jpjm1
146            DO ji = 1, fs_jpim1   ! vector opt.
147               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
148               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
149            END DO
150         END DO
151      CASE ( 2 )                                               != max of the 2 neighbour MLDs
152         DO jj = 1, jpjm1
153            DO ji = 1, fs_jpim1   ! vector opt.
154               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
155               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
156            END DO
157         END DO
158      END SELECT
159      !                                                ! convert density into buoyancy
160      zbm(:,:) = + grav * zbm(:,:) / MAX( fse3t(:,:,1), zmld(:,:) )
161      !
162      !
163      !                                      !==  Magnitude of the MLE stream function  ==!
164      !
165      !                 di[bm]  Ds
166      ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
167      !                  e1u   Lf fu                      and the e2u for the "transport"
168      !                                                      (not *e3u as divided by e3u at the end)
169      !
170      IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
171         DO jj = 1, jpjm1
172            DO ji = 1, fs_jpim1   ! vector opt.
173               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
174                  &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp          , e1u(ji,jj)                 )   &
175                  &           / (         e1u(ji,jj)          * MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
176                  !
177               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            &
178                  &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp          , e2v(ji,jj)                 )   &
179                  &           / (         e2v(ji,jj)          * MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
180            END DO
181         END DO
182         !
183      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
184         DO jj = 1, jpjm1
185            DO ji = 1, fs_jpim1   ! vector opt.
186               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj) / e1u(ji,jj)          &
187                  &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
188                  !
189               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj) / e2v(ji,jj)          &
190                  &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
191            END DO
192         END DO
193      ENDIF
194      !
195      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
196         DO jj = 1, jpjm1
197            DO ji = 1, fs_jpim1   ! vector opt.
198               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
199               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
200            END DO
201         END DO
202      ENDIF
203      !
204      !                                      !==  structure function value at uw- and vw-points  ==!
205      zhu(:,:) = 1._wp / zhu(:,:)                   ! hu --> 1/hu
206      zhv(:,:) = 1._wp / zhv(:,:)
207      zpsi_uw(:,:,:) = 0._wp
208      zpsi_vw(:,:,:) = 0._wp
209      !
210      DO jk = 2, ikmax                                ! start from 2 : surface value = 0
211         DO jj = 1, jpjm1
212            DO ji = 1, fs_jpim1   ! vector opt.
213               zcuw = 1._wp - ( fsdepw(ji+1,jj,jk) + fsdepw(ji,jj,jk) ) * zhu(ji,jj)
214               zcvw = 1._wp - ( fsdepw(ji,jj+1,jk) + fsdepw(ji,jj,jk) ) * zhv(ji,jj)
215               zcuw = zcuw * zcuw
216               zcvw = zcvw * zcvw
217               zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
218               zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
219               !
220               zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
221               zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
222            END DO
223         END DO
224      END DO
225      !
226      !                                      !==  transport increased by the MLE induced transport ==!
227      DO jk = 1, ikmax
228         DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
229            DO ji = 1, fs_jpim1   ! vector opt.
230               pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
231               pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
232            END DO
233         END DO
234         DO jj = 2, jpjm1
235            DO ji = fs_2, fs_jpim1   ! vector opt.
236               pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
237                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
238            END DO
239         END DO
240      END DO
241
242      IF( cdtype == 'TRA') THEN              !==  outputs  ==!
243         !
244         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
245         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
246         !
247         ! divide by cross distance to give streamfunction with dimensions m^2/s
248         DO jk = 1, ikmax+1
249            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk)/e2u(:,:)
250            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk)/e1v(:,:)
251         END DO
252         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
253         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
254      ENDIF
255      CALL wrk_dealloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH)
256      CALL wrk_dealloc( jpi, jpj, jpk, zpsi_uw, zpsi_vw)
257      CALL wrk_dealloc( jpi, jpj, inml_mle)
258
259      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mle')
260      !
261   END SUBROUTINE tra_adv_mle
262
263
264   SUBROUTINE tra_adv_mle_init
265      !!---------------------------------------------------------------------
266      !!                  ***  ROUTINE tra_adv_mle_init  ***
267      !!
268      !! ** Purpose :   Control the consistency between namelist options for
269      !!              tracer advection schemes and set nadv
270      !!----------------------------------------------------------------------
271      INTEGER  ::   ji, jj, jk   ! dummy loop indices
272      INTEGER  ::   ierr
273      INTEGER ::    ios                 ! Local integer output status for namelist read
274      REAL(wp) ::   z1_t2, zfu, zfv                                !    -         -
275      !
276      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
277      !!----------------------------------------------------------------------
278
279
280      REWIND( numnam_ref )              ! Namelist namtra_adv_mle in reference namelist : Tracer advection scheme
281      READ  ( numnam_ref, namtra_adv_mle, IOSTAT = ios, ERR = 901)
282901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv_mle in reference namelist', lwp )
283
284      REWIND( numnam_cfg )              ! Namelist namtra_adv_mle in configuration namelist : Tracer advection scheme
285      READ  ( numnam_cfg, namtra_adv_mle, IOSTAT = ios, ERR = 902 )
286902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv_mle in configuration namelist', lwp )
287      IF(lwm) WRITE ( numond, namtra_adv_mle )
288
289      IF(lwp) THEN                     ! Namelist print
290         WRITE(numout,*)
291         WRITE(numout,*) 'tra_adv_mle_init : mixed layer eddy (MLE) advection acting on tracers'
292         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
293         WRITE(numout,*) '   Namelist namtra_adv_mle : mixed layer eddy advection on tracers'
294         WRITE(numout,*) '      use mixed layer eddy (MLE, i.e. Fox-Kemper param) (T/F)      ln_mle    = ', ln_mle
295         WRITE(numout,*) '      MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_mle    = ', nn_mle
296         WRITE(numout,*) '      magnitude of the MLE (typical value: 0.06 to 0.08)           rn_ce     = ', rn_ce
297         WRITE(numout,*) '      scale of ML front (ML radius of deformation) (rn_mle=0)      rn_lf     = ', rn_lf, 'm'
298         WRITE(numout,*) '      maximum time scale of MLE                    (rn_mle=0)      rn_time   = ', rn_time, 's'
299         WRITE(numout,*) '      reference latitude (degrees) of MLE coef.    (rn_mle=1)      rn_lat    = ', rn_lat, 'deg'
300         WRITE(numout,*) '      space interp. of MLD at u-(v-)pts (0=min,1=averaged,2=max)   nn_mld_uv = ', nn_mld_uv
301         WRITE(numout,*) '      =1 no MLE in case of convection ; =0 always MLE              nn_conv   = ', nn_conv
302         WRITE(numout,*) '      Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle
303      ENDIF
304      !
305      IF(lwp) THEN
306         WRITE(numout,*)
307         IF( ln_mle ) THEN
308            WRITE(numout,*) '   Mixed Layer Eddy induced transport added to tracer advection'
309            IF( nn_mle == 0 )   WRITE(numout,*) '   Fox-Kemper et al 2010 formulation'
310            IF( nn_mle == 1 )   WRITE(numout,*) '   New formulation'
311         ELSE
312            WRITE(numout,*) '   Mixed Layer Eddy parametrisation NOT used'
313         ENDIF
314      ENDIF
315      !
316      IF( ln_mle ) THEN                ! MLE initialisation
317         !
318         rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria
319         IF(lwp) WRITE(numout,*)
320         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
321         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rho_c, 'kg/m3'
322         !
323         IF( nn_mle == 0 ) THEN           ! MLE array allocation & initialisation
324            ALLOCATE( rfu(jpi,jpj) , rfv(jpi,jpj) , STAT= ierr )
325            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
326            z1_t2 = 1._wp / ( rn_time * rn_time )
327            DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
328               DO ji = fs_2, jpi   ! vector opt.
329                  zfu = ( ff(ji,jj) + ff(ji,jj-1) ) * 0.5_wp
330                  zfv = ( ff(ji,jj) + ff(ji-1,jj) ) * 0.5_wp
331                  rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 )
332                  rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
333               END DO
334            END DO
335            CALL lbc_lnk( rfu, 'U', 1. )   ;   CALL lbc_lnk( rfv, 'V', 1. )
336            !
337         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation
338            rc_f = rn_ce / (  5.e3_wp * 2._wp * omega * SIN( rad * rn_lat )  )
339            !
340         ENDIF
341         !
342         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_mle case)
343         ALLOCATE( r1_ft(jpi,jpj) , STAT= ierr )
344         IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate r1_ft array' )
345         !
346         z1_t2 = 1._wp / ( rn_time * rn_time )
347         r1_ft(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )
348         r1_ft(:,:) = 1._wp / SQRT(  r1_ft(:,:) * r1_ft(:,:) + z1_t2 )
349         !
350      ENDIF
351      !
352   END SUBROUTINE tra_adv_mle_init
353
354   !!==============================================================================
355END MODULE traadv_mle
Note: See TracBrowser for help on using the repository browser.