Context Navigation

isftbl.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.

File size: 12.3 KB

Line
1	MODULE isftbl
2	!!======================================================================
3	!! * MODULE isftbl *
4	!! isftbl module : compute properties of top boundary layer
5	!!======================================================================
6	!! History : 4.1 ! 2019-09 (P. Mathiot) original code
7	!!----------------------------------------------------------------------
8
9	!!----------------------------------------------------------------------
10	!! isftbl : routine to compute :
11	!! - geometry of the ice shelf tbl (isf_tbl_lvl, isftbl_ktop, isftbl_kbot)
12	!! (top and bottom level, thickness and fraction of deepest level affected)
13	!! - tbl averaged properties (isf_tbl, isf_tbl_avg)
14	!!----------------------------------------------------------------------
15
16	USE isf_oce ! ice shelf variables
17
18	USE dom_oce ! vertical scale factor and depth
19
20	IMPLICIT NONE
21
22	PRIVATE
23
24	PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
25	!! * Substitutions
26	# include "do_loop_substitute.h90"
27
28	CONTAINS
29
30	SUBROUTINE isf_tbl( Kmm, pvarin, pvarout, cd_ptin, ktop, phtbl, kbot, pfrac )
31	!!--------------------------------------------------------------------
32	!! * SUBROUTINE isf_tbl *
33	!!
34	!! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
35	!!
36	!! ** Method : Average properties over a specific thickness
37	!!
38	!! ** Reference : inspired from : Losch, Modeling ice shelf cavities in a z coordinate ocean general circulation model
39	!! https://doi.org/10.1029/2007JC004368 , 2008
40	!!
41	!!-------------------------- OUT -------------------------------------
42	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pvarout ! 2d average of pvarin
43	!!-------------------------- IN -------------------------------------
44	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
45	CHARACTER(len=1) , INTENT(in ) :: cd_ptin ! point of variable in/out
46	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pvarin ! 3d variable to average over the tbl
47	INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop ! top level
48	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl ! tbl thickness
49	!!-------------------------- IN OPTIONAL -----------------------------
50	INTEGER, DIMENSION(jpi,jpj), OPTIONAL, INTENT(in ) :: kbot ! bottom level
51	REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in ) :: pfrac ! fraction of bottom cell affected by tbl
52	!!--------------------------------------------------------------------
53	INTEGER :: ji, jj ! loop index
54	INTEGER , DIMENSION(jpi,jpj) :: ikbot ! bottom level of the tbl
55	REAL(wp), DIMENSION(jpi,jpj) :: zvarout ! 2d average of pvarin
56	REAL(wp), DIMENSION(jpi,jpj) :: zhtbl ! thickness of the tbl
57	REAL(wp), DIMENSION(jpi,jpj) :: zfrac ! thickness of the tbl
58	!!--------------------------------------------------------------------
59	!
60	SELECT CASE ( cd_ptin )
61	CASE ( 'U' )
62	!
63	! copy phtbl (phtbl is INTENT in as we don't want to change it)
64	zhtbl = phtbl
65	!
66	! compute tbl lvl and thickness
67	CALL isf_tbl_lvl( hu(:,:,Kmm), e3u(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
68	!
69	! compute tbl property at U point
70	CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u(:,:,:,Kmm), pvarin, zvarout )
71	!
72	! compute tbl property at T point
73	pvarout(1,:) = 0._wp
74	DO_2D_11_01
75	pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
76	END_2D
77	! lbclnk not needed as a final communication is done after the computation of fwf
78	!
79	CASE ( 'V' )
80	!
81	! copy phtbl (phtbl is INTENT in as we don't want to change it)
82	zhtbl = phtbl
83	!
84	! compute tbl lvl and thickness
85	CALL isf_tbl_lvl( hv(:,:,Kmm), e3v(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
86	!
87	! compute tbl property at V point
88	CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v(:,:,:,Kmm), pvarin, zvarout )
89	!
90	! pvarout is an averaging of wet point
91	pvarout(:,1) = 0._wp
92	DO_2D_01_11
93	pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
94	END_2D
95	! lbclnk not needed as a final communication is done after the computation of fwf
96	!
97	CASE ( 'T' )
98	!
99	! compute tbl property at T point
100	CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, e3t(:,:,:,Kmm), pvarin, pvarout )
101	!
102	END SELECT
103	!
104	END SUBROUTINE isf_tbl
105
106	SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
107	!!--------------------------------------------------------------------
108	!! * ROUTINE isf_tbl_avg *
109	!!
110	!! ** Purpose : compute mean property in the boundary layer
111	!!
112	!! ** Method : Depth average is made between the top level ktop and the bottom level kbot
113	!! over a thickness phtbl. The bottom level is partially counted (pfrac).
114	!!
115	!!-------------------------- OUT -------------------------------------
116	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pvarout ! tbl property averaged over phtbl between level ktop and kbot
117	!!-------------------------- IN -------------------------------------
118	INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop, kbot ! top and bottom level of the top boundary layer
119	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
120	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3 ! vertical scale factor
121	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pvarin ! tbl property to average between ktop, kbot over phtbl
122	!!--------------------------------------------------------------------
123	INTEGER :: ji,jj,jk ! loop indices
124	INTEGER :: ikt, ikb ! top and bottom levels
125	!!--------------------------------------------------------------------
126	!
127	! compute tbl top.bottom level and thickness
128	DO_2D_11_11
129	!
130	! tbl top/bottom indices initialisation
131	ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
132	!
133	! level fully include in the ice shelf boundary layer
134	pvarout(ji,jj) = SUM( pvarin(ji,jj,ikt:ikb-1) * pe3(ji,jj,ikt:ikb-1) ) / phtbl(ji,jj)
135	!
136	! level partially include in ice shelf boundary layer
137	pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * pe3(ji,jj,ikb) / phtbl(ji,jj) * pfrac(ji,jj)
138	!
139	END_2D
140
141	END SUBROUTINE isf_tbl_avg
142
143	SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
144	!!--------------------------------------------------------------------
145	!! * ROUTINE isf_tbl_lvl *
146	!!
147	!! ** Purpose : - compute bottom level off the top boundary layer
148	!! - thickness of the top boundary layer
149	!! - fraction of the bottom level affected by the tbl
150	!!
151	!!-------------------------- OUT --------------------------------------
152	INTEGER, DIMENSION(jpi,jpj) , INTENT( out) :: kbot ! bottom level of the top boundary layer
153	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pfrac ! fraction of bottom level in the tbl
154	!!-------------------------- IN --------------------------------------
155	INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop ! top level of the top boundary layer
156	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phw ! water column thickness
157	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3 ! vertical scale factor
158	!!-------------------------- INOUT ------------------------------------
159	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: phtbl ! top boundary layer thickness
160	!!---------------------------------------------------------------------
161	INTEGER :: ji,jj,jk
162	INTEGER :: ikt, ikb
163	!!---------------------------------------------------------------------
164	!
165	! get htbl
166	DO_2D_11_11
167	!
168	! tbl top/bottom indices initialisation
169	ikt = ktop(ji,jj)
170	!
171	! limit the tbl to water thickness.
172	phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) )
173	!
174	! thickness of boundary layer must be at least the top level thickness
175	phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) )
176	!
177	END_2D
178	!
179	! get ktbl
180	CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
181	!
182	! get pfrac
183	DO_2D_11_11
184	!
185	! tbl top/bottom indices initialisation
186	ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
187	!
188	! proportion of the bottom cell included in ice shelf boundary layer
189	pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb)
190	!
191	END_2D
192	!
193	END SUBROUTINE isf_tbl_lvl
194	!
195	SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
196	!!--------------------------------------------------------------------
197	!! * ROUTINE isf_tbl_bot *
198	!!
199	!! ** Purpose : compute bottom level of the isf top boundary layer
200	!!
201	!!-------------------------- OUT -------------------------------------
202	INTEGER, DIMENSION(jpi,jpj) , INTENT( out) :: kbot ! bottom level of the top boundary layer
203	!!-------------------------- IN -------------------------------------
204	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl ! top boundary layer thickness
205	INTEGER, DIMENSION(jpi,jpj) , INTENT(in ) :: ktop ! top level of the top boundary layer
206	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3 ! vertical scale factor
207	!!--------------------------------------------------------------------
208	INTEGER :: ji, jj
209	INTEGER :: ikt, ikb
210	!!--------------------------------------------------------------------
211	!
212	! phtbl need to be bounded by water column thickness before
213	! test: if htbl = water column thickness, should return mbathy
214	! test: if htbl = 0 should return ktop (phtbl cap to e3t(ji,jj,1))
215	!
216	! get ktbl
217	DO_2D_11_11
218	!
219	! determine the deepest level influenced by the boundary layer
220	ikt = ktop(ji,jj)
221	ikb = ikt
222	DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ; ikb = ikb + 1 ; END DO
223	kbot(ji,jj) = ikb - 1
224	!
225	END_2D
226	!
227	END SUBROUTINE isf_tbl_kbot
228	!
229	SUBROUTINE isf_tbl_ktop(pdep, ktop)
230	!!--------------------------------------------------------------------
231	!! * ROUTINE isf_tbl_top *
232	!!
233	!! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
234	!!
235	!!-------------------------- OUT -------------------------------------
236	INTEGER, DIMENSION(jpi,jpj), INTENT( out) :: ktop ! top level affected by the ice shelf parametrisation
237	!!-------------------------- IN -------------------------------------
238	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdep ! top depth of the parametrisation influence
239	!!--------------------------------------------------------------------
240	INTEGER :: ji,jj
241	INTEGER :: ikt
242	!!--------------------------------------------------------------------
243	!
244	! if we need to recompute the top level at every time stepcompute top level (z*, z~)
245	! in case of weak ht variation we can assume the top level of htbl to be constant
246	! => only done using gdepw_0
247	! be sure pdep is already correctly bounded
248	! test: this routine run on isfdraft should return mikt
249	! test: this routine run with pdep = 0 should return 1
250	!
251	DO_2D_11_11
252	! comput ktop
253	ikt = 2
254	DO WHILE ( gdepw_0(ji,jj,ikt) <= pdep(ji,jj ) ) ; ikt = ikt + 1 ; END DO
255	ktop(ji,jj) = ikt - 1
256	!
257	! update pdep
258	pdep(ji,jj) = gdepw_0(ji,jj,ktop(ji,jj))
259	END_2D
260	!
261	END SUBROUTINE isf_tbl_ktop
262
263	END MODULE isftbl

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/ISF/isftbl.F90 @ 12340

Download in other formats: