Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

iscplhsb.F90 @ 8314

Last change on this file since 8314 was 8186, checked in by acc, 7 years ago

Branch 2017/dev_r8126_ROBUST08_no_ghost. Incorporation of re-written lbc routines. This introduces generic routines for: lbc_lnk, lbc_lnk_multi, lbc_nfd, mpp_bdy, mpp_lnk and mpp_nfd in .h90 files which are pre-processor included multiple times (with different arguments) to recreate equivalences to all the original variants from a much smaller code base (more than 2000 lines shorter). These changes have been SETTE tested and shown to reproduce identical results to the branch base revision. There are a few caveats: the ice cavity routine: iscplhsb.F90, needs to be rewritten to avoid sums over the overlap regions; this will be done elsewhere and has merely been disabled on this branch. The work is not yet complete for the nogather option for the north-fold. The default MPI ALLGATHER option is working but do not activate ln_nogather until further notice.

File size: 15.8 KB

Line
1	MODULE iscplhsb
2	!!======================================================================
3	!! * MODULE iscplhsb *
4	!! Ocean forcing: ice sheet/ocean coupling (conservation)
5	!!=====================================================================
6	!! History : NEMO ! 2015-01 P. Mathiot: original
7	!!----------------------------------------------------------------------
8
9	!!----------------------------------------------------------------------
10	!! iscpl_alloc : variable allocation
11	!! iscpl_hsb : compute and store the input of heat/salt/volume
12	!! into the system due to the coupling process
13	!! iscpl_div : correction of divergence to keep volume conservation
14	!!----------------------------------------------------------------------
15	USE dom_oce ! ocean space and time domain
16	USE domwri ! ocean space and time domain
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 lib_mpp ! MPP library
22	USE lib_fortran ! MPP library
23	USE wrk_nemo ! Memory allocation
24	USE lbclnk !
25	USE domngb !
26	USE iscplini
27
28	IMPLICIT NONE
29	PRIVATE
30
31	PUBLIC iscpl_div
32	PUBLIC iscpl_cons
33	!! * Substitutions
34	# include "vectopt_loop_substitute.h90"
35	!!----------------------------------------------------------------------
36	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
37	!! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
39	!!----------------------------------------------------------------------
40	CONTAINS
41
42	SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
43	!!----------------------------------------------------------------------
44	!! * ROUTINE iscpl_cons *
45	!!
46	!! ** Purpose : compute input into the system during the coupling step
47	!! compute the correction term
48	!! compute where the correction have to be applied
49	!!
50	!! ** Method : compute tsne3t-tsbe3tb and e3t-e3t_b
51	!!----------------------------------------------------------------------
52	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b !! mask before
53	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b !! scale factor before
54	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: psmask_b !! mask before
55	REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx !! corrective flux to have tracer conservation
56	REAL(wp), DIMENSION(:,:,: ), INTENT(out) :: pvol_flx !! corrective flux to have volume conservation
57	REAL(wp), INTENT(in ) :: prdt_iscpl !! coupling period
58	!!
59	INTEGER :: ji, jj, jk !! loop index
60	INTEGER :: jip1, jim1, jjp1, jjm1
61	!!
62	REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
63	REAL(wp):: r1_rdtiscpl
64	REAL(wp):: zjip1_ratio , zjim1_ratio , zjjp1_ratio , zjjm1_ratio
65	!!
66	REAL(wp):: zde3t, zdtem, zdsal
67	REAL(wp), DIMENSION(:,:), POINTER :: zdssh
68	!!
69	REAL(wp), DIMENSION(:), ALLOCATABLE :: zlon, zlat
70	REAL(wp), DIMENSION(:), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
71	INTEGER , DIMENSION(:), ALLOCATABLE :: ixpts, iypts, izpts, inpts
72	INTEGER :: jpts, npts
73
74	CALL wrk_alloc(jpi,jpj, zdssh )
75
76	! get imbalance (volume heat and salt)
77	! initialisation difference
78	zde3t = 0.0_wp; zdsal = 0.0_wp ; zdtem = 0.0_wp
79
80	! initialisation correction term
81	pvol_flx(:,:,: ) = 0.0_wp
82	pts_flx (:,:,:,:) = 0.0_wp
83
84	r1_rdtiscpl = 1._wp / prdt_iscpl
85
86	! mask tsn and tsb
87	tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)tmask(:,:,:);
88	tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)tmask(:,:,:);
89
90	!==============================================================================
91	! diagnose the heat, salt and volume input and compute the correction variable
92	!==============================================================================
93
94	!
95	zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
96	IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
97
98	DO jk = 1,jpk-1
99	DO jj = 2,jpj-1
100	DO ji = fs_2,fs_jpim1
101	IF (tmask_h(ji,jj) == 1._wp) THEN
102
103	! volume differences
104	zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
105
106	! heat diff
107	zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) * tmask (ji,jj,jk) &
108	- tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
109	! salt diff
110	zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) * tmask (ji,jj,jk) &
111	- tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
112
113	! shh changes
114	IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN
115	zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
116	zdssh(ji,jj) = 0._wp
117	END IF
118
119	! volume, heat and salt differences in each cell
120	pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk) + zde3t * r1_rdtiscpl
121	pts_flx (ji,jj,jk,jp_sal)= pts_flx (ji,jj,jk,jp_sal) + zdsal * r1_rdtiscpl
122	pts_flx (ji,jj,jk,jp_tem)= pts_flx (ji,jj,jk,jp_tem) + zdtem * r1_rdtiscpl
123
124	! case where we close a cell: check if the neighbour cells are wet
125	IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
126
127	jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
128
129	zsum = e1e2t(ji ,jjp1) * tmask(ji ,jjp1,jk) + e1e2t(ji ,jjm1) * tmask(ji ,jjm1,jk) &
130	& + e1e2t(jim1,jj ) * tmask(jim1,jj ,jk) + e1e2t(jip1,jj ) * tmask(jip1,jj ,jk)
131
132	IF ( zsum /= 0._wp ) THEN
133	zjip1_ratio = e1e2t(jip1,jj ) * tmask(jip1,jj ,jk) / zsum
134	zjim1_ratio = e1e2t(jim1,jj ) * tmask(jim1,jj ,jk) / zsum
135	zjjp1_ratio = e1e2t(ji ,jjp1) * tmask(ji ,jjp1,jk) / zsum
136	zjjm1_ratio = e1e2t(ji ,jjm1) * tmask(ji ,jjm1,jk) / zsum
137
138	pvol_flx(ji ,jjp1,jk ) = pvol_flx(ji ,jjp1,jk ) + pvol_flx(ji,jj,jk ) * zjjp1_ratio
139	pvol_flx(ji ,jjm1,jk ) = pvol_flx(ji ,jjm1,jk ) + pvol_flx(ji,jj,jk ) * zjjm1_ratio
140	pvol_flx(jip1,jj ,jk ) = pvol_flx(jip1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjip1_ratio
141	pvol_flx(jim1,jj ,jk ) = pvol_flx(jim1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjim1_ratio
142	pts_flx (ji ,jjp1,jk,jp_sal) = pts_flx (ji ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
143	pts_flx (ji ,jjm1,jk,jp_sal) = pts_flx (ji ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
144	pts_flx (jip1,jj ,jk,jp_sal) = pts_flx (jip1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
145	pts_flx (jim1,jj ,jk,jp_sal) = pts_flx (jim1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
146	pts_flx (ji ,jjp1,jk,jp_tem) = pts_flx (ji ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
147	pts_flx (ji ,jjm1,jk,jp_tem) = pts_flx (ji ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
148	pts_flx (jip1,jj ,jk,jp_tem) = pts_flx (jip1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
149	pts_flx (jim1,jj ,jk,jp_tem) = pts_flx (jim1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
150
151	! set to 0 the cell we distributed over neigbourg cells
152	pvol_flx(ji,jj,jk ) = 0._wp
153	pts_flx (ji,jj,jk,jp_sal) = 0._wp
154	pts_flx (ji,jj,jk,jp_tem) = 0._wp
155
156	ELSE IF (zsum == 0._wp ) THEN
157	! case where we close a cell and no adjacent cell open
158	! check if the cell beneath is wet
159	IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
160	pvol_flx(ji,jj,jk+1) = pvol_flx(ji,jj,jk+1) + pvol_flx(ji,jj,jk)
161	pts_flx (ji,jj,jk+1,jp_sal)= pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
162	pts_flx (ji,jj,jk+1,jp_tem)= pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
163
164	! set to 0 the cell we distributed over neigbourg cells
165	pvol_flx(ji,jj,jk ) = 0._wp
166	pts_flx (ji,jj,jk,jp_sal) = 0._wp
167	pts_flx (ji,jj,jk,jp_tem) = 0._wp
168	ELSE
169	! case no adjacent cell on the horizontal and on the vertical
170	IF ( lwp ) THEN ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
171	! open for narea== 1 (lwp=T)
172	WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
173	WRITE(numout,*) ' ',mig(ji),' ',mjg(jj),' ',jk
174	WRITE(numout,*) ' ',ji,' ',jj,' ',jk,' ',narea
175	WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
176	ENDIF
177	! We deal with these points later.
178	END IF
179	END IF
180	END IF
181	END IF
182	END DO
183	END DO
184	END DO
185
186	!!gm ERROR !!!!
187	!! juste use tmask_i or in case of ISF smask_i (to be created to compute the sum without halos)
188	!
189	! CALL lbc_sum(pvol_flx(:,:,: ),'T',1.)
190	! CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
191	! CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
192	STOP ' iscpl_cons: please modify this module !'
193	!!gm end
194	! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point
195	! allocation and initialisation of the list of problematic point
196	ALLOCATE(inpts(jpnij))
197	inpts(:)=0
198
199	! fill narea location with the number of problematic point
200	DO jk = 1,jpk-1
201	DO jj = 2,jpj-1
202	DO ji = fs_2,fs_jpim1
203	IF ( ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1) == 0._wp .AND. tmask_h(ji,jj) == 1._wp &
204	.AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
205	inpts(narea) = inpts(narea) + 1
206	END IF
207	END DO
208	END DO
209	END DO
210
211	! build array of total problematic point on each cpu (share to each cpu)
212	CALL mpp_max(inpts,jpnij)
213
214	! size of the new variable
215	npts = SUM(inpts)
216
217	! allocation of the coordinates, correction, index vector for the problematic points
218	ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
219	ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
220	zcorr_vol(:) = -1.0e20_wp
221	zcorr_sal(:) = -1.0e20_wp
222	zcorr_tem(:) = -1.0e20_wp
223
224	! fill new variable
225	jpts = SUM(inpts(1:narea-1))
226	DO jk = 1,jpk-1
227	DO jj = 2,jpj-1
228	DO ji = fs_2,fs_jpim1
229	IF ( ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1) == 0._wp .AND. tmask_h(ji,jj) == 1._wp &
230	.AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
231	jpts = jpts + 1 ! positioning in the inpts vector for the area narea
232	ixpts(jpts) = ji ; iypts(jpts) = jj ; izpts(jpts) = jk
233	zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
234	zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
235	zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
236	zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
237
238	! set flx to 0 (safer)
239	pvol_flx(ji,jj,jk ) = 0.0_wp
240	pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
241	pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
242	END IF
243	END DO
244	END DO
245	END DO
246
247	! build array of total problematic point on each cpu (share to each cpu)
248	! point coordinates
249	CALL mpp_max(zlat ,npts)
250	CALL mpp_max(zlon ,npts)
251	CALL mpp_max(izpts,npts)
252
253	! correction values
254	CALL mpp_max(zcorr_vol,npts)
255	CALL mpp_max(zcorr_sal,npts)
256	CALL mpp_max(zcorr_tem,npts)
257
258	! put correction term in the closest cell
259	DO jpts = 1,npts
260	CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
261	DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
262	DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
263	jk = izpts(jpts)
264
265	IF (tmask_h(ji,jj) == 1._wp) THEN
266	! correct the vol_flx in the closest cell
267	pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk ) + zcorr_vol(jpts)
268	pts_flx (ji,jj,jk,jp_sal) = pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
269	pts_flx (ji,jj,jk,jp_tem) = pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
270
271	! set correction to 0
272	zcorr_vol(jpts) = 0.0_wp
273	zcorr_sal(jpts) = 0.0_wp
274	zcorr_tem(jpts) = 0.0_wp
275	END IF
276	END DO
277	END DO
278	END DO
279
280	! deallocate variables
281	DEALLOCATE(inpts)
282	DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
283
284	! add contribution store on the hallo (lbclnk remove one of the contribution)
285	pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:)
286	pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
287	pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
288
289	!!gm ERROR !!!!
290	!! juste use tmask_i or in case of ISF smask_i (to be created to compute the sum without halos)
291	!
292	! ! compute sum over the halo and set it to 0.
293	! CALL lbc_sum(pvol_flx(:,:,: ),'T',1._wp)
294	! CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
295	! CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
296	!!gm end
297
298	! deallocate variables
299	CALL wrk_dealloc(jpi,jpj, zdssh )
300
301	END SUBROUTINE iscpl_cons
302
303	SUBROUTINE iscpl_div( phdivn )
304	!!----------------------------------------------------------------------
305	!! * ROUTINE iscpl_div *
306	!!
307	!! ** Purpose : update the horizontal divergenc
308	!!
309	!! ** Method :
310	!! CAUTION : iscpl is positive (inflow) and expressed in m/s
311	!!
312	!! ** Action : phdivn increase by the iscpl correction term
313	!!----------------------------------------------------------------------
314	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
315	!!
316	INTEGER :: ji, jj, jk ! dummy loop indices
317	!!----------------------------------------------------------------------
318	!
319	DO jk = 1, jpk
320	DO jj = 1, jpj
321	DO ji = 1, jpi
322	phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
323	END DO
324	END DO
325	END DO
326	!
327	END SUBROUTINE iscpl_div
328
329	END MODULE iscplhsb

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90 @ 8314

Download in other formats: