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_r12377_KERNEL-06_techene_e3/src/OCE/ISF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90 @ 12724

Last change on this file since 12724 was 12724, checked in by techene, 4 years ago
branch KERNEL-06 : merge with trunk@12698 #2385 - in duplcated files : changes to comply to the new trunk variables and some loop bug fixes
File size: 36.3 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 domqe , ONLY: dom_qe_zgr ! vertical scale factor interpolation
21	#endif
22	USE domngb , 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 ! e3t , e3u, e3v
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	END DO
152	!
153	IF( lwxios ) CALL iom_swap( cwxios_context )
154	CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask , ldxios = lwxios )
155	CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
156	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n' , ze3t , ldxios = lwxios )
157	CALL iom_rstput( kt, nitrst, numrow, 'e3u_n' , ze3u , ldxios = lwxios )
158	CALL iom_rstput( kt, nitrst, numrow, 'e3v_n' , ze3v , ldxios = lwxios )
159	CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw(:,:,:,Kmm) , ldxios = lwxios )
160	IF( lwxios ) CALL iom_swap( cxios_context )
161	!
162	END SUBROUTINE isfcpl_rst_write
163
164	SUBROUTINE isfcpl_ssh(Kbb, Kmm, Kaa)
165	!!----------------------------------------------------------------------
166	!! * ROUTINE iscpl_ssh *
167	!!
168	!! ** Purpose : basic guess of ssh in new wet cell
169	!!
170	!! ** Method : basic extrapolation from neigbourg cells
171	!!
172	!!----------------------------------------------------------------------
173	!!
174	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
175	!!----------------------------------------------------------------------
176	INTEGER :: ji, jj, jd, jk !! loop index
177	INTEGER :: jip1, jim1, jjp1, jjm1
178	!!
179	REAL(wp):: zsummsk
180	REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
181	!!----------------------------------------------------------------------
182	!
183	CALL iom_get( numror, jpdom_autoglo, 'ssmask' , zssmask_b, ldxios = lrxios ) ! need to extrapolate T/S
184
185	! compute new ssh if we open a full water column
186	! rude average of the closest neigbourgs (e1e2t not taking into account)
187	!
188	zssh(:,:) = ssh(:,:,Kmm)
189	zssmask0(:,:) = zssmask_b(:,:)
190	!
191	DO jd = 1, nn_drown
192	!
193	zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
194	DO_2D_00_00
195	jip1=ji+1; jim1=ji-1;
196	jjp1=jj+1; jjm1=jj-1;
197	!
198	zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
199	!
200	IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
201	ssh(ji,jj,Kmm)=( zssh(jip1,jj)*zssmask0(jip1,jj) &
202	& + zssh(jim1,jj)*zssmask0(jim1,jj) &
203	& + zssh(ji,jjp1)*zssmask0(ji,jjp1) &
204	& + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
205	zssmask_b(ji,jj) = 1._wp
206	ENDIF
207	END_2D
208	!
209	zssh(:,:) = ssh(:,:,Kmm)
210	zssmask0(:,:) = zssmask_b(:,:)
211	!
212	CALL lbc_lnk_multi( 'iscplrst', zssh, 'T', 1., zssmask0, 'T', 1. )
213	!
214	END DO
215	!
216	! update ssh(:,:,Kmm)
217	ssh(:,:,Kmm) = zssh(:,:) * ssmask(:,:)
218	!
219	ssh(:,:,Kbb) = ssh(:,:,Kmm)
220	!
221	IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',ssh(:,:,Kmm))
222	!
223	! recompute the vertical scale factor, depth and water thickness
224	IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
225	IF(lwp) write(numout,*) '~~~~~~~~~~~'
226	#if ! defined key_qco
227	DO jk = 1, jpk
228	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + (ht_0(:,:) + ssh(:,:,Kmm)) * r1_ht_0(:,:) )
229	END DO
230	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
231	CALL dom_vvl_zgr(Kbb, Kmm, Kaa)
232	#else
233	CALL dom_qe_zgr(Kbb, Kmm, Kaa)
234	#endif
235	!
236	END SUBROUTINE isfcpl_ssh
237
238	SUBROUTINE isfcpl_tra(Kmm)
239	!!----------------------------------------------------------------------
240	!! * ROUTINE iscpl_tra *
241	!!
242	!! ** Purpose : compute new tn, sn in case of evolving geometry of ice shelves
243	!!
244	!! ** Method : tn, sn : basic extrapolation from neigbourg cells
245	!!
246	!!----------------------------------------------------------------------
247	INTEGER, INTENT(in) :: Kmm ! ocean time level index
248	!!----------------------------------------------------------------------
249	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
250	!REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b !! depth w before
251	!!
252	INTEGER :: ji, jj, jk, jd !! loop index
253	INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
254	!!
255	REAL(wp):: zsummsk
256	REAL(wp):: zdz, zdzm1, zdzp1
257	!!
258	REAL(wp), DIMENSION(jpi,jpj) :: zdmask
259	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn
260	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d
261	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
262	!!----------------------------------------------------------------------
263	!
264	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b, ldxios = lrxios ) ! need to extrapolate T/S
265	!CALL iom_get( numror, jpdom_autoglo, 'wmask' , zwmask_b, ldxios = lrxios ) ! need to extrapolate T/S
266	!CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
267	!
268	!
269	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
270	!PM: Is this IF needed since change to VVL by default
271	!bugged : to be corrected (PM)
272	! back up original t/s/mask
273	!tsb (:,:,:,:) = ts(:,:,:,:,Kmm)
274	!
275	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
276
277	! IF (.NOT.ln_linssh) THEN
278	! DO jk = 2,jpk-1
279	! DO jj = 1,jpj
280	! DO ji = 1,jpi
281	! 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
282	!
283	! !compute weight
284	! zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw(ji,jj,jk+1,Kmm))
285	! zdzm1 = MAX(0._wp,gdepw(ji,jj,jk ,Kmm) - pdepw_b(ji,jj,jk ))
286	! zdz = e3t(ji,jj,jk,Kmm) - zdzp1 - zdzm1 ! if isf : e3t = gdepw(ji,jj,jk+1,Kmm)- gdepw(ji,jj,jk,Kmm)
287	!
288	! IF (zdz .LT. 0._wp) THEN
289	! CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
290	! END IF
291	!
292	! ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb) &
293	! & + zdz *ts(ji,jj,jk ,jp_tem,Kbb) &
294	! & + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb) )/e3t(ji,jj,jk,Kmm)
295	!
296	! ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb) &
297	! & + zdz *ts(ji,jj,jk ,jp_sal,Kbb) &
298	! & + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb) )/e3t(ji,jj,jk,Kmm)
299	!
300	! END IF
301	! END DO
302	! END DO
303	! END DO
304	! END IF
305
306	zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
307	ztmask0(:,:,:) = ztmask_b(:,:,:)
308	ztmask1(:,:,:) = ztmask_b(:,:,:)
309	!
310	! iterate the extrapolation processes nn_drown times
311	DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
312	DO jk = 1,jpk-1
313	!
314	! define new wet cell
315	zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
316	!
317	DO_2D_00_00
318	jip1=ji+1; jim1=ji-1;
319	jjp1=jj+1; jjm1=jj-1;
320	!
321	! check if a wet neigbourg cell is present
322	zsummsk = ztmask0(jip1,jj ,jk) + ztmask0(jim1,jj ,jk) &
323	+ ztmask0(ji ,jjp1,jk) + ztmask0(ji ,jjm1,jk)
324	!
325	! if neigbourg wet cell available at the same level
326	IF ( zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
327	!
328	! horizontal basic extrapolation
329	ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj ,jk,1) * ztmask0(jip1,jj ,jk) &
330	& + zts0(jim1,jj ,jk,1) * ztmask0(jim1,jj ,jk) &
331	& + zts0(ji ,jjp1,jk,1) * ztmask0(ji ,jjp1,jk) &
332	& + zts0(ji ,jjm1,jk,1) * ztmask0(ji ,jjm1,jk) ) / zsummsk
333	ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj ,jk,2) * ztmask0(jip1,jj ,jk) &
334	& + zts0(jim1,jj ,jk,2) * ztmask0(jim1,jj ,jk) &
335	& + zts0(ji ,jjp1,jk,2) * ztmask0(ji ,jjp1,jk) &
336	& + zts0(ji ,jjm1,jk,2) * ztmask0(ji ,jjm1,jk) ) / zsummsk
337	!
338	! update mask for next pass
339	ztmask1(ji,jj,jk)=1
340	!
341	! in case no neigbourg wet cell available at the same level
342	! check if a wet cell is available below
343	ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
344	!
345	! vertical extrapolation if horizontal extrapolation failed
346	jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
347	!
348	! check if a wet neigbourg cell is present
349	zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
350	IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
351	ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) &
352	& + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
353	ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) &
354	& + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
355	!
356	! update mask for next pass
357	ztmask1(ji,jj,jk)=1._wp
358	END IF
359	END IF
360	END_2D
361	END DO
362	!
363	! update temperature and salinity and mask
364	zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
365	ztmask0(:,:,:) = ztmask1(:,:,:)
366	!
367	CALL lbc_lnk_multi( 'iscplrst', zts0(:,:,:,jp_tem), 'T', 1., zts0(:,:,:,jp_sal), 'T', 1., ztmask0, 'T', 1.)
368	!
369	END DO ! nn_drown
370	!
371	! mask new ts(:,:,:,:,Kmm) field
372	ts(:,:,:,jp_tem,Kmm) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
373	ts(:,:,:,jp_sal,Kmm) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
374	!
375	! sanity check
376	! -----------------------------------------------------------------------------------------
377	! case we open a cell but no neigbour cells available to get an estimate of T and S
378	DO_3D_11_11( 1,jpk-1 )
379	IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp) &
380	& CALL ctl_stop('STOP', 'failing to fill all new weet cell, &
381	& try increase nn_drown or activate XXXX &
382	& in your domain cfg computation' )
383	END_3D
384	!
385	END SUBROUTINE isfcpl_tra
386
387	SUBROUTINE isfcpl_vol(Kmm)
388	!!----------------------------------------------------------------------
389	!! * ROUTINE iscpl_vol *
390	!!
391	!! ** Purpose : compute the correction of the local divergence to apply
392	!! during the first time step after the coupling.
393	!!
394	!! ** Method : - compute horizontal vol div. before/after coupling
395	!! - compute vertical input
396	!! - compute correction
397	!!
398	!!----------------------------------------------------------------------
399	!!
400	INTEGER, INTENT(in) :: Kmm ! ocean time level index
401	!!----------------------------------------------------------------------
402	INTEGER :: ji, jj, jk
403	INTEGER :: ikb, ikt
404	!!
405	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln ! vol flux div. before/after coupling
406	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b ! vertical scale factor before/after coupling
407	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before coupling
408	!!----------------------------------------------------------------------
409	!
410	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b, ldxios = lrxios )
411	CALL iom_get( numror, jpdom_autoglo, 'e3u_n' , ze3u_b , ldxios = lrxios )
412	CALL iom_get( numror, jpdom_autoglo, 'e3v_n' , ze3v_b , ldxios = lrxios )
413	!
414	! 1.0: compute horizontal volume flux divergence difference before-after coupling
415	!
416	DO jk = 1, jpk ! Horizontal slab
417	! 1.1: get volume flux before coupling (>0 out)
418	DO_2D_00_00
419	zqvolb(ji,jj,jk) = &
420	& ( e2u(ji ,jj ) * ze3u_b(ji ,jj ,jk) * uu(ji ,jj ,jk,Kmm) &
421	& - e2u(ji-1,jj ) * ze3u_b(ji-1,jj ,jk) * uu(ji-1,jj ,jk,Kmm) &
422	& + e1v(ji ,jj ) * ze3v_b(ji ,jj ,jk) * vv(ji ,jj ,jk,Kmm) &
423	& - e1v(ji ,jj-1) * ze3v_b(ji ,jj-1,jk) * vv(ji ,jj-1,jk,Kmm) ) &
424	& * ztmask_b(ji,jj,jk)
425	END_2D
426	!
427	! 1.2: get volume flux after coupling (>0 out)
428	! properly mask velocity
429	! (velocity are still mask with old mask at this stage)
430	uu(:,:,jk,Kmm) = uu(:,:,jk,Kmm) * umask(:,:,jk)
431	vv(:,:,jk,Kmm) = vv(:,:,jk,Kmm) * vmask(:,:,jk)
432	! compute volume flux divergence after coupling
433	DO_2D_00_00
434	zqvoln(ji,jj,jk) = &
435	& ( e2u(ji ,jj ) * e3u(ji ,jj ,jk,Kmm) * uu(ji ,jj ,jk,Kmm) &
436	& - e2u(ji-1,jj ) * e3u(ji-1,jj ,jk,Kmm) * uu(ji-1,jj ,jk,Kmm) &
437	& + e1v(ji ,jj ) * e3v(ji ,jj ,jk,Kmm) * vv(ji ,jj ,jk,Kmm) &
438	& - e1v(ji ,jj-1) * e3v(ji ,jj-1,jk,Kmm) * vv(ji ,jj-1,jk,Kmm) ) &
439	& * tmask(ji,jj,jk)
440	END_2D
441	!
442	! 1.3: get 3d volume flux difference (before - after cpl) (>0 out)
443	! correction to add is _b - _n
444	risfcpl_vol(:,:,jk) = zqvolb(:,:,jk) - zqvoln(:,:,jk)
445	END DO
446	!
447	! 2.0: include the contribution of the vertical velocity in the volume flux correction
448	!
449	DO_2D_00_00
450	!
451	ikt = mikt(ji,jj)
452	IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN
453	risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1)) ! test sign
454	ENDIF
455	!
456	END_2D
457	!
458	CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
459	!
460	! 3.0: set total correction (div, tr(:,:,:,:,Krhs), ssh)
461	!
462	! 3.1: mask volume flux divergence correction
463	risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
464	!
465	! 3.2: get 3d tr(:,:,:,:,Krhs) increment to apply at the first time step
466	! temperature and salt content flux computed using local ts(:,:,:,:,Kmm)
467	! (very simple advection scheme)
468	! (>0 out)
469	risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_tem,Kmm)
470	risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_sal,Kmm)
471	!
472	! 3.3: ssh correction (for dynspg_ts)
473	risfcpl_ssh(:,:) = 0.0
474	DO jk = 1,jpk
475	risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
476	END DO
477
478	END SUBROUTINE isfcpl_vol
479
480	SUBROUTINE isfcpl_cons(Kmm)
481	!!----------------------------------------------------------------------
482	!! * ROUTINE iscpl_cons *
483	!!
484	!! ** Purpose : compute the corrective increment in volume/salt/heat to put back the vol/heat/salt
485	!! removed or added during the coupling processes (wet or dry new cell)
486	!!
487	!! ** Method : - compare volume/heat/salt before and after
488	!! - look for the closest wet cells (share amoung neigbourgs if there are)
489	!! - build the correction increment to applied at each time step
490	!!
491	!!----------------------------------------------------------------------
492	!
493	TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
494	!
495	!!----------------------------------------------------------------------
496	INTEGER, INTENT(in) :: Kmm ! ocean time level index
497	!!----------------------------------------------------------------------
498	INTEGER :: ji , jj , jk , jproc ! loop index
499	INTEGER :: jip1 , jim1, jjp1, jjm1 ! dummy indices
500	INTEGER :: iig , ijg, ik ! dummy indices
501	INTEGER :: jisf ! start, end and current position in the increment array
502	INTEGER :: ingb, ifind ! 0/1 target found or need to be found
503	INTEGER :: nisfl_area ! global number of cell concerned by the wet->dry case
504	INTEGER, DIMENSION(jpnij) :: nisfl ! local number of cell concerned by the wet->dry case
505	!
506	REAL(wp) :: z1_sum, z1_rdtiscpl
507	REAL(wp) :: zdtem, zdsal, zdvol, zratio ! tem, sal, vol increment
508	REAL(wp) :: zlon , zlat ! target location
509	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before
510	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b ! scale factor before
511	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b ! scale factor before
512	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b ! scale factor before
513	!!----------------------------------------------------------------------
514
515	!==============================================================================
516	! 1.0: initialisation
517	!==============================================================================
518
519	! get restart variable
520	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b(:,:,:), ldxios = lrxios ) ! need to extrapolate T/S
521	CALL iom_get( numror, jpdom_autoglo, 'e3t_n' , ze3t_b(:,:,:) , ldxios = lrxios )
522	CALL iom_get( numror, jpdom_autoglo, 'tn' , zt_b(:,:,:) , ldxios = lrxios )
523	CALL iom_get( numror, jpdom_autoglo, 'sn' , zs_b(:,:,:) , ldxios = lrxios )
524
525	! compute run length
526	nstp_iscpl = nitend - nit000 + 1
527	rdt_iscpl = nstp_iscpl * rn_Dt
528	z1_rdtiscpl = 1._wp / rdt_iscpl
529
530	IF (lwp) WRITE(numout,*) ' nb of stp for cons = ', nstp_iscpl
531	IF (lwp) WRITE(numout,*) ' coupling time step = ', rdt_iscpl
532
533	! initialisation correction
534	risfcpl_cons_vol = 0.0
535	risfcpl_cons_ssh = 0.0
536	risfcpl_cons_tsc = 0.0
537
538	!==============================================================================
539	! 2.0: diagnose the heat, salt and volume input and compute the correction variable
540	! for case where we wet a cell or cell still wet (no change in cell status)
541	!==============================================================================
542
543	DO jk = 1,jpk-1
544	DO jj = nldj,nlej
545	DO ji = nldi,nlei
546
547	! volume diff
548	zdvol = e3t (ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
549	& - ze3t_b(ji,jj,jk ) * ztmask_b(ji,jj,jk)
550
551	! heat diff
552	zdtem = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
553	- zt_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
554
555	! salt diff
556	zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
557	- zs_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
558
559	! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
560	! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
561	risfcpl_cons_vol(ji,jj,jk) = ( zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk) ) * z1_rdtiscpl
562	risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
563	risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
564
565	END DO
566	END DO
567	END DO
568	!
569	!==============================================================================
570	! 3.0: diagnose the heat, salt and volume input and compute the correction variable
571	! for case where we close a cell
572	!==============================================================================
573	!
574	! compute the total number of point receiving a correction increment for each processor
575	! local
576	nisfl(:)=0
577	DO jk = 1,jpk-1
578	DO jj = nldj,nlej
579	DO ji = nldi,nlei
580	jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
581	IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
582	nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
583	ENDIF
584	ENDDO
585	ENDDO
586	ENDDO
587	!
588	! global
589	CALL mpp_sum('isfcpl',nisfl )
590	!
591	! allocate list of point receiving correction
592	ALLOCATE(zisfpts(nisfl(narea)))
593	!
594	zisfpts(:) = isfcons(0,0,0,-HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), 0)
595	!
596	! start computing the correction and fill zisfpts
597	! local
598	jisf = 0
599	DO jk = 1,jpk-1
600	DO jj = nldj,nlej
601	DO ji = nldi,nlei
602	IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
603
604	jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
605
606	zdvol = risfcpl_cons_vol(ji,jj,jk )
607	zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
608	zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
609
610	IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
611	! spread correction amoung neigbourg wet cells (horizontal direction first)
612	! as it is a rude correction corner and lateral cell have the same weight
613	!
614	z1_sum = 1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
615	!
616	! lateral cells
617	IF (tmask(jip1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
618	IF (tmask(jim1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
619	IF (tmask(ji ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
620	IF (tmask(ji ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
621	!
622	! corner cells
623	IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
624	IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
625	IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
626	IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
627	!
628	ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
629	! spread correction amoung neigbourg wet cells (vertical direction)
630	CALL update_isfpts(zisfpts, jisf, ji , jj , jk+1, zdvol, zdsal, zdtem, 1., 0)
631	ELSE
632	! need to find where to put correction in later on
633	CALL update_isfpts(zisfpts, jisf, ji , jj , jk , zdvol, zdsal, zdtem, 1., 1)
634	END IF
635	END IF
636	END DO
637	END DO
638	END DO
639	!
640	! share data among all processes because for some point we need to find the closest wet point (could be on other process)
641	DO jproc=1,jpnij
642	!
643	! share total number of isf point treated for proc jproc
644	IF (jproc==narea) THEN
645	nisfl_area=nisfl(jproc)
646	ELSE
647	nisfl_area=0
648	END IF
649	CALL mpp_max('isfcpl',nisfl_area)
650	!
651	DO jisf = 1,nisfl_area
652	!
653	IF (jproc==narea) THEN
654	! indices (conversion to global indices and sharing)
655	iig = zisfpts(jisf)%ii ; ijg = zisfpts(jisf)%jj ; ik = zisfpts(jisf)%kk
656	!
657	! data
658	zdvol = zisfpts(jisf)%dvol ; zdsal = zisfpts(jisf)%dsal ; zdtem = zisfpts(jisf)%dtem
659	!
660	! location
661	zlat = zisfpts(jisf)%lat ; zlon = zisfpts(jisf)%lon
662	!
663	! find flag
664	ingb = zisfpts(jisf)%ngb
665	ELSE
666	iig =0 ; ijg =0 ; ik =0
667	zdvol=-HUGE(1.0) ; zdsal=-HUGE(1.0) ; zdtem=-HUGE(1.0)
668	zlat =-HUGE(1.0) ; zlon =-HUGE(1.0)
669	ingb = 0
670	END IF
671	!
672	! share data (need synchronisation of data as get_correction call a global com)
673	CALL mpp_max('isfcpl',iig) ; CALL mpp_max('isfcpl',ijg) ; CALL mpp_max('isfcpl',ik)
674	CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
675	CALL mpp_max('isfcpl',zlat) ; CALL mpp_max('isfcpl',zlon)
676	CALL mpp_max('isfcpl',ingb)
677	!
678	! fill the 3d correction array
679	CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
680	END DO
681	END DO
682	!
683	!==============================================================================
684	! 4.0: finalisation and compute ssh equivalent of the volume correction
685	!==============================================================================
686	!
687	! mask (>0 out)
688	risfcpl_cons_vol(:,:,: ) = risfcpl_cons_vol(:,:,: ) * tmask(:,:,:)
689	risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
690	risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
691	!
692	! add lbclnk
693	CALL lbc_lnk_multi( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1., risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1., &
694	& risfcpl_cons_vol(:,:,:) , 'T', 1.)
695	!
696	! ssh correction (for dynspg_ts)
697	DO jk = 1,jpk
698	risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
699	END DO
700	risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
701	!
702	END SUBROUTINE isfcpl_cons
703	!
704	SUBROUTINE update_isfpts(sisfpts, kpts, ki, kj, kk, pdvol, pdsal, pdtem, pratio, kfind)
705	!!---------------------------------------------------------------------
706	!! * ROUTINE update_isfpts *
707	!!
708	!! ** Purpose : if a cell become dry, we need to put the corrective increment elsewhere
709	!!
710	!! ** Action : update the list of point
711	!!
712	!!----------------------------------------------------------------------
713	!!----------------------------------------------------------------------
714	TYPE(isfcons), DIMENSION(:), INTENT(inout) :: sisfpts
715	INTEGER, INTENT(inout) :: kpts
716	!!----------------------------------------------------------------------
717	INTEGER, INTENT(in ) :: ki, kj, kk ! target location (kfind=0)
718	! ! or source location (kfind=1)
719	INTEGER, INTENT(in ), OPTIONAL :: kfind ! 0 target cell already found
720	! ! 1 target to be determined
721	REAL(wp), INTENT(in ) :: pdvol, pdsal, pdtem, pratio ! vol/sal/tem increment
722	! ! and ratio in case increment span over multiple cells.
723	!!----------------------------------------------------------------------
724	INTEGER :: ifind
725	!!----------------------------------------------------------------------
726	!
727	! increment position
728	kpts = kpts + 1
729	!
730	! define if we need to look for closest valid wet cell (no neighbours or neigbourg on halo)
731	IF ( PRESENT(kfind) ) THEN
732	ifind = kfind
733	ELSE
734	ifind = ( 1 - tmask_h(ki,kj) ) * tmask(ki,kj,kk)
735	END IF
736	!
737	! update isfpts structure
738	sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
739	!
740	END SUBROUTINE update_isfpts
741	!
742	SUBROUTINE get_correction( ki, kj, kk, plon, plat, pvolinc, psalinc, pteminc, kfind)
743	!!---------------------------------------------------------------------
744	!! * ROUTINE get_correction *
745	!!
746	!! ** Action : - Find the closest valid cell if needed (wet and not on the halo)
747	!! - Scale the correction depending of pratio (case where multiple wet neigbourgs)
748	!! - Fill the correction array
749	!!
750	!!----------------------------------------------------------------------
751	INTEGER , INTENT(in) :: ki, kj, kk, kfind ! target point indices
752	REAL(wp), INTENT(in) :: plon, plat ! target point lon/lat
753	REAL(wp), INTENT(in) :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
754	!!----------------------------------------------------------------------
755	INTEGER :: jj, ji, iig, ijg
756	!!----------------------------------------------------------------------
757	!
758	! define global indice of correction location
759	iig = ki ; ijg = kj
760	IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
761	!
762	! fill the correction array
763	DO jj = mj0(ijg),mj1(ijg)
764	DO ji = mi0(iig),mi1(iig)
765	! correct the vol_flx and corresponding heat/salt flx in the closest cell
766	risfcpl_cons_vol(ji,jj,kk) = risfcpl_cons_vol(ji,jj,kk ) + pvolinc
767	risfcpl_cons_tsc(ji,jj,kk,jp_sal) = risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
768	risfcpl_cons_tsc(ji,jj,kk,jp_tem) = risfcpl_cons_tsc(ji,jj,kk,jp_tem) + pteminc
769	END DO
770	END DO
771
772	END SUBROUTINE get_correction
773
774	END MODULE isfcpl

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

Download in other formats: