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/2020/r12377_ticket2386/src/OCE/TRA – NEMO

source: NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/tramle.F90 @ 12451

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

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

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

svn merge --ignore-ancestry \

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

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

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

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