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/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/TRA/tramle.F90 @ 11954

Last change on this file since 11954 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 18.1 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 "vectopt_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 )
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  : - (pun,pvn,pwn) 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      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
83      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
84      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
85      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
86      !
87      INTEGER  ::   ji, jj, jk          ! dummy loop indices
88      INTEGER  ::   ii, ij, ik, ikmax   ! local integers
89      REAL(wp) ::   zcuw, zmuw, zc      ! local scalar
90      REAL(wp) ::   zcvw, zmvw          !   -      -
91      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
92      REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
93      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
94      !!----------------------------------------------------------------------
95      !
96      !                                      !==  MLD used for MLE  ==!
97      !                                                ! compute from the 10m density to deal with the diurnal cycle
98      inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
99      IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
100         DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m)
101            DO jj = 1, jpj
102               DO ji = 1, jpi                          ! index of the w-level at the ML based
103                  IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
104               END DO
105            END DO
106         END DO
107      ENDIF
108      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
109      !
110      !
111      zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
112      zbm (:,:) = 0._wp
113      zn2 (:,:) = 0._wp
114      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
115         DO jj = 1, jpj
116            DO ji = 1, jpi
117               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
118               zmld(ji,jj) = zmld(ji,jj) + zc
119               zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
120               zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
121            END DO
122         END DO
123      END DO
124
125      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
126      CASE ( 0 )                                               != min of the 2 neighbour MLDs
127         DO jj = 1, jpjm1
128            DO ji = 1, fs_jpim1   ! vector opt.
129               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
130               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
131            END DO
132         END DO
133      CASE ( 1 )                                               != average of the 2 neighbour MLDs
134         DO jj = 1, jpjm1
135            DO ji = 1, fs_jpim1   ! vector opt.
136               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
137               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
138            END DO
139         END DO
140      CASE ( 2 )                                               != max of the 2 neighbour MLDs
141         DO jj = 1, jpjm1
142            DO ji = 1, fs_jpim1   ! vector opt.
143               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
144               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
145            END DO
146         END DO
147      END SELECT
148      !                                                ! convert density into buoyancy
149      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
150      !
151      !
152      !                                      !==  Magnitude of the MLE stream function  ==!
153      !
154      !                 di[bm]  Ds
155      ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
156      !                  e1u   Lf fu                      and the e2u for the "transport"
157      !                                                      (not *e3u as divided by e3u at the end)
158      !
159      IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
160         DO jj = 1, jpjm1
161            DO ji = 1, fs_jpim1   ! vector opt.
162               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
163                  &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
164                  &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
165                  !
166               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            &
167                  &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   &
168                  &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
169            END DO
170         END DO
171         !
172      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
173         DO jj = 1, jpjm1
174            DO ji = 1, fs_jpim1   ! vector opt.
175               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
176                  &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
177                  !
178               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               &
179                  &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
180            END DO
181         END DO
182      ENDIF
183      !
184      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
185         DO jj = 1, jpjm1
186            DO ji = 1, fs_jpim1   ! vector opt.
187               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
188               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
189            END DO
190         END DO
191      ENDIF
192      !
193      !                                      !==  structure function value at uw- and vw-points  ==!
194      DO jj = 1, jpjm1
195         DO ji = 1, fs_jpim1   ! vector opt.
196            zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
197            zhv(ji,jj) = 1._wp / zhv(ji,jj)
198         END DO
199      END DO
200      !
201      zpsi_uw(:,:,:) = 0._wp
202      zpsi_vw(:,:,:) = 0._wp
203      !
204      DO jk = 2, ikmax                                ! start from 2 : surface value = 0
205         DO jj = 1, jpjm1
206            DO ji = 1, fs_jpim1   ! vector opt.
207               zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj)
208               zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj)
209               zcuw = zcuw * zcuw
210               zcvw = zcvw * zcvw
211               zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
212               zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
213               !
214               zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
215               zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
216            END DO
217         END DO
218      END DO
219      !
220      !                                      !==  transport increased by the MLE induced transport ==!
221      DO jk = 1, ikmax
222         DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
223            DO ji = 1, fs_jpim1   ! vector opt.
224               pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
225               pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
226            END DO
227         END DO
228         DO jj = 2, jpjm1
229            DO ji = fs_2, fs_jpim1   ! vector opt.
230               pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
231                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
232            END DO
233         END DO
234      END DO
235
236      IF( cdtype == 'TRA') THEN              !==  outputs  ==!
237         !
238         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
239         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
240         !
241         ! divide by cross distance to give streamfunction with dimensions m^2/s
242         DO jk = 1, ikmax+1
243            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
244            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
245         END DO
246         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
247         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
248      ENDIF
249      !
250   END SUBROUTINE tra_mle_trp
251
252
253   SUBROUTINE tra_mle_init
254      !!---------------------------------------------------------------------
255      !!                  ***  ROUTINE tra_mle_init  ***
256      !!
257      !! ** Purpose :   Control the consistency between namelist options for
258      !!              tracer advection schemes and set nadv
259      !!----------------------------------------------------------------------
260      INTEGER  ::   ji, jj, jk   ! dummy loop indices
261      INTEGER  ::   ierr
262      INTEGER ::    ios                 ! Local integer output status for namelist read
263      REAL(wp) ::   z1_t2, zfu, zfv                                !    -         -
264      !
265      NAMELIST/namtra_mle/ ln_mle , nn_mle, rn_ce, rn_lf, rn_time, rn_lat, nn_mld_uv, nn_conv, rn_rho_c_mle
266      !!----------------------------------------------------------------------
267
268      READ  ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901)
269901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist' )
270
271      READ  ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 )
272902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' )
273      IF(lwm) WRITE ( numond, namtra_mle )
274
275      IF(lwp) THEN                     ! Namelist print
276         WRITE(numout,*)
277         WRITE(numout,*) 'tra_mle_init : mixed layer eddy (MLE) advection acting on tracers'
278         WRITE(numout,*) '~~~~~~~~~~~~~'
279         WRITE(numout,*) '   Namelist namtra_mle : mixed layer eddy advection applied on tracers'
280         WRITE(numout,*) '      use mixed layer eddy (MLE, i.e. Fox-Kemper param) (T/F)      ln_mle       = ', ln_mle
281         WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_mle    = ', nn_mle
282         WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_ce     = ', rn_ce
283         WRITE(numout,*) '         scale of ML front (ML radius of deformation) (rn_mle=0)      rn_lf     = ', rn_lf, 'm'
284         WRITE(numout,*) '         maximum time scale of MLE                    (rn_mle=0)      rn_time   = ', rn_time, 's'
285         WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (rn_mle=1)      rn_lat    = ', rn_lat, 'deg'
286         WRITE(numout,*) '         space interp. of MLD at u-(v-)pts (0=min,1=averaged,2=max)   nn_mld_uv = ', nn_mld_uv
287         WRITE(numout,*) '         =1 no MLE in case of convection ; =0 always MLE              nn_conv   = ', nn_conv
288         WRITE(numout,*) '         Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle
289      ENDIF
290      !
291      IF(lwp) THEN
292         WRITE(numout,*)
293         IF( ln_mle ) THEN
294            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to tracer advection'
295            IF( nn_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
296            IF( nn_mle == 1 )   WRITE(numout,*) '              New formulation'
297         ELSE
298            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy parametrisation NOT used'
299         ENDIF
300      ENDIF
301      !
302      IF( ln_mle ) THEN                ! MLE initialisation
303         !
304         rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria
305         IF(lwp) WRITE(numout,*)
306         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
307         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rho_c, 'kg/m3'
308         !
309         IF( nn_mle == 0 ) THEN           ! MLE array allocation & initialisation
310            ALLOCATE( rfu(jpi,jpj) , rfv(jpi,jpj) , STAT= ierr )
311            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
312            z1_t2 = 1._wp / ( rn_time * rn_time )
313            DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
314               DO ji = fs_2, jpi   ! vector opt.
315                  zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp
316                  zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp
317                  rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 )
318                  rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
319               END DO
320            END DO
321            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. )
322            !
323         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation
324            rc_f = rn_ce / (  5.e3_wp * 2._wp * omega * SIN( rad * rn_lat )  )
325            !
326         ENDIF
327         !
328         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_mle case)
329         ALLOCATE( r1_ft(jpi,jpj) , STAT= ierr )
330         IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate r1_ft array' )
331         !
332         z1_t2 = 1._wp / ( rn_time * rn_time )
333         r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
334         !
335      ENDIF
336      !
337   END SUBROUTINE tra_mle_init
338
339   !!==============================================================================
340END MODULE tramle
Note: See TracBrowser for help on using the repository browser.