Context Navigation

tramle.F90 @ 12340

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

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 16.7 KB

Line
1	MODULE 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	# include "do_loop_substitute.h90"
52	!!----------------------------------------------------------------------
53	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
54	!! $Id$
55	!! Software governed by the CeCILL license (see ./LICENSE)
56	!!----------------------------------------------------------------------
57	CONTAINS
58
59	SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype, Kmm )
60	!!----------------------------------------------------------------------
61	!! * ROUTINE tra_mle_trp *
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	!!
74	!! ** Action : - (pu,pv,pw) increased by the mle transport
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
83	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
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	!
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
96	!!----------------------------------------------------------------------
97	!
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)
101	IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m
102	DO_3DS_11_11( jpkm1, nlb10, -1 )
103	IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer
104	END_3D
105	ENDIF
106	ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation
107	!
108	!
109	zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==!
110	zbm (:,:) = 0._wp
111	zn2 (:,:) = 0._wp
112	DO_3D_11_11( 1, ikmax )
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
115	zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
116	zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
117	END_3D
118
119	SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts
120	CASE ( 0 ) != min of the 2 neighbour MLDs
121	DO_2D_10_10
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
125	CASE ( 1 ) != average of the 2 neighbour MLDs
126	DO_2D_10_10
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
130	CASE ( 2 ) != max of the 2 neighbour MLDs
131	DO_2D_10_10
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
135	END SELECT
136	! ! convert density into buoyancy
137	zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )
138	!
139	!
140	! !== Magnitude of the MLE stream function ==!
141	!
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
148	DO_2D_10_10
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
157	!
158	ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
159	DO_2D_10_10
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
166	ENDIF
167	!
168	IF( nn_conv == 1 ) THEN ! No MLE in case of convection
169	DO_2D_10_10
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
173	ENDIF
174	!
175	! !== structure function value at uw- and vw-points ==!
176	DO_2D_10_10
177	zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu
178	zhv(ji,jj) = 1._wp / zhv(ji,jj)
179	END_2D
180	!
181	zpsi_uw(:,:,:) = 0._wp
182	zpsi_vw(:,:,:) = 0._wp
183	!
184	DO_3D_10_10( 2, ikmax )
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
195	!
196	! !== transport increased by the MLE induced transport ==!
197	DO jk = 1, ikmax
198	DO_2D_10_10
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
202	DO_2D_00_00
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
206	END DO
207
208	IF( cdtype == 'TRA') THEN !== outputs ==!
209	!
210	zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f
211	CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f
212	!
213	! divide by cross distance to give streamfunction with dimensions m^2/s
214	DO jk = 1, ikmax+1
215	zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
216	zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
217	END DO
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	!
222	END SUBROUTINE tra_mle_trp
223
224
225	SUBROUTINE tra_mle_init
226	!!---------------------------------------------------------------------
227	!! * ROUTINE tra_mle_init *
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
234	INTEGER :: ios ! Local integer output status for namelist read
235	REAL(wp) :: z1_t2, zfu, zfv ! - -
236	!
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
238	!!----------------------------------------------------------------------
239
240	READ ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901)
241	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_mle in reference namelist' )
242
243	READ ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 )
244	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' )
245	IF(lwm) WRITE ( numond, namtra_mle )
246
247	IF(lwp) THEN ! Namelist print
248	WRITE(numout,*)
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'
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
261	ENDIF
262	!
263	IF(lwp) THEN
264	WRITE(numout,*)
265	IF( ln_mle ) THEN
266	WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to tracer advection'
267	IF( nn_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
268	IF( nn_mle == 1 ) WRITE(numout,*) ' New formulation'
269	ELSE
270	WRITE(numout,*) ' ==>>> Mixed Layer Eddy parametrisation NOT used'
271	ENDIF
272	ENDIF
273	!
274	IF( ln_mle ) THEN ! MLE initialisation
275	!
276	rb_c = grav * rn_rho_c_mle /rau0 ! Mixed Layer buoyancy criteria
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 )
285	DO_2D_01_01
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
291	CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. )
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 )
303	r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
304	!
305	ENDIF
306	!
307	END SUBROUTINE tra_mle_init
308
309	!!==============================================================================
310	END MODULE tramle

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/tramle.F90 @ 12340

Download in other formats: