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.

isfcpl.F90 in NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/ISF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/ISF/isfcpl.F90 @ 13630

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

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

Download in other formats: