Context Navigation

isfcpl.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: 35.8 KB

Line
1	MODULE isfcpl
2	!!======================================================================
3	!! * MODULE isfcpl *
4	!!
5	!! iceshelf coupling module : module managing the coupling between NEMO and an ice sheet model
6	!!
7	!!======================================================================
8	!! History : 4.1 ! 2019-07 (P. Mathiot) Original code
9	!!----------------------------------------------------------------------
10
11	!!----------------------------------------------------------------------
12	!! isfrst : read/write iceshelf variables in/from restart
13	!!----------------------------------------------------------------------
14	USE isf_oce ! ice shelf variable
15	USE isfutils, ONLY : debug
16	USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
17	USE domvvl , ONLY: dom_vvl_zgr ! vertical scale factor interpolation
18	!
19	USE oce ! ocean dynamics and tracers
20	USE in_out_manager ! I/O manager
21	USE iom ! I/O library
22	!
23	IMPLICIT NONE
24
25	PRIVATE
26
27	PUBLIC isfcpl_rst_write, isfcpl_init ! iceshelf restart read and write
28	PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons ! iceshelf correction for ssh, tra, dyn and conservation
29
30	TYPE isfcons
31	INTEGER :: ii ! i global
32	INTEGER :: jj ! j global
33	INTEGER :: kk ! k level
34	REAL(wp):: dvol ! volume increment
35	REAL(wp):: dsal ! salt increment
36	REAL(wp):: dtem ! heat increment
37	REAL(wp):: lon ! lon
38	REAL(wp):: lat ! lat
39	INTEGER :: ngb ! 0/1 (valid location or not (ie on halo or no neigbourg))
40	END TYPE
41	!
42	!! * Substitutions
43	# include "do_loop_substitute.h90"
44	!!----------------------------------------------------------------------
45	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
46	!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
47	!! Software governed by the CeCILL license (see ./LICENSE)
48	!!----------------------------------------------------------------------
49	CONTAINS
50	SUBROUTINE isfcpl_init(Kbb, Kmm, Kaa)
51	!!---------------------------------------------------------------------
52	!! * ROUTINE iscpl_init *
53	!!
54	!! ** Purpose : correct ocean state for new wet cell and horizontal divergence
55	!! correction for the dynamical adjustement
56	!!
57	!! ** Action : - compute ssh on new wet cell
58	!! - compute T/S on new wet cell
59	!! - compute horizontal divergence correction as a volume flux
60	!! - compute the T/S/vol correction increment to keep trend to 0
61	!!
62	!!---------------------------------------------------------------------
63	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
64	!!---------------------------------------------------------------------
65	INTEGER :: id
66	!!----------------------------------------------------------------------
67	!
68	! start on an euler time step
69	neuler = 0
70	!
71	! allocation and initialisation to 0
72	CALL isf_alloc_cpl()
73	!
74	! check presence of variable needed for coupling
75	! iom_varid return 0 if not found
76	id = 1
77	id = id * iom_varid(numror, 'ssmask', ldstop = .false.)
78	id = id * iom_varid(numror, 'tmask' , ldstop = .false.)
79	id = id * iom_varid(numror, 'e3t_n' , ldstop = .false.)
80	id = id * iom_varid(numror, 'e3u_n' , ldstop = .false.)
81	id = id * iom_varid(numror, 'e3v_n' , ldstop = .false.)
82	IF(lwp) WRITE(numout,*) ' isfcpl_init:', id
83	IF (id == 0) THEN
84	IF(lwp) WRITE(numout,*) ' isfcpl_init: restart variables for ice sheet coupling are missing, skip coupling for this leg '
85	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
86	IF(lwp) WRITE(numout,*) ''
87	ELSE
88	! extrapolation ssh
89	CALL isfcpl_ssh(Kbb, Kmm, Kaa)
90	!
91	! extrapolation tracer properties
92	CALL isfcpl_tra(Kmm)
93	!
94	! correction of the horizontal divergence and associated temp. and salt content flux
95	! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
96	! - decide how to manage thickness level change in conservation
97	CALL isfcpl_vol(Kmm)
98	!
99	! apply the 'conservation' method
100	IF ( ln_isfcpl_cons ) CALL isfcpl_cons(Kmm)
101	!
102	END IF
103	!
104	! mask velocity properly (mask used in restart not compatible with new mask)
105	uu(:,:,:,Kmm) = uu(:,:,:,Kmm) * umask(:,:,:)
106	vv(:,:,:,Kmm) = vv(:,:,:,Kmm) * vmask(:,:,:)
107	!
108	! all before fields set to now values
109	ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm)
110	uu (:,:,:,Kbb) = uu (:,:,:,Kmm)
111	vv (:,:,:,Kbb) = vv (:,:,:,Kmm)
112	ssh (:,:,Kbb) = ssh (:,:,Kmm)
113	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
114
115	! prepare writing restart
116	IF( lwxios ) THEN
117	CALL iom_set_rstw_var_active('ssmask')
118	CALL iom_set_rstw_var_active('tmask')
119	CALL iom_set_rstw_var_active('e3t_n')
120	CALL iom_set_rstw_var_active('e3u_n')
121	CALL iom_set_rstw_var_active('e3v_n')
122	END IF
123	!
124	END SUBROUTINE isfcpl_init
125	!
126	SUBROUTINE isfcpl_rst_write(kt, Kmm)
127	!!---------------------------------------------------------------------
128	!! * ROUTINE iscpl_rst_write *
129	!!
130	!! ** Purpose : write icesheet coupling variables in restart
131	!!
132	!!-------------------------- IN --------------------------------------
133	INTEGER, INTENT(in) :: kt
134	INTEGER, INTENT(in) :: Kmm ! ocean time level index
135	!!----------------------------------------------------------------------
136	!
137	IF( lwxios ) CALL iom_swap( cwxios_context )
138	CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask , ldxios = lwxios )
139	CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
140	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n' , e3t(:,:,:,Kmm) , ldxios = lwxios )
141	CALL iom_rstput( kt, nitrst, numrow, 'e3u_n' , e3u(:,:,:,Kmm) , ldxios = lwxios )
142	CALL iom_rstput( kt, nitrst, numrow, 'e3v_n' , e3v(:,:,:,Kmm) , ldxios = lwxios )
143	CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw(:,:,:,Kmm) , ldxios = lwxios )
144	IF( lwxios ) CALL iom_swap( cxios_context )
145	!
146	END SUBROUTINE isfcpl_rst_write
147
148	SUBROUTINE isfcpl_ssh(Kbb, Kmm, Kaa)
149	!!----------------------------------------------------------------------
150	!! * ROUTINE iscpl_ssh *
151	!!
152	!! ** Purpose : basic guess of ssh in new wet cell
153	!!
154	!! ** Method : basic extrapolation from neigbourg cells
155	!!
156	!!----------------------------------------------------------------------
157	!!
158	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
159	!!----------------------------------------------------------------------
160	INTEGER :: ji, jj, jd, jk !! loop index
161	INTEGER :: jip1, jim1, jjp1, jjm1
162	!!
163	REAL(wp):: zsummsk
164	REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
165	!!----------------------------------------------------------------------
166	!
167	CALL iom_get( numror, jpdom_autoglo, 'ssmask' , zssmask_b, ldxios = lrxios ) ! need to extrapolate T/S
168
169	! compute new ssh if we open a full water column
170	! rude average of the closest neigbourgs (e1e2t not taking into account)
171	!
172	zssh(:,:) = ssh(:,:,Kmm)
173	zssmask0(:,:) = zssmask_b(:,:)
174	!
175	DO jd = 1, nn_drown
176	!
177	zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
178	DO jj = 2,jpj-1
179	DO ji = 2,jpi-1
180	jip1=ji+1; jim1=ji-1;
181	jjp1=jj+1; jjm1=jj-1;
182	!
183	zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
184	!
185	IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
186	ssh(ji,jj,Kmm)=( zssh(jip1,jj)*zssmask0(jip1,jj) &
187	& + zssh(jim1,jj)*zssmask0(jim1,jj) &
188	& + zssh(ji,jjp1)*zssmask0(ji,jjp1) &
189	& + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
190	zssmask_b(ji,jj) = 1._wp
191	ENDIF
192	END DO
193	END DO
194	!
195	zssh(:,:) = ssh(:,:,Kmm)
196	zssmask0(:,:) = zssmask_b(:,:)
197	!
198	CALL lbc_lnk_multi( 'iscplrst', zssh, 'T', 1., zssmask0, 'T', 1. )
199	!
200	END DO
201	!
202	! update ssh(:,:,Kmm)
203	ssh(:,:,Kmm) = zssh(:,:) * ssmask(:,:)
204	!
205	ssh(:,:,Kbb) = ssh(:,:,Kmm)
206	!
207	IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',ssh(:,:,Kmm))
208	!
209	! recompute the vertical scale factor, depth and water thickness
210	IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
211	IF(lwp) write(numout,*) '~~~~~~~~~~~'
212	DO jk = 1, jpk
213	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
214	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
215	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
216	END DO
217	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
218	CALL dom_vvl_zgr(Kbb, Kmm, Kaa)
219	!
220	END SUBROUTINE isfcpl_ssh
221
222	SUBROUTINE isfcpl_tra(Kmm)
223	!!----------------------------------------------------------------------
224	!! * ROUTINE iscpl_tra *
225	!!
226	!! ** Purpose : compute new tn, sn in case of evolving geometry of ice shelves
227	!!
228	!! ** Method : tn, sn : basic extrapolation from neigbourg cells
229	!!
230	!!----------------------------------------------------------------------
231	INTEGER, INTENT(in) :: Kmm ! ocean time level index
232	!!----------------------------------------------------------------------
233	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
234	!REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b !! depth w before
235	!!
236	INTEGER :: ji, jj, jk, jd !! loop index
237	INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
238	!!
239	REAL(wp):: zsummsk
240	REAL(wp):: zdz, zdzm1, zdzp1
241	!!
242	REAL(wp), DIMENSION(jpi,jpj) :: zdmask
243	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn
244	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d
245	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
246	!!----------------------------------------------------------------------
247	!
248	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b, ldxios = lrxios ) ! need to extrapolate T/S
249	!CALL iom_get( numror, jpdom_autoglo, 'wmask' , zwmask_b, ldxios = lrxios ) ! need to extrapolate T/S
250	!CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
251	!
252	!
253	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
254	!PM: Is this IF needed since change to VVL by default
255	!bugged : to be corrected (PM)
256	! back up original t/s/mask
257	!tsb (:,:,:,:) = ts(:,:,:,:,Kmm)
258	!
259	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
260
261	! IF (.NOT.ln_linssh) THEN
262	! DO jk = 2,jpk-1
263	! DO jj = 1,jpj
264	! DO ji = 1,jpi
265	! IF (wmask(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ztmask_b(ji,jj,1)==0._wp) ) THEN
266	!
267	! !compute weight
268	! zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw(ji,jj,jk+1,Kmm))
269	! zdzm1 = MAX(0._wp,gdepw(ji,jj,jk ,Kmm) - pdepw_b(ji,jj,jk ))
270	! zdz = e3t(ji,jj,jk,Kmm) - zdzp1 - zdzm1 ! if isf : e3t = gdepw(ji,jj,jk+1,Kmm)- gdepw(ji,jj,jk,Kmm)
271	!
272	! IF (zdz .LT. 0._wp) THEN
273	! CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
274	! END IF
275	!
276	! ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb) &
277	! & + zdz *ts(ji,jj,jk ,jp_tem,Kbb) &
278	! & + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb) )/e3t(ji,jj,jk,Kmm)
279	!
280	! ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb) &
281	! & + zdz *ts(ji,jj,jk ,jp_sal,Kbb) &
282	! & + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb) )/e3t(ji,jj,jk,Kmm)
283	!
284	! END IF
285	! END DO
286	! END DO
287	! END DO
288	! END IF
289
290	zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
291	ztmask0(:,:,:) = ztmask_b(:,:,:)
292	ztmask1(:,:,:) = ztmask_b(:,:,:)
293	!
294	! iterate the extrapolation processes nn_drown times
295	DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
296	DO jk = 1,jpk-1
297	!
298	! define new wet cell
299	zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
300	!
301	DO jj = 2,jpj-1
302	DO ji = 2,jpi-1
303	jip1=ji+1; jim1=ji-1;
304	jjp1=jj+1; jjm1=jj-1;
305	!
306	! check if a wet neigbourg cell is present
307	zsummsk = ztmask0(jip1,jj ,jk) + ztmask0(jim1,jj ,jk) &
308	+ ztmask0(ji ,jjp1,jk) + ztmask0(ji ,jjm1,jk)
309	!
310	! if neigbourg wet cell available at the same level
311	IF ( zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
312	!
313	! horizontal basic extrapolation
314	ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj ,jk,1) * ztmask0(jip1,jj ,jk) &
315	& + zts0(jim1,jj ,jk,1) * ztmask0(jim1,jj ,jk) &
316	& + zts0(ji ,jjp1,jk,1) * ztmask0(ji ,jjp1,jk) &
317	& + zts0(ji ,jjm1,jk,1) * ztmask0(ji ,jjm1,jk) ) / zsummsk
318	ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj ,jk,2) * ztmask0(jip1,jj ,jk) &
319	& + zts0(jim1,jj ,jk,2) * ztmask0(jim1,jj ,jk) &
320	& + zts0(ji ,jjp1,jk,2) * ztmask0(ji ,jjp1,jk) &
321	& + zts0(ji ,jjm1,jk,2) * ztmask0(ji ,jjm1,jk) ) / zsummsk
322	!
323	! update mask for next pass
324	ztmask1(ji,jj,jk)=1
325	!
326	! in case no neigbourg wet cell available at the same level
327	! check if a wet cell is available below
328	ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
329	!
330	! vertical extrapolation if horizontal extrapolation failed
331	jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
332	!
333	! check if a wet neigbourg cell is present
334	zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
335	IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
336	ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) &
337	& + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
338	ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) &
339	& + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
340	!
341	! update mask for next pass
342	ztmask1(ji,jj,jk)=1._wp
343	END IF
344	END IF
345	END DO
346	END DO
347	END DO
348	!
349	! update temperature and salinity and mask
350	zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
351	ztmask0(:,:,:) = ztmask1(:,:,:)
352	!
353	CALL lbc_lnk_multi( 'iscplrst', zts0(:,:,:,jp_tem), 'T', 1., zts0(:,:,:,jp_sal), 'T', 1., ztmask0, 'T', 1.)
354	!
355	END DO ! nn_drown
356	!
357	! mask new ts(:,:,:,:,Kmm) field
358	ts(:,:,:,jp_tem,Kmm) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
359	ts(:,:,:,jp_sal,Kmm) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
360	!
361	! sanity check
362	! -----------------------------------------------------------------------------------------
363	! case we open a cell but no neigbour cells available to get an estimate of T and S
364	DO_3D_11_11( 1,jpk-1 )
365	IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp) &
366	& CALL ctl_stop('STOP', 'failing to fill all new weet cell, &
367	& try increase nn_drown or activate XXXX &
368	& in your domain cfg computation' )
369	END_3D
370	!
371	END SUBROUTINE isfcpl_tra
372
373	SUBROUTINE isfcpl_vol(Kmm)
374	!!----------------------------------------------------------------------
375	!! * ROUTINE iscpl_vol *
376	!!
377	!! ** Purpose : compute the correction of the local divergence to apply
378	!! during the first time step after the coupling.
379	!!
380	!! ** Method : - compute horizontal vol div. before/after coupling
381	!! - compute vertical input
382	!! - compute correction
383	!!
384	!!----------------------------------------------------------------------
385	!!
386	INTEGER, INTENT(in) :: Kmm ! ocean time level index
387	!!----------------------------------------------------------------------
388	INTEGER :: ji, jj, jk
389	INTEGER :: ikb, ikt
390	!!
391	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln ! vol flux div. before/after coupling
392	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b ! vertical scale factor before/after coupling
393	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before coupling
394	!!----------------------------------------------------------------------
395	!
396	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b, ldxios = lrxios )
397	CALL iom_get( numror, jpdom_autoglo, 'e3u_n' , ze3u_b , ldxios = lrxios )
398	CALL iom_get( numror, jpdom_autoglo, 'e3v_n' , ze3v_b , ldxios = lrxios )
399	!
400	! 1.0: compute horizontal volume flux divergence difference before-after coupling
401	!
402	DO jk = 1, jpk ! Horizontal slab
403	! 1.1: get volume flux before coupling (>0 out)
404	DO_2D_00_00
405	zqvolb(ji,jj,jk) = ( e2u(ji,jj) * ze3u_b(ji,jj,jk) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj ) * ze3u_b(ji-1,jj ,jk) * uu(ji-1,jj ,jk,Kmm) &
406	& + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vv(ji,jj,jk,Kmm) - e1v(ji ,jj-1) * ze3v_b(ji ,jj-1,jk) * vv(ji ,jj-1,jk,Kmm) ) &
407	& * ztmask_b(ji,jj,jk)
408	END_2D
409	!
410	! 1.2: get volume flux after coupling (>0 out)
411	! properly mask velocity
412	! (velocity are still mask with old mask at this stage)
413	uu(:,:,jk,Kmm) = uu(:,:,jk,Kmm) * umask(:,:,jk)
414	vv(:,:,jk,Kmm) = vv(:,:,jk,Kmm) * vmask(:,:,jk)
415	! compute volume flux divergence after coupling
416	DO_2D_00_00
417	zqvoln(ji,jj,jk) = ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) - e2u(ji-1,jj ) * e3u(ji-1,jj ,jk,Kmm) * uu(ji-1,jj ,jk,Kmm) &
418	& + e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) - e1v(ji ,jj-1) * e3v(ji ,jj-1,jk,Kmm) * vv(ji ,jj-1,jk,Kmm) ) &
419	& * tmask(ji,jj,jk)
420	END_2D
421	!
422	! 1.3: get 3d volume flux difference (before - after cpl) (>0 out)
423	! correction to add is _b - _n
424	risfcpl_vol(:,:,jk) = zqvolb(:,:,jk) - zqvoln(:,:,jk)
425	END DO
426	!
427	! 2.0: include the contribution of the vertical velocity in the volume flux correction
428	!
429	DO_2D_00_00
430	!
431	ikt = mikt(ji,jj)
432	IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN
433	risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1)) ! test sign
434	ENDIF
435	!
436	END_2D
437	!
438	CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
439	!
440	! 3.0: set total correction (div, tr(:,:,:,:,Krhs), ssh)
441	!
442	! 3.1: mask volume flux divergence correction
443	risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
444	!
445	! 3.2: get 3d tr(:,:,:,:,Krhs) increment to apply at the first time step
446	! temperature and salt content flux computed using local ts(:,:,:,:,Kmm)
447	! (very simple advection scheme)
448	! (>0 out)
449	risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_tem,Kmm)
450	risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_sal,Kmm)
451	!
452	! 3.3: ssh correction (for dynspg_ts)
453	risfcpl_ssh(:,:) = 0.0
454	DO jk = 1,jpk
455	risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
456	END DO
457
458	END SUBROUTINE isfcpl_vol
459
460	SUBROUTINE isfcpl_cons(Kmm)
461	!!----------------------------------------------------------------------
462	!! * ROUTINE iscpl_cons *
463	!!
464	!! ** Purpose : compute the corrective increment in volume/salt/heat to put back the vol/heat/salt
465	!! removed or added during the coupling processes (wet or dry new cell)
466	!!
467	!! ** Method : - compare volume/heat/salt before and after
468	!! - look for the closest wet cells (share amoung neigbourgs if there are)
469	!! - build the correction increment to applied at each time step
470	!!
471	!!----------------------------------------------------------------------
472	!
473	TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
474	!
475	!!----------------------------------------------------------------------
476	INTEGER, INTENT(in) :: Kmm ! ocean time level index
477	!!----------------------------------------------------------------------
478	INTEGER :: ji , jj , jk , jproc ! loop index
479	INTEGER :: jip1 , jim1, jjp1, jjm1 ! dummy indices
480	INTEGER :: iig , ijg, ik ! dummy indices
481	INTEGER :: jisf ! start, end and current position in the increment array
482	INTEGER :: ingb, ifind ! 0/1 target found or need to be found
483	INTEGER :: nisfl_area ! global number of cell concerned by the wet->dry case
484	INTEGER, DIMENSION(jpnij) :: nisfl ! local number of cell concerned by the wet->dry case
485	!
486	REAL(wp) :: z1_sum, z1_rdtiscpl
487	REAL(wp) :: zdtem, zdsal, zdvol, zratio ! tem, sal, vol increment
488	REAL(wp) :: zlon , zlat ! target location
489	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before
490	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b ! scale factor before
491	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b ! scale factor before
492	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b ! scale factor before
493	!!----------------------------------------------------------------------
494
495	!==============================================================================
496	! 1.0: initialisation
497	!==============================================================================
498
499	! get restart variable
500	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b(:,:,:), ldxios = lrxios ) ! need to extrapolate T/S
501	CALL iom_get( numror, jpdom_autoglo, 'e3t_n' , ze3t_b(:,:,:) , ldxios = lrxios )
502	CALL iom_get( numror, jpdom_autoglo, 'tn' , zt_b(:,:,:) , ldxios = lrxios )
503	CALL iom_get( numror, jpdom_autoglo, 'sn' , zs_b(:,:,:) , ldxios = lrxios )
504
505	! compute run length
506	nstp_iscpl = nitend - nit000 + 1
507	rdt_iscpl = nstp_iscpl * rn_rdt
508	z1_rdtiscpl = 1._wp / rdt_iscpl
509
510	IF (lwp) WRITE(numout,*) ' nb of stp for cons = ', nstp_iscpl
511	IF (lwp) WRITE(numout,*) ' coupling time step = ', rdt_iscpl
512
513	! initialisation correction
514	risfcpl_cons_vol = 0.0
515	risfcpl_cons_ssh = 0.0
516	risfcpl_cons_tsc = 0.0
517
518	!==============================================================================
519	! 2.0: diagnose the heat, salt and volume input and compute the correction variable
520	! for case where we wet a cell or cell still wet (no change in cell status)
521	!==============================================================================
522
523	DO jk = 1,jpk-1
524	DO jj = nldj,nlej
525	DO ji = nldi,nlei
526
527	! volume diff
528	zdvol = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) - ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
529
530	! heat diff
531	zdtem = ts (ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
532	- zt_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
533
534	! salt diff
535	zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
536	- zs_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
537
538	! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
539	! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
540	risfcpl_cons_vol(ji,jj,jk) = ( zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk) ) * z1_rdtiscpl
541	risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
542	risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
543
544	END DO
545	END DO
546	END DO
547	!
548	!==============================================================================
549	! 3.0: diagnose the heat, salt and volume input and compute the correction variable
550	! for case where we close a cell
551	!==============================================================================
552	!
553	! compute the total number of point receiving a correction increment for each processor
554	! local
555	nisfl(:)=0
556	DO jk = 1,jpk-1
557	DO jj = nldj,nlej
558	DO ji = nldi,nlei
559	jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
560	IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
561	ENDDO
562	ENDDO
563	ENDDO
564	!
565	! global
566	CALL mpp_sum('isfcpl',nisfl )
567	!
568	! allocate list of point receiving correction
569	ALLOCATE(zisfpts(nisfl(narea)))
570	!
571	zisfpts(:) = isfcons(0,0,0,-HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), 0)
572	!
573	! start computing the correction and fill zisfpts
574	! local
575	jisf = 0
576	DO jk = 1,jpk-1
577	DO jj = nldj,nlej
578	DO ji = nldi,nlei
579	IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
580
581	jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
582
583	zdvol = risfcpl_cons_vol(ji,jj,jk )
584	zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
585	zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
586
587	IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
588	! spread correction amoung neigbourg wet cells (horizontal direction first)
589	! as it is a rude correction corner and lateral cell have the same weight
590	!
591	z1_sum = 1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
592	!
593	! lateral cells
594	IF (tmask(jip1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
595	IF (tmask(jim1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
596	IF (tmask(ji ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
597	IF (tmask(ji ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
598	!
599	! corner cells
600	IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
601	IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
602	IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
603	IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
604	!
605	ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
606	! spread correction amoung neigbourg wet cells (vertical direction)
607	CALL update_isfpts(zisfpts, jisf, ji , jj , jk+1, zdvol, zdsal, zdtem, 1., 0)
608	ELSE
609	! need to find where to put correction in later on
610	CALL update_isfpts(zisfpts, jisf, ji , jj , jk , zdvol, zdsal, zdtem, 1., 1)
611	END IF
612	END IF
613	END DO
614	END DO
615	END DO
616	!
617	! share data among all processes because for some point we need to find the closest wet point (could be on other process)
618	DO jproc=1,jpnij
619	!
620	! share total number of isf point treated for proc jproc
621	IF (jproc==narea) THEN
622	nisfl_area=nisfl(jproc)
623	ELSE
624	nisfl_area=0
625	END IF
626	CALL mpp_max('isfcpl',nisfl_area)
627	!
628	DO jisf = 1,nisfl_area
629	!
630	IF (jproc==narea) THEN
631	! indices (conversion to global indices and sharing)
632	iig = zisfpts(jisf)%ii ; ijg = zisfpts(jisf)%jj ; ik = zisfpts(jisf)%kk
633	!
634	! data
635	zdvol = zisfpts(jisf)%dvol ; zdsal = zisfpts(jisf)%dsal ; zdtem = zisfpts(jisf)%dtem
636	!
637	! location
638	zlat = zisfpts(jisf)%lat ; zlon = zisfpts(jisf)%lon
639	!
640	! find flag
641	ingb = zisfpts(jisf)%ngb
642	ELSE
643	iig =0 ; ijg =0 ; ik =0
644	zdvol=-HUGE(1.0) ; zdsal=-HUGE(1.0) ; zdtem=-HUGE(1.0)
645	zlat =-HUGE(1.0) ; zlon =-HUGE(1.0)
646	ingb = 0
647	END IF
648	!
649	! share data (need synchronisation of data as get_correction call a global com)
650	CALL mpp_max('isfcpl',iig) ; CALL mpp_max('isfcpl',ijg) ; CALL mpp_max('isfcpl',ik)
651	CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
652	CALL mpp_max('isfcpl',zlat) ; CALL mpp_max('isfcpl',zlon)
653	CALL mpp_max('isfcpl',ingb)
654	!
655	! fill the 3d correction array
656	CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
657	END DO
658	END DO
659	!
660	!==============================================================================
661	! 4.0: finalisation and compute ssh equivalent of the volume correction
662	!==============================================================================
663	!
664	! mask (>0 out)
665	risfcpl_cons_vol(:,:,: ) = risfcpl_cons_vol(:,:,: ) * tmask(:,:,:)
666	risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
667	risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
668	!
669	! add lbclnk
670	CALL lbc_lnk_multi( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1., risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1., &
671	& risfcpl_cons_vol(:,:,:) , 'T', 1.)
672	!
673	! ssh correction (for dynspg_ts)
674	DO jk = 1,jpk
675	risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
676	END DO
677	risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
678	!
679	END SUBROUTINE isfcpl_cons
680	!
681	SUBROUTINE update_isfpts(sisfpts, kpts, ki, kj, kk, pdvol, pdsal, pdtem, pratio, kfind)
682	!!---------------------------------------------------------------------
683	!! * ROUTINE update_isfpts *
684	!!
685	!! ** Purpose : if a cell become dry, we need to put the corrective increment elsewhere
686	!!
687	!! ** Action : update the list of point
688	!!
689	!!----------------------------------------------------------------------
690	!!----------------------------------------------------------------------
691	TYPE(isfcons), DIMENSION(:), INTENT(inout) :: sisfpts
692	INTEGER, INTENT(inout) :: kpts
693	!!----------------------------------------------------------------------
694	INTEGER, INTENT(in ) :: ki, kj, kk ! target location (kfind=0)
695	! ! or source location (kfind=1)
696	INTEGER, INTENT(in ), OPTIONAL :: kfind ! 0 target cell already found
697	! ! 1 target to be determined
698	REAL(wp), INTENT(in ) :: pdvol, pdsal, pdtem, pratio ! vol/sal/tem increment
699	! ! and ratio in case increment span over multiple cells.
700	!!----------------------------------------------------------------------
701	INTEGER :: ifind
702	!!----------------------------------------------------------------------
703	!
704	! increment position
705	kpts = kpts + 1
706	!
707	! define if we need to look for closest valid wet cell (no neighbours or neigbourg on halo)
708	IF ( PRESENT(kfind) ) THEN
709	ifind = kfind
710	ELSE
711	ifind = ( 1 - tmask_h(ki,kj) ) * tmask(ki,kj,kk)
712	END IF
713	!
714	! update isfpts structure
715	sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
716	!
717	END SUBROUTINE update_isfpts
718	!
719	SUBROUTINE get_correction( ki, kj, kk, plon, plat, pvolinc, psalinc, pteminc, kfind)
720	!!---------------------------------------------------------------------
721	!! * ROUTINE get_correction *
722	!!
723	!! ** Action : - Find the closest valid cell if needed (wet and not on the halo)
724	!! - Scale the correction depending of pratio (case where multiple wet neigbourgs)
725	!! - Fill the correction array
726	!!
727	!!----------------------------------------------------------------------
728	INTEGER , INTENT(in) :: ki, kj, kk, kfind ! target point indices
729	REAL(wp), INTENT(in) :: plon, plat ! target point lon/lat
730	REAL(wp), INTENT(in) :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
731	!!----------------------------------------------------------------------
732	INTEGER :: jj, ji, iig, ijg
733	!!----------------------------------------------------------------------
734	!
735	! define global indice of correction location
736	iig = ki ; ijg = kj
737	IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
738	!
739	! fill the correction array
740	DO jj = mj0(ijg),mj1(ijg)
741	DO ji = mi0(iig),mi1(iig)
742	! correct the vol_flx and corresponding heat/salt flx in the closest cell
743	risfcpl_cons_vol(ji,jj,kk) = risfcpl_cons_vol(ji,jj,kk ) + pvolinc
744	risfcpl_cons_tsc(ji,jj,kk,jp_sal) = risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
745	risfcpl_cons_tsc(ji,jj,kk,jp_tem) = risfcpl_cons_tsc(ji,jj,kk,jp_tem) + pteminc
746	END DO
747	END DO
748
749	END SUBROUTINE get_correction
750
751	END MODULE isfcpl

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

Download in other formats: