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

source: NEMO/trunk/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
RevLine 
[9531]1MODULE tramle
[3959]2   !!======================================================================
[9531]3   !!                    ***  MODULE  tramle  ***
[3959]4   !! Ocean tracers: Mixed Layer Eddy induced transport
5   !!======================================================================
6   !! History :  3.3  !  2010-08  (G. Madec)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
[9531]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
[3959]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
[9019]17   !
[3959]18   USE in_out_manager ! I/O manager
19   USE iom            ! IOM library
20   USE lib_mpp        ! MPP library
[9124]21   USE lbclnk         ! lateral boundary condition / mpp link
[3959]22
23   IMPLICIT NONE
24   PRIVATE
25
[9531]26   PUBLIC   tra_mle_trp        ! routine called in traadv.F90
27   PUBLIC   tra_mle_init   ! routine called in traadv.F90
[3959]28
[9531]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
[5836]35   !                                        ! parameters used in nn_mle = 0 case
[9531]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
[5836]38   !                                        ! parameters used in nn_mle = 1 case
[9531]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
[3959]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
[12377]50#  include "do_loop_substitute.h90"
[3959]51   !!----------------------------------------------------------------------
[9598]52   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]53   !! $Id$
[10068]54   !! Software governed by the CeCILL license (see ./LICENSE)
[3959]55   !!----------------------------------------------------------------------
56CONTAINS
57
[12377]58   SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype, Kmm )
[3959]59      !!----------------------------------------------------------------------
[9531]60      !!                  ***  ROUTINE tra_mle_trp  ***
[3959]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      !!
[12377]73      !! ** Action  : - (pu,pv,pw) increased by the mle transport
[3959]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
[12377]82      INTEGER                         , INTENT(in   ) ::   Kmm        ! ocean time level index
[3959]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      !
[9019]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
[3959]95      !!----------------------------------------------------------------------
[5836]96      !
[3995]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)
[10105]100      IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
[12377]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
[10105]104      ENDIF
[4325]105      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
[3959]106      !
107      !
[3995]108      zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
109      zbm (:,:) = 0._wp
110      zn2 (:,:) = 0._wp
[12377]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
[3959]117
[3995]118      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
[3959]119      CASE ( 0 )                                               != min of the 2 neighbour MLDs
[12377]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
[3959]124      CASE ( 1 )                                               != average of the 2 neighbour MLDs
[12377]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
[3959]129      CASE ( 2 )                                               != max of the 2 neighbour MLDs
[12377]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
[3959]134      END SELECT
[3995]135      !                                                ! convert density into buoyancy
[12377]136      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )
[3959]137      !
138      !
[3995]139      !                                      !==  Magnitude of the MLE stream function  ==!
140      !
[3959]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
[12377]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
[3959]156         !
157      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
[12377]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
[3959]165      ENDIF
166      !
167      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
[12377]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
[3959]172      ENDIF
173      !
[3995]174      !                                      !==  structure function value at uw- and vw-points  ==!
[12377]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
[4835]179      !
[3959]180      zpsi_uw(:,:,:) = 0._wp
181      zpsi_vw(:,:,:) = 0._wp
182      !
[12377]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
[3995]194      !
195      !                                      !==  transport increased by the MLE induced transport ==!
[3959]196      DO jk = 1, ikmax
[12377]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
[3959]205      END DO
206
[3995]207      IF( cdtype == 'TRA') THEN              !==  outputs  ==!
[3959]208         !
209         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
[3995]210         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
[3959]211         !
[4188]212         ! divide by cross distance to give streamfunction with dimensions m^2/s
213         DO jk = 1, ikmax+1
[5836]214            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
215            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
[4188]216         END DO
[3959]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      !
[9531]221   END SUBROUTINE tra_mle_trp
[3959]222
223
[9531]224   SUBROUTINE tra_mle_init
[3959]225      !!---------------------------------------------------------------------
[9531]226      !!                  ***  ROUTINE tra_mle_init  ***
[3959]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
[4245]233      INTEGER ::    ios                 ! Local integer output status for namelist read
[3959]234      REAL(wp) ::   z1_t2, zfu, zfv                                !    -         -
235      !
[9531]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
[3959]237      !!----------------------------------------------------------------------
238
[9531]239      READ  ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901)
[11536]240901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist' )
[4245]241
[9531]242      READ  ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 )
[11536]243902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' )
[9531]244      IF(lwm) WRITE ( numond, namtra_mle )
[4245]245
[3959]246      IF(lwp) THEN                     ! Namelist print
247         WRITE(numout,*)
[9531]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'
[9168]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
[3959]260      ENDIF
261      !
262      IF(lwp) THEN
263         WRITE(numout,*)
264         IF( ln_mle ) THEN
[9190]265            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to tracer advection'
[7646]266            IF( nn_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
267            IF( nn_mle == 1 )   WRITE(numout,*) '              New formulation'
[3959]268         ELSE
[9190]269            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy parametrisation NOT used'
[3959]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 )
[12377]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
[10425]290            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. )
[3959]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 )
[7753]302         r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
[3959]303         !
304      ENDIF
305      !
[9531]306   END SUBROUTINE tra_mle_init
[3959]307
308   !!==============================================================================
[9531]309END MODULE tramle
Note: See TracBrowser for help on using the repository browser.