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 @ 13295

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

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision

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