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 branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90 @ 13320

Last change on this file since 13320 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 19.3 KB
RevLine 
[3959]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 *
[4245]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
[3959]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)
[5215]55   !! $Id$
[3959]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)
[3995]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)
[3959]113         DO jj = 1, jpj
[3995]114            DO ji = 1, jpi                             ! index of the w-level at the ML based
[3959]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
[4325]119      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
[3959]120      !
121      !
[3995]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
[3959]126         DO jj = 1, jpj
127            DO ji = 1, jpi
[3995]128               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
[3959]129               zmld(ji,jj) = zmld(ji,jj) + zc
[3995]130               zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
[3959]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
[3995]136      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
[3959]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
[3995]159      !                                                ! convert density into buoyancy
[3959]160      zbm(:,:) = + grav * zbm(:,:) / MAX( fse3t(:,:,1), zmld(:,:) )
161      !
162      !
[3995]163      !                                      !==  Magnitude of the MLE stream function  ==!
164      !
[3959]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)                 )   &
[3995]179                  &           / (         e2v(ji,jj)          * MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
[3959]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                  !
[3995]189               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj) / e2v(ji,jj)          &
[3959]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      !
[3995]204      !                                      !==  structure function value at uw- and vw-points  ==!
[4835]205      DO jj = 1, jpjm1
206         DO ji = 1, fs_jpim1   ! vector opt.
207            zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
208            zhv(ji,jj) = 1._wp / zhv(ji,jj)
209         END DO
210      END DO
211      !
[3959]212      zpsi_uw(:,:,:) = 0._wp
213      zpsi_vw(:,:,:) = 0._wp
214      !
215      DO jk = 2, ikmax                                ! start from 2 : surface value = 0
216         DO jj = 1, jpjm1
217            DO ji = 1, fs_jpim1   ! vector opt.
218               zcuw = 1._wp - ( fsdepw(ji+1,jj,jk) + fsdepw(ji,jj,jk) ) * zhu(ji,jj)
219               zcvw = 1._wp - ( fsdepw(ji,jj+1,jk) + fsdepw(ji,jj,jk) ) * zhv(ji,jj)
220               zcuw = zcuw * zcuw
221               zcvw = zcvw * zcvw
222               zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
223               zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
224               !
225               zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
226               zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
227            END DO
228         END DO
229      END DO
[3995]230      !
231      !                                      !==  transport increased by the MLE induced transport ==!
[3959]232      DO jk = 1, ikmax
233         DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
234            DO ji = 1, fs_jpim1   ! vector opt.
235               pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
236               pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
237            END DO
238         END DO
239         DO jj = 2, jpjm1
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
242                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
243            END DO
244         END DO
245      END DO
246
[3995]247      IF( cdtype == 'TRA') THEN              !==  outputs  ==!
[3959]248         !
249         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
[3995]250         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
[3959]251         !
[4188]252         ! divide by cross distance to give streamfunction with dimensions m^2/s
253         DO jk = 1, ikmax+1
254            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk)/e2u(:,:)
255            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk)/e1v(:,:)
256         END DO
[3959]257         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
258         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
259      ENDIF
260      CALL wrk_dealloc( jpi, jpj, zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH)
261      CALL wrk_dealloc( jpi, jpj, jpk, zpsi_uw, zpsi_vw)
262      CALL wrk_dealloc( jpi, jpj, inml_mle)
263
264      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mle')
265      !
266   END SUBROUTINE tra_adv_mle
267
268
269   SUBROUTINE tra_adv_mle_init
270      !!---------------------------------------------------------------------
271      !!                  ***  ROUTINE tra_adv_mle_init  ***
272      !!
273      !! ** Purpose :   Control the consistency between namelist options for
274      !!              tracer advection schemes and set nadv
275      !!----------------------------------------------------------------------
276      INTEGER  ::   ji, jj, jk   ! dummy loop indices
277      INTEGER  ::   ierr
[4245]278      INTEGER ::    ios                 ! Local integer output status for namelist read
[3959]279      REAL(wp) ::   z1_t2, zfu, zfv                                !    -         -
280      !
281      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
282      !!----------------------------------------------------------------------
283
284
[4245]285      REWIND( numnam_ref )              ! Namelist namtra_adv_mle in reference namelist : Tracer advection scheme
286      READ  ( numnam_ref, namtra_adv_mle, IOSTAT = ios, ERR = 901)
287901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv_mle in reference namelist', lwp )
288
289      REWIND( numnam_cfg )              ! Namelist namtra_adv_mle in configuration namelist : Tracer advection scheme
290      READ  ( numnam_cfg, namtra_adv_mle, IOSTAT = ios, ERR = 902 )
291902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv_mle in configuration namelist', lwp )
[11101]292      IF(lwm .AND. nprint > 2) WRITE ( numond, namtra_adv_mle )
[4245]293
[3959]294      IF(lwp) THEN                     ! Namelist print
295         WRITE(numout,*)
296         WRITE(numout,*) 'tra_adv_mle_init : mixed layer eddy (MLE) advection acting on tracers'
297         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
298         WRITE(numout,*) '   Namelist namtra_adv_mle : mixed layer eddy advection on tracers'
299         WRITE(numout,*) '      use mixed layer eddy (MLE, i.e. Fox-Kemper param) (T/F)      ln_mle    = ', ln_mle
300         WRITE(numout,*) '      MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_mle    = ', nn_mle
301         WRITE(numout,*) '      magnitude of the MLE (typical value: 0.06 to 0.08)           rn_ce     = ', rn_ce
302         WRITE(numout,*) '      scale of ML front (ML radius of deformation) (rn_mle=0)      rn_lf     = ', rn_lf, 'm'
303         WRITE(numout,*) '      maximum time scale of MLE                    (rn_mle=0)      rn_time   = ', rn_time, 's'
304         WRITE(numout,*) '      reference latitude (degrees) of MLE coef.    (rn_mle=1)      rn_lat    = ', rn_lat, 'deg'
305         WRITE(numout,*) '      space interp. of MLD at u-(v-)pts (0=min,1=averaged,2=max)   nn_mld_uv = ', nn_mld_uv
306         WRITE(numout,*) '      =1 no MLE in case of convection ; =0 always MLE              nn_conv   = ', nn_conv
307         WRITE(numout,*) '      Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle
[11101]308         IF(lflush) CALL flush(numout)
[3959]309      ENDIF
310      !
311      IF(lwp) THEN
312         WRITE(numout,*)
313         IF( ln_mle ) THEN
314            WRITE(numout,*) '   Mixed Layer Eddy induced transport added to tracer advection'
315            IF( nn_mle == 0 )   WRITE(numout,*) '   Fox-Kemper et al 2010 formulation'
316            IF( nn_mle == 1 )   WRITE(numout,*) '   New formulation'
317         ELSE
318            WRITE(numout,*) '   Mixed Layer Eddy parametrisation NOT used'
319         ENDIF
[11101]320
321         IF(lflush) CALL flush(numout)
322
[3959]323      ENDIF
324      !
325      IF( ln_mle ) THEN                ! MLE initialisation
326         !
327         rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria
328         IF(lwp) WRITE(numout,*)
329         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
330         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rho_c, 'kg/m3'
[11101]331         IF(lwp .AND. lflush) CALL flush(numout)
[3959]332         !
333         IF( nn_mle == 0 ) THEN           ! MLE array allocation & initialisation
334            ALLOCATE( rfu(jpi,jpj) , rfv(jpi,jpj) , STAT= ierr )
335            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
336            z1_t2 = 1._wp / ( rn_time * rn_time )
337            DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points
338               DO ji = fs_2, jpi   ! vector opt.
339                  zfu = ( ff(ji,jj) + ff(ji,jj-1) ) * 0.5_wp
340                  zfv = ( ff(ji,jj) + ff(ji-1,jj) ) * 0.5_wp
341                  rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 )
342                  rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 )
343               END DO
344            END DO
345            CALL lbc_lnk( rfu, 'U', 1. )   ;   CALL lbc_lnk( rfv, 'V', 1. )
346            !
347         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation
348            rc_f = rn_ce / (  5.e3_wp * 2._wp * omega * SIN( rad * rn_lat )  )
349            !
350         ENDIF
351         !
352         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_mle case)
353         ALLOCATE( r1_ft(jpi,jpj) , STAT= ierr )
354         IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate r1_ft array' )
355         !
356         z1_t2 = 1._wp / ( rn_time * rn_time )
357         r1_ft(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) )
358         r1_ft(:,:) = 1._wp / SQRT(  r1_ft(:,:) * r1_ft(:,:) + z1_t2 )
359         !
360      ENDIF
361      !
362   END SUBROUTINE tra_adv_mle_init
363
364   !!==============================================================================
365END MODULE traadv_mle
Note: See TracBrowser for help on using the repository browser.