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.

sbciscpl.F90 in branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbciscpl.F90 @ 5619

Last change on this file since 5619 was 5619, checked in by mathiot, 9 years ago
ocean/ice sheet coupling: initial commit
File size: 42.2 KB

Line
1	MODULE sbc_iscpl
2	!!======================================================================
3	!! * MODULE sbciscpl*
4	!! Ocean forcing: river runoff
5	!!=====================================================================
6	!! History : NEMO ! 2015-01 P. Mathiot: original
7	!!----------------------------------------------------------------------
8
9	!!----------------------------------------------------------------------
10	!! sbc_iscpl_alloc: variable allocation
11	!! rst_iscpl : restart correction in case of coupling with ice sheet
12	!! sbc_iscpl_div : correction of divergence to keep volume conservation
13	!!----------------------------------------------------------------------
14	USE dom_oce ! ocean space and time domain
15	USE domwri ! ocean space and time domain
16	USE domvvl, ONLY : dom_vvl_interpol
17	USE phycst ! physical constants
18	USE sbc_oce ! surface boundary condition variables
19	USE oce ! global tra/dyn variable
20	USE in_out_manager ! I/O manager
21	USE iom ! I/O module
22	USE lib_mpp ! MPP library
23	USE lib_fortran ! MPP library
24	USE wrk_nemo ! Memory allocation
25	USE lbclnk
26	USE domngb
27	USE iom
28	USE sbc_ice, ONLY : lk_lim3
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC sbc_iscpl_div
34	PUBLIC rst_iscpl
35	PUBLIC rst_iscpl_interpol ! routine to wet and dry
36	PUBLIC rst_iscpl_alloc
37	!! !!* namsbc_iscpl namelist *
38	LOGICAL , PUBLIC :: ln_hfb
39	REAL(wp), PUBLIC :: rn_fiscpl
40	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,: ) :: hdiv_iscpl
41	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: htsc_iscpl
42	!! * Substitutions
43	# include "domzgr_substitute.h90"
44	!!----------------------------------------------------------------------
45	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
46	!! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
48	!!----------------------------------------------------------------------
49	CONTAINS
50
51	INTEGER FUNCTION rst_iscpl_alloc()
52	!!----------------------------------------------------------------------
53	!! * ROUTINE sbc_iscpl_alloc *
54	!!----------------------------------------------------------------------
55	ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=rst_iscpl_alloc )
56	!
57	IF( lk_mpp ) CALL mpp_sum ( rst_iscpl_alloc )
58	IF( rst_iscpl_alloc > 0 ) CALL ctl_warn('rst_iscpl_alloc: allocation of arrays failed')
59	END FUNCTION rst_iscpl_alloc
60
61	SUBROUTINE rst_iscpl
62	REAL(wp), DIMENSION(:,: ), POINTER :: zsmask_b
63	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b
64	REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b , ze3u_b , ze3v_b
65	REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b
66	INTEGER :: jk
67	LOGICAL :: llok
68	!!----------------------------------------------------------------------
69	INTEGER :: inum0
70	CHARACTER(20) :: cfile
71
72
73	WRITE(numout,*) 'rst_read ln_iscpl : ', ln_iscpl
74	IF ( ln_iscpl ) THEN
75	CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! need this new variable for evoving isf cavity
76	CALL wrk_alloc(jpi,jpj,jpk, ze3t_b , ze3u_b , ze3v_b ) ! need of this variable afterwards (only used for interpolation)
77	CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
78	CALL wrk_alloc(jpi,jpj, zsmask_b )
79
80	IF( rst_iscpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
81
82	WRITE(numout,*) 'rst_read : modify restart to fit new geometry'
83	WRITE(numout,*) '~~~~~~~~'
84
85	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b ) ! need to extrapolate T/S
86	CALL iom_get( numror, jpdom_autoglo, 'umask' , zumask_b ) ! need to correct barotropic velocity
87	CALL iom_get( numror, jpdom_autoglo, 'vmask' , zvmask_b ) ! need to correct barotropic velocity
88	CALL iom_get( numror, jpdom_autoglo, 'smask' , zsmask_b ) ! need to correct barotropic velocity
89	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', ze3t_b(:,:,:) ) ! need to compute temperature correction
90	CALL iom_get( numror, jpdom_autoglo, 'fse3u_n', ze3u_b(:,:,:) ) ! need to compute volume correction ????
91	CALL iom_get( numror, jpdom_autoglo, 'fse3v_n', ze3v_b(:,:,:) ) ! need to compute volume correction ????
92	CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to compute volume correction ????
93	!! ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
94	CALL rst_iscpl_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b, htsc_iscpl, hdiv_iscpl )
95	IF( nmsh /= 0 .AND. ln_iscpl ) CALL dom_wri ! Create a domain file
96
97	cfile='correction'
98	cfile = TRIM( cfile )
99	CALL iom_open ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
100	CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
101	CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
102	CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
103	CALL iom_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)tmask(:,:,:) - ze3t_b(:,:,:)ztmask_b(:,:,:))*tmask(:,:,:) )
104	CALL iom_close ( inum0 )
105
106	CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b ) ! no need of this variable afterwards (only used for interpolation)
107	CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b ,ze3u_b ,ze3v_b ) ! no need of this variable afterwards (only used for interpolation)
108	CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b )
109	CALL wrk_dealloc(jpi,jpj, zsmask_b )
110	END IF
111	IF( neuler == 0 ) THEN ! Euler restart (neuler=0)
112	tsb (:,:,:,:) = tsn (:,:,:,:) ! all before fields set to now values
113	ub (:,:,:) = un (:,:,:)
114	vb (:,:,:) = vn (:,:,:)
115	rotb (:,:,:) = rotn (:,:,:)
116	hdivb(:,:,:) = hdivn(:,:,:)
117	sshb (:,:) = sshn (:,:)
118
119	IF( lk_vvl ) THEN
120	DO jk = 1, jpk
121	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
122	END DO
123	ENDIF
124
125	IF( lk_lim3 .AND. .NOT. lk_vvl ) THEN
126	DO jk = 1, jpk
127	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
128	END DO
129	ENDIF
130
131	ENDIF
132	!
133	IF( lk_lim3 ) THEN
134	CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev )
135	ENDIF
136
137	END SUBROUTINE rst_iscpl
138
139	SUBROUTINE rst_iscpl_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b, pts_flx, pvol_flx)
140	!!----------------------------------------------------------------------
141	!! * ROUTINE rst_iscpl *
142	!!
143	!! ** Purpose : compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
144	!! compute 2d fields of heat, salt and volume correction
145	!!
146	!! ** Method : tn, sn : extrapolation from neigbourg cells
147	!! un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
148	!!----------------------------------------------------------------------
149	!!######## ~ cp/paste pyton script for this routine
150	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b !! mask before
151	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b , pe3u_b , pe3v_b !! scale factor before
152	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b
153	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: psmask_b !! mask before
154	REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx !! corrective flux to have tracer conservation
155	REAL(wp), DIMENSION(:,:,: ), INTENT(out) :: pvol_flx !! corrective flux to have volume conservation
156	!!
157	INTEGER :: ji, jj, ii, ij, jk, iz !! loop index
158	INTEGER :: iip1, iim1, ijp1, ijm1, jkp1, jkm1
159	INTEGER :: ierror, inum ! temporary integer
160	INTEGER :: ios ! Local integer output status for namelist read
161	!!
162	REAL(wp):: summsk, r1_tiscpl, zsum, zsum1, zarea, zsumn, zsumb
163	REAL(wp):: ziip1_ratio, ziim1_ratio, zijp1_ratio, zijm1_ratio
164	REAL(wp):: zdz, zdzm1, zdzp1
165	!!
166	REAL(wp), DIMENSION(:,: ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
167	REAL(wp), DIMENSION(:,: ), POINTER :: zbub , zbvb , zbun , zbvn
168	REAL(wp), DIMENSION(:,: ), POINTER :: zssh0 , zssh1, hu1, hv1
169	REAL(wp), DIMENSION(:,: ), POINTER :: zsmask0, zsmask1
170	REAL(wp), DIMENSION(:,:,: ), POINTER :: ztmask0, ztmask1, ztrp
171	REAL(wp), DIMENSION(:,:,: ), POINTER :: zwmaskn, zwmaskb, ztmp3d
172	REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
173	!
174	REAL(wp), DIMENSION(: ), ALLOCATABLE :: zlon, zlat
175	REAL(wp), DIMENSION(: ), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
176	INTEGER , DIMENSION(: ), ALLOCATABLE :: ixpts, iypts, izpts, vnpts
177	INTEGER :: jpts, npts
178	!
179	NAMELIST/namsbc_iscpl/rn_fiscpl,ln_hfb
180	!!----------------------------------------------------------------------
181	! ! ============
182	! ! Namelist
183	! ! ============
184	!
185	rn_fiscpl = 2480
186	ln_hfb = .FALSE.
187	REWIND( numnam_ref ) ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling
188	READ ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901)
189	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in reference namelist', lwp )
190
191	REWIND( numnam_cfg ) ! Namelist namsbc_iscpl in configuration namelist : Ice Sheet coupling
192	READ ( numnam_cfg, namsbc_iscpl, IOSTAT = ios, ERR = 902 )
193	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in configuration namelist', lwp )
194	IF(lwm) WRITE ( numond, namsbc_iscpl )
195	!
196	PRINT *, 'rn_iscpl .. =',rn_fiscpl, ln_hfb
197
198	CALL wrk_alloc(jpi,jpj,jpk,2, zts0 )
199	CALL wrk_alloc(jpi,jpj,jpk, ztmask0, ztmask1 , ztrp, ztmp3d )
200	CALL wrk_alloc(jpi,jpj,jpk, zwmaskn, zwmaskb )
201	CALL wrk_alloc(jpi,jpj, zsmask0, zsmask1 )
202	CALL wrk_alloc(jpi,jpj, zdmask , zdsmask, zvcorr, zucorr, zde3t)
203	CALL wrk_alloc(jpi,jpj, zbub , zbvb , zbun , zbvn )
204	CALL wrk_alloc(jpi,jpj, zssh0 , zssh1, hu1, hv1 )
205	! mask value to be sure
206	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
207	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
208	!=============================================================================
209
210	! compute wmask
211	zwmaskn(:,:,1) = tmask (:,:,1)
212	zwmaskb(:,:,1) = ptmask_b(:,:,1)
213	DO jk = 2,jpk
214	zwmaskn(:,:,jk) = tmask (:,:,jk) * tmask (:,:,jk-1)
215	zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
216	END DO
217
218
219	!=============================================================================
220
221	! compute ssh
222	zsmask0(:,:) = psmask_b(:,:)
223	zsmask1(:,:) = psmask_b(:,:)
224
225	! compute new ssh if we open a full water column (mean of the closest neigbourgs)
226	sshb (:,:)=sshn(:,:)
227	zssh0(:,:)=sshn(:,:)
228	DO iz = 1,10 ! need to be tuned (configuration dependent)
229	zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
230	DO ii = 2,jpi-1
231	DO ij = 2,jpj-1
232	iip1=ii+1; iim1=ii-1;
233	ijp1=ij+1; ijm1=ij-1;
234	summsk=(zsmask0(iip1,ij)+zsmask0(iim1,ij)+zsmask0(ii,ijp1)+zsmask0(ii,ijm1))
235	IF (zdsmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
236	PRINT *, 'add ssh to 1 cell',ii,ij,narea
237	sshn(ii,ij)=(zssh0(iip1,ij)*zsmask0(iip1,ij) &
238	& +zssh0(iim1,ij)*zsmask0(iim1,ij) &
239	& +zssh0(ii,ijp1)*zsmask0(ii,ijp1) &
240	& +zssh0(ii,ijm1)*zsmask0(ii,ijm1))/summsk
241	zsmask1(ii,ij)=1
242	END IF
243	END DO
244	END DO
245	CALL lbc_lnk(sshn,'T',1.)
246	CALL lbc_lnk(zsmask1,'T',1.)
247	zssh0 = sshn
248	zsmask0 = zsmask1
249	END DO
250	sshn(:,:) = sshn(:,:) * ssmask(:,:)
251
252	!=============================================================================
253	! WARNING we cannot used glob_sum for before time variable (because glob_sum use tmask_i). Need to mask the halo and glob_sum_full
254	ztmp3d(:,:,:) = 0.0
255	ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
256	& - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
257	zsum = glob_sum_full(ztmp3d)
258	IF (lwp) PRINT *, 'total volume correction 00 = ',zsum
259	IF ( lk_vvl ) THEN
260	! compute fse3t_n
261	DO jk = 1,jpk
262	fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1 + sshn(:,:) / ( ht_0(:,:) + 1. - ssmask(:,:) ) * tmask(:,:,jk) )
263	END DO
264	ztmp3d(:,:,:) = 0.0
265	ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
266	& - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
267	zsum = glob_sum_full(ztmp3d)
268	IF (lwp) PRINT *, 'total volume correction 01 = ',zsum
269
270	! compute fse3u/v ... (call interpolation vvl)
271	! Reconstruction of all vertical scale factors at now and before time steps
272	! =============================================================================
273	! Horizontal scale factor interpolations
274	! --------------------------------------
275	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
276	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
277	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
278	! Vertical scale factor interpolations
279	! ------------------------------------
280	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
281	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
282	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
283	! t- and w- points depth
284	! ----------------------
285	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
286	fsdepw_n(:,:,1) = 0.0_wp
287	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
288	DO jj = 1,jpj
289	DO ji = 1,jpi
290	DO jk = 2,mikt(ji,jj)-1
291	fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
292	fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
293	fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
294	END DO
295	IF (mikt(ji,jj) .GT. 1) THEN
296	jk = mikt(ji,jj)
297	fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
298	fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
299	fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj)
300	END IF
301	DO jk = mikt(ji,jj)+1, jpk
302	fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
303	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
304	fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj)
305	END DO
306	END DO
307	END DO
308
309	hu(:,:) = 0._wp ! Ocean depth at U-points
310	hv(:,:) = 0._wp ! Ocean depth at V-points
311	ht(:,:) = 0._wp ! Ocean depth at T-points
312	DO jk = 1, jpkm1
313	hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
314	hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
315	ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
316	END DO
317	! ! Inverse of the local depth
318	hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
319	hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
320
321	END IF
322	ztmp3d(:,:,:) = 0.0
323	ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
324	& - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
325	zsum = glob_sum_full(ztmp3d)
326	IF (lwp) PRINT *, 'total volume correction 02 = ',zsum
327
328	!=============================================================================
329
330	! compute velocity
331	! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
332	ub(:,:,:)=un(:,:,:)
333	vb(:,:,:)=vn(:,:,:)
334	DO jk = 1,jpk
335	DO ii = 1,jpi
336	DO ij = 1,jpj
337	un(ii,ij,jk) = ub(ii,ij,jk)pe3u_b(ii,ij,jk)pumask_b(ii,ij,jk)/fse3u_n(ii,ij,jk)*umask(ii,ij,jk);
338	vn(ii,ij,jk) = vb(ii,ij,jk)pe3v_b(ii,ij,jk)pvmask_b(ii,ij,jk)/fse3v_n(ii,ij,jk)*vmask(ii,ij,jk);
339	!IF ( umask(ii,ij,1) == 1 .AND. pumask(ii,ij,1) == 0 ) un(ii,ij,jk) = ub(ii,ij,jk)
340	!IF ( vmask(ii,ij,1) == 1 .AND. pvmask(ii,ij,1) == 0 ) vn(ii,ij,jk) = vb(ii,ij,jk)
341	END DO
342	END DO
343	END DO
344	! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
345	! compute barotropic velocity now and after
346	ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:);
347	zbub(:,:) = SUM(ztrp,DIM=3)
348	ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:);
349	zbvb(:,:) = SUM(ztrp,DIM=3)
350	ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:);
351	zbun(:,:) = SUM(ztrp,DIM=3)
352	ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:);
353	zbvn(:,:) = SUM(ztrp,DIM=3)
354	! Already know ????????
355	hu1=0.0_wp ;
356	hv1=0.0_wp ;
357	DO jk = 1,jpk
358	hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
359	hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
360	END DO
361	zucorr = 0._wp
362	zvcorr = 0._wp
363	DO ii = 1,jpi
364	DO ij = 1,jpj
365	IF (zbun(ii,ij) .NE. zbub(ii,ij) .AND. hu1(ii,ij) .NE. 0._wp ) THEN
366	zucorr(ii,ij) = (zbun(ii,ij) - zbub(ii,ij))/hu1(ii,ij)
367	END IF
368	IF (zbvn(ii,ij) .NE. zbvb(ii,ij) .AND. hv1(ii,ij) .NE. 0._wp ) THEN
369	zvcorr(ii,ij) = (zbvn(ii,ij) - zbvb(ii,ij))/hv1(ii,ij)
370	END IF
371	END DO
372	END DO
373	DO jk = 1,jpk
374	un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
375	vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
376	END DO
377
378	!=============================================================================
379	! compute temp and salt
380	! compute new tn and sn if we open a new cell
381	tsb (:,:,:,:) = tsn(:,:,:,:)
382	zts0(:,:,:,:) = tsn(:,:,:,:)
383	ztmask1(:,:,:) = ptmask_b(:,:,:)
384	ztmask0(:,:,:) = ptmask_b(:,:,:)
385	DO iz = 1,10
386	DO jk = 1,jpk-1
387	zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
388	DO ii = 2,jpi-1
389	DO ij = 2,jpj-1
390	iip1=ii+1; iim1=ii-1;
391	ijp1=ij+1; ijm1=ij-1;
392	summsk=(ztmask0(iip1,ij,jk)+ztmask0(iim1,ij,jk)+ztmask0(ii,ijp1,jk)+ztmask0(ii,ijm1,jk))
393	!! horizontal extrapolation
394	IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
395	tsn(ii,ij,jk,1)=( zts0(iip1,ij,jk,1)*ztmask0(iip1,ij,jk) &
396	& +zts0(iim1,ij,jk,1)*ztmask0(iim1,ij,jk) &
397	& +zts0(ii,ijp1,jk,1)*ztmask0(ii,ijp1,jk) &
398	& +zts0(ii,ijm1,jk,1)*ztmask0(ii,ijm1,jk))/summsk
399	tsn(ii,ij,jk,2)=( zts0(iip1,ij,jk,2)*ztmask0(iip1,ij,jk) &
400	& +zts0(iim1,ij,jk,2)*ztmask0(iim1,ij,jk) &
401	& +zts0(ii,ijp1,jk,2)*ztmask0(ii,ijp1,jk) &
402	& +zts0(ii,ijm1,jk,2)*ztmask0(ii,ijm1,jk))/summsk
403	ztmask1(ii,ij,jk)=1
404	END IF
405	!! vertical extrapolation if horizontal extra failed
406	IF (zdmask(ii,ij)==1 .AND. summsk==0) THEN
407	jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
408	summsk=(ztmask0(ii,ij,jkm1)+ztmask0(ii,ij,jkp1))
409	IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0 ) THEN
410	tsn(ii,ij,jk,1)=( zts0(ii,ij,jkp1,1)*ztmask0(ii,ij,jkp1) &
411	& +zts0(ii,ij,jkm1,1)*ztmask0(ii,ij,jkm1))/summsk
412	tsn(ii,ij,jk,2)=( zts0(ii,ij,jkp1,2)*ztmask0(ii,ij,jkp1) &
413	& +zts0(ii,ij,jkm1,2)*ztmask0(ii,ij,jkm1))/summsk
414	ztmask1(ii,ij,jk)=1
415	END IF
416	END IF
417	!! horizontal corner extrapolation if horizontal and vertical failed
418	IF (zdmask(ii,ij)==1 .AND. summsk==0) THEN
419	iip1=ii+1; iim1=ii-1;
420	ijp1=ij+1; ijm1=ij-1;
421	summsk=(ztmask0(iip1,ijp1,jk)+ztmask0(iim1,ijm1,jk)+ztmask0(iim1,ijp1,jk)+ztmask0(iip1,ijm1,jk))
422	IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
423	tsn(ii,ij,jk,1)=( zts0(iip1,ijp1,jk,1)*ztmask0(iip1,ijp1,jk) &
424	& +zts0(iim1,ijm1,jk,1)*ztmask0(iim1,ijm1,jk) &
425	& +zts0(iim1,ijp1,jk,1)*ztmask0(iim1,ijp1,jk) &
426	& +zts0(iip1,ijm1,jk,1)*ztmask0(iip1,ijm1,jk))/summsk
427	tsn(ii,ij,jk,2)=( zts0(iip1,ijp1,jk,2)*ztmask0(iip1,ijp1,jk) &
428	& +zts0(iim1,ijm1,jk,2)*ztmask0(iim1,ijm1,jk) &
429	& +zts0(iim1,ijp1,jk,2)*ztmask0(iim1,ijp1,jk) &
430	& +zts0(iip1,ijm1,jk,2)*ztmask0(iip1,ijm1,jk))/summsk
431	ztmask1(ii,ij,jk)=1
432	END IF
433	END IF
434	END DO
435	END DO
436	END DO
437	CALL lbc_lnk(tsn(:,:,:,1),'T',1.)
438	CALL lbc_lnk(tsn(:,:,:,2),'T',1.)
439	CALL lbc_lnk(ztmask1,'T',1.)
440	zts0(:,:,:,:) = tsn(:,:,:,:)
441	ztmask0 = ztmask1
442	END DO
443	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
444	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
445
446	! CHECK heat
447	zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
448	zsumb = glob_sum( tsb(:,:,:,jp_tem) * pe3t_b(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
449	IF (lwp) PRINT *, 'CHECK tsn = ',zsumn, zsumb
450
451	! compute new T/S (interpolation) because of vvl
452	IF ( lk_vvl .AND. .FALSE. ) THEN
453	!IF ( lk_vvl ) THEN
454	DO jk = 2,jpk-1
455	DO jj = 1,jpj
456	DO ji = 1,jpi
457	IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1.0_wp ) THEN
458	!compute weight
459	zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
460	zdz = fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk )
461	zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk ) - fsdepw_n(ji,jj,jk ))
462	IF (zdz .LT. 0.0_wp) THEN
463	PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
464	PRINT *, 'ERROR dz n = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
465	PRINT *, 'ERROR dz b = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
466	PRINT *, 'ERROR dz b = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
467	PRINT *, 'ERROR dz 0 = ', e3t_0 (ji,jj,jk+1), e3t_0 (ji,jj,jk), e3t_0 (ji,jj,jk-1)
468	PRINT *, 'ERROR dz n = ', tmask (ji,jj,jk+1), tmask (ji,jj,jk), tmask (ji,jj,jk-1)
469	PRINT *, 'ERROR dz n = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
470	PRINT *, 'ERROR dz b = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
471	PRINT *, 'ERROR dz b = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
472	PRINT *, 'ERROR dz b = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
473	STOP
474	END IF
475	tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
476	& + zdz *tsb(ji,jj,jk ,jp_tem) &
477	& + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/fse3t_n(ji,jj,jk)
478	tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
479	& + zdz *tsb(ji,jj,jk ,jp_sal) &
480	& + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/fse3t_n(ji,jj,jk)
481	END IF
482	END DO
483	END DO
484	END DO
485	END IF
486
487	! CHECK heat
488	ztmp3d(:,:,:) = 0.0
489	ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
490	& - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
491	zsum = glob_sum_full(ztmp3d)
492	IF (lwp) PRINT *, 'total volume correction 03 = ',zsum
493
494	WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp)
495	tsn(:,:,:,2)= -99._wp
496	tmask(:,:,:) = 0.0
497	umask(:,:,:) = 0.0
498	vmask(:,:,:) = 0.0
499	END WHERE
500
501	WHERE (SUM(tmask,dim=3) == 0)
502	mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
503	mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
504	END WHERE
505	! compute new tn and sn if we close cell
506	! nothing to do
507
508	!=============================================================================
509	! put this stuff in a subroutine ????
510	!IF ( ln_hfb ) CALL scb_iscpl_cons(pvol_flx, pts_flx, pe3t_b, zssh0)
511	! initialisation
512	zde3t (:,:) = 0.0_wp
513	pvol_flx(:,:,: ) = 0.0_wp
514	pts_flx (:,:,:,:) = 0.0_wp
515	IF ( ln_hfb ) THEN
516	zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt
517	IF (lwp) PRINT *, 'total volume correction 0 = ',zsum
518	zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
519	IF (lwp) PRINT *, 'total heat correction 0 = ',zsum
520	zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
521	IF (lwp) PRINT *, 'total salt correction 0 = ',zsum
522
523	! mask tsn and tsb (should be useless)
524	tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)tmask(:,:,:);
525	tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)tmask(:,:,:);
526	! why tsn = tsb ????? (bug ???????)
527	! diagnose non conservation of heat, salt and volume
528	r1_tiscpl = 1._wp / (rn_fiscpl * rn_rdt)
529	zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
530	IF ( lk_vvl ) zssh0 = 0.0_wp
531	DO jk = 1,jpk-1
532	DO ii = 2,jpi-1
533	DO ij = 2,jpj-1
534	! volume differences
535	zde3t(ii,ij) = fse3t_n(ii,ij,jk) * tmask(ii,ij,jk) - pe3t_b(ii,ij,jk) * ptmask_b(ii,ij,jk);
536
537	! shh changes
538	IF ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) THEN
539	zde3t(ii,ij) = zde3t(ii,ij) + zssh0(ii,ij)
540	zssh0(ii,ij) = 0._wp
541	END IF
542
543	! ocean cell now
544	! IF (tmask(ii,ij,jk) == 1._wp) THEN
545	! case where we open, enlarge or thin a cell :
546	pvol_flx(ii,ij,jk) = zde3t(ii,ij) * r1_tiscpl
547	pts_flx (ii,ij,jk,jp_sal)= tsn(ii,ij,jk,jp_sal) * zde3t(ii,ij) * r1_tiscpl
548	pts_flx (ii,ij,jk,jp_tem)= tsn(ii,ij,jk,jp_tem) * zde3t(ii,ij) * r1_tiscpl
549	! END IF
550	END DO
551	END DO
552	END DO
553	! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0
554	PRINT , 'test ', narea, SUM(pvol_flx(:,:,:)) rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt
555	zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt
556	IF (lwp) PRINT *, 'total volume correction 1 = ',zsum
557	zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
558	IF (lwp) PRINT *, 'total heat correction 1 = ',zsum
559	zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
560	IF (lwp) PRINT *, 'total salt correction 1 = ',zsum
561
562	zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
563	IF ( lk_vvl ) zssh0 = 0.0_wp
564	DO jk = 1,jpk-1
565	DO ii = 2,jpi-1
566	DO ij = 2,jpj-1
567	! volume differences
568	zde3t(ii,ij) = fse3t_n(ii,ij,jk) * tmask(ii,ij,jk) - pe3t_b(ii,ij,jk) * ptmask_b(ii,ij,jk);
569
570	! shh changes
571	!IF ( ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) .AND. zssh0(ii,ij) .NE. 0._wp) THEN
572	IF ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) THEN
573	zde3t(ii,ij) = zde3t(ii,ij) + zssh0(ii,ij)
574	zssh0(ii,ij) = 0._wp
575	END IF
576
577	! ocean cell before and mask cell now
578	IF ( tmask(ii,ij,jk) == 0._wp .AND. ptmask_b(ii,ij,jk) == 1._wp ) THEN
579	! case where we close a cell and adjacent cell open
580	pvol_flx(ii,ij,jk) = zde3t(ii,ij) * r1_tiscpl
581	pts_flx (ii,ij,jk,jp_sal)= tsb(ii,ij,jk,jp_sal) * zde3t(ii,ij) * r1_tiscpl
582	pts_flx (ii,ij,jk,jp_tem)= tsb(ii,ij,jk,jp_tem) * zde3t(ii,ij) * r1_tiscpl
583
584	iip1=ii+1 ; iim1=ii-1 ; ijp1=ij+1 ; ijm1=ij-1 ;
585
586	zsum = e12t(ii ,ijp1) * tmask(ii ,ijp1,jk) + e12t(ii ,ijm1) * tmask(ii ,ijm1,jk) &
587	& + e12t(iim1,ij ) * tmask(iim1,ij ,jk) + e12t(iip1,ij ) * tmask(iip1,ij ,jk)
588
589	IF ( zsum .NE. 0._wp ) THEN
590	ziip1_ratio = e12t(iip1,ij ) * tmask(iip1,ij ,jk) / zsum
591	ziim1_ratio = e12t(iim1,ij ) * tmask(iim1,ij ,jk) / zsum
592	zijp1_ratio = e12t(ii ,ijp1) * tmask(ii ,ijp1,jk) / zsum
593	zijm1_ratio = e12t(ii ,ijm1) * tmask(ii ,ijm1,jk) / zsum
594
595	pvol_flx(ii ,ijp1,jk ) = pvol_flx(ii ,ijp1,jk ) + pvol_flx(ii,ij,jk ) * zijp1_ratio
596	pvol_flx(ii ,ijm1,jk ) = pvol_flx(ii ,ijm1,jk ) + pvol_flx(ii,ij,jk ) * zijm1_ratio
597	pvol_flx(iip1,ij ,jk ) = pvol_flx(iip1,ij ,jk ) + pvol_flx(ii,ij,jk ) * ziip1_ratio
598	pvol_flx(iim1,ij ,jk ) = pvol_flx(iim1,ij ,jk ) + pvol_flx(ii,ij,jk ) * ziim1_ratio
599	pts_flx (ii ,ijp1,jk,jp_sal) = pts_flx (ii ,ijp1,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * zijp1_ratio
600	pts_flx (ii ,ijm1,jk,jp_sal) = pts_flx (ii ,ijm1,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * zijm1_ratio
601	pts_flx (iip1,ij ,jk,jp_sal) = pts_flx (iip1,ij ,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * ziip1_ratio
602	pts_flx (iim1,ij ,jk,jp_sal) = pts_flx (iim1,ij ,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * ziim1_ratio
603	pts_flx (ii ,ijp1,jk,jp_tem) = pts_flx (ii ,ijp1,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * zijp1_ratio
604	pts_flx (ii ,ijm1,jk,jp_tem) = pts_flx (ii ,ijm1,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * zijm1_ratio
605	pts_flx (iip1,ij ,jk,jp_tem) = pts_flx (iip1,ij ,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * ziip1_ratio
606	pts_flx (iim1,ij ,jk,jp_tem) = pts_flx (iim1,ij ,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * ziim1_ratio
607
608	! set to 0 the cell we distributed over neigbourg cells
609	pvol_flx(ii,ij,jk ) = 0._wp
610	pts_flx (ii,ij,jk,jp_sal) = 0._wp
611	pts_flx (ii,ij,jk,jp_tem) = 0._wp
612	ELSE IF (zsum .EQ. 0._wp ) THEN
613	! case where we close a cell and no adjacent cell open
614	IF ( tmask(ii,ij,jk+1) .EQ. 1._wp ) THEN
615	pvol_flx(ii,ij,jk+1) = pvol_flx(ii,ij,jk+1) + pvol_flx(ii,ij,jk)
616	pts_flx (ii,ij,jk+1,jp_sal)= pts_flx (ii,ij,jk+1,jp_sal) + pts_flx (ii,ij,jk,jp_sal)
617	pts_flx (ii,ij,jk+1,jp_tem)= pts_flx (ii,ij,jk+1,jp_tem) + pts_flx (ii,ij,jk,jp_tem)
618
619	! set to 0 the cell we distributed over neigbourg cells
620	pvol_flx(ii,ij,jk ) = 0._wp
621	pts_flx (ii,ij,jk,jp_sal) = 0._wp
622	pts_flx (ii,ij,jk,jp_tem) = 0._wp
623	ELSE
624	! case no adjacent cell on the horizontal and on the vertical
625	PRINT *, 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
626	PRINT *, ' ',mig(ii),' ',mjg(ij),' ',jk
627	PRINT *, ' ',ii,' ',ij,' ',jk,' ',narea
628	PRINT *, ' conservation broken '
629	! We deal with this point later.
630	END IF
631	END IF
632	END IF
633	END DO
634	END DO
635	END DO
636
637	zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt
638	IF (lwp) PRINT *, 'total volume correction 2 = ',zsum
639	zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
640	IF (lwp) PRINT *, 'total heat correction 2 = ',zsum
641	zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
642	IF (lwp) PRINT *, 'total salt correction 2 = ',zsum
643
644	! allocation and initialisation of the list of problematic point
645	ALLOCATE(vnpts(jpnij))
646	vnpts(:)=0
647
648	! fill narea location with the number of problematic point
649	DO jk = 1,jpk-1
650	DO ii = 2,jpi-1
651	DO ij = 2,jpj-1
652	IF ( ptmask_b(ii,ij,jk) == 1 .AND. SUM(tmask(ii-1:ii+1,ij,jk)) == 0 &
653	& .AND. SUM(tmask(ii,ij-1:ij+1,jk)) == 0 .AND. tmask(ii,ij,jk+1) == 0 ) THEN
654	vnpts(narea) = vnpts(narea) + 1
655	END IF
656	END DO
657	END DO
658	END DO
659
660	! build array of total problematic point on each cpu (share to each cpu)
661	CALL mpp_max(vnpts,jpnij)
662
663	! size of the new variable
664	npts = SUM(vnpts)
665
666	! allocation of the coordiantes, correction, index vector for the problematic points
667	ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
668	ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20
669	zcorr_vol(:) = 0.0_wp
670	zcorr_sal(:) = 0.0_wp
671	zcorr_tem(:) = 0.0_wp
672
673	! fill new variable
674	jpts = SUM(vnpts(1:narea-1))
675	DO jk = 1,jpk-1
676	DO ii = 2,jpi-1
677	DO ij = 2,jpj-1
678	IF ( ptmask_b(ii,ij,jk) == 1 .AND. SUM(tmask(ii-1:ii+1,ij,jk)) == 0 &
679	& .AND. SUM(tmask(ii,ij-1:ij+1,jk)) == 0 .AND. tmask(ii,ij,jk+1) == 0 ) THEN
680	jpts = jpts + 1
681	PRINT *, 'corrected point ', narea, ii, ij, jk, jpts
682	ixpts(jpts) = ii ; iypts(jpts) = ij ; izpts(jpts) = jk
683	zlon (jpts) = glamt(ii,ij) ; zlat (jpts) = gphit(ii,ij)
684	zcorr_vol(jpts) = pvol_flx(ii,ij,jk)
685	zcorr_sal(jpts) = pts_flx (ii,ij,jk,jp_sal)
686	zcorr_tem(jpts) = pts_flx (ii,ij,jk,jp_tem)
687	! set flx to 0 (safer)
688	pvol_flx(ii,ij,jk ) = 0.0_wp
689	pts_flx (ii,ij,jk,jp_sal) = 0.0_wp
690	pts_flx (ii,ij,jk,jp_tem) = 0.0_wp
691	PRINT , zcorr_vol(jpts)rn_fiscplrn_rdt, zcorr_sal(jpts)rn_fiscplrn_rdt, zcorr_tem(jpts)rn_fiscpl*rn_rdt
692	END IF
693	END DO
694	END DO
695	END DO
696
697	! build array of total problematic point on each cpu (share to each cpu)
698	CALL mpp_max(zlat ,npts)
699	CALL mpp_max(zlon ,npts)
700	CALL mpp_max(izpts,npts)
701
702	! put correction term in the closest cell
703	PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts
704	DO jpts = 1,npts
705	CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
706	PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts)
707	DO ij = mj0(iypts(jpts)),mj1(iypts(jpts))
708	DO ii = mi0(ixpts(jpts)),mi1(ixpts(jpts))
709	jk = izpts(jpts)
710	pvol_flx(ii,ij,jk) = pvol_flx(ii,ij,jk ) + zcorr_vol(jpts)
711	pts_flx (ii,ij,jk,jp_sal) = pts_flx (ii,ij,jk,jp_sal) + zcorr_sal(jpts)
712	pts_flx (ii,ij,jk,jp_tem) = pts_flx (ii,ij,jk,jp_tem) + zcorr_tem(jpts)
713	END DO
714	END DO
715	END DO
716	! deallocate variable
717	DEALLOCATE(vnpts)
718	DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
719
720	! add contribution store on the allo (lbclnk remove one of the contribution)
721	pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:)
722	pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
723	pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
724
725	CALL lbc_sum(pvol_flx(:,:,: ),'T',1.)
726	CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
727	CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
728
729	! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!!
730	zsumn = glob_sum ( fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt
731	ztmp3d(:,:,:) = 0.0
732	ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
733	zsumb = glob_sum_full(ztmp3d)
734	zsum = glob_sum ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt)
735	IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum
736	! CHECK salt
737	zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
738	zsumb = glob_sum( tsb(:,:,:,jp_sal) * pe3t_b(:,:,:) * ptmask_b(:,:,:))
739	zsum = glob_sum( pts_flx(:,:,:,jp_sal)rn_fiscpl rn_rdt)
740	IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum
741	! CHECK heat
742	zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
743	zsumb = glob_sum( tsb(:,:,:,jp_tem) * pe3t_b(:,:,:) * ptmask_b(:,:,:))
744	zsum = glob_sum( pts_flx(:,:,:,jp_tem)rn_fiscpl rn_rdt)
745	IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum
746
747	END IF
748	!=============================================================================
749	!! now T/S/U are compute in agreement with the new grid.
750	!! need to interpolate on vvl grid
751	!! next step is an euler time step
752	neuler = 0
753	!! set _b and _n variables equal
754	tsb (:,:,:,:)=tsn (:,:,:,:)
755	ub (:,:,: )=un (:,:,: )
756	vb (:,:,: )=vn (:,:,: )
757	sshb(:,: )=sshn(:,:)
758	!! set _b and _n vertical scale factor equal
759	fse3t_b (:,:,:)=fse3t_n (:,:,:)
760	fse3u_b (:,:,:)=fse3u_n (:,:,:)
761	fse3v_b (:,:,:)=fse3v_n (:,:,:)
762	IF ( lk_vvl ) THEN
763	fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
764	fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
765	!! set depth_b = depth_n
766	fsdept_b(:,:,:) = fsdept_n(:,:,:)
767	fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
768	!! set hu_b
769	hu_b(:,:) = hu(:,:)
770	hv_b(:,:) = hv(:,:)
771	hur_b(:,:) = hur(:,:)
772	hvr_b(:,:) = hvr(:,:)
773	END IF
774	!!
775	CALL wrk_dealloc(jpi,jpj,jpk,2, zts0 )
776	CALL wrk_dealloc(jpi,jpj,jpk, ztmask0, ztmask1 , ztrp )
777	CALL wrk_dealloc(jpi,jpj,jpk, zwmaskn, zwmaskb , ztmp3d )
778	CALL wrk_dealloc(jpi,jpj, zsmask0, zsmask1 )
779	CALL wrk_dealloc(jpi,jpj, zdmask , zdsmask, zvcorr, zucorr, zde3t)
780	CALL wrk_dealloc(jpi,jpj, zbub , zbvb , zbun , zbvn )
781	CALL wrk_dealloc(jpi,jpj, zssh0 , zssh1 , hu1 , hv1 )
782	END SUBROUTINE rst_iscpl_interpol
783
784	SUBROUTINE sbc_iscpl_div( phdivn )
785	!!----------------------------------------------------------------------
786	!! * ROUTINE sbc_iscpl_div *
787	!!
788	!! ** Purpose : update the horizontal divergence with the iscpl imbalance
789	!!
790	!! ** Method :
791	!! CAUTION : iscpl is positive (inflow) increase the
792	!! divergence and expressed in m/s
793	!!
794	!! ** Action : phdivn decreased by the runoff inflow
795	!!----------------------------------------------------------------------
796	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
797	!!
798	INTEGER :: ji, jj, jk ! dummy loop indices
799	!!----------------------------------------------------------------------
800	!
801	DO jk = 1, jpk
802	DO jj = 1, jpj
803	DO ji = 1, jpi
804	phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
805	END DO
806	END DO
807	END DO
808	!
809	END SUBROUTINE sbc_iscpl_div
810
811	END MODULE sbc_iscpl

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

Download in other formats: