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_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ISF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ISF/isfcpl.F90 @ 14644

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

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

Download in other formats: