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.

agrif_opa_update.F90 in branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 7824

Last change on this file since 7824 was 7824, checked in by timgraham, 7 years ago
Commit changes to interp & update routines for U, V, T&S. This version now works correctly for simple GYRE test case with 41 levels on the child grid.
Property svn:keywords set to `Id`
File size: 31.0 KB

Line
1	#define TWO_WAY /* TWO WAY NESTING */
2	#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
3
4	MODULE agrif_opa_update
5	#if defined key_agrif && ! defined key_offline
6	USE par_oce
7	USE oce
8	USE dom_oce
9	USE agrif_oce
10	USE in_out_manager ! I/O manager
11	USE lib_mpp
12	USE wrk_nemo
13	USE zdf_oce ! vertical physics: ocean variables
14
15	IMPLICIT NONE
16	PRIVATE
17
18	PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales
19
20	! VERTICAL REFINEMENT BEGIN
21	PUBLIC Agrif_Init_UpdateScales
22	REAL, DIMENSION(:,:,:), ALLOCATABLE :: update_scales_t, update_scales_u, update_scales_v
23	!$AGRIF_DO_NOT_TREAT
24	LOGICAL :: scaleT, scaleU, scaleV = .FALSE.
25	!$AGRIF_END_DO_NOT_TREAT
26	! VERTICAL REFINEMENT END
27
28	# if defined key_zdftke
29	PUBLIC Agrif_Update_Tke
30	# endif
31	!!----------------------------------------------------------------------
32	!! NEMO/NST 3.6 , NEMO Consortium (2010)
33	!! $Id$
34	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
35	!!----------------------------------------------------------------------
36
37	CONTAINS
38
39	! VERTICAL REFINEMENT BEGIN
40	SUBROUTINE Agrif_Init_UpdateScales
41
42	scaleT = .TRUE.
43	Agrif_UseSpecialValueInUpdate = .FALSE.
44	Agrif_SpecialValueFineGrid = 0.
45	Call Agrif_Update_Variable(scales_t_id,procname=updatescales)
46	Agrif_UseSpecialValueInUpdate = .FALSE.
47	scaleT = .FALSE.
48
49	scaleU = .TRUE.
50	Call Agrif_Update_Variable(scales_u_id,procname=updatescales)
51	scaleU = .FALSE.
52
53	scaleV = .TRUE.
54	Call Agrif_Update_Variable(scales_v_id,procname=updatescales)
55	scaleV = .FALSE.
56
57	END SUBROUTINE Agrif_Init_UpdateScales
58
59	SUBROUTINE updatescales( ptab, i1, i2, j1, j2, k1, k2, before )
60
61	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
62	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
63	LOGICAL, iNTENT(in) :: before
64
65	INTEGER :: ji,jj,jk
66
67	IF (before) THEN
68	IF (scaleT) THEN
69	DO jk=k1,k2
70	DO jj=j1,j2
71	DO ji=i1,i2
72	ptab(ji,jj,jk) = e3t_n(ji,jj,jk)*tmask(ji,jj,jk)
73	END DO
74	END DO
75	END DO
76	ELSEIF (scaleU) THEN
77	DO jk=k1,k2
78	DO jj=j1,j2
79	DO ji=i1,i2
80	ptab(ji,jj,jk) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
81	END DO
82	END DO
83	END DO
84	ELSEIF (scaleV) THEN
85	DO jk=k1,k2
86	DO jj=j1,j2
87	DO ji=i1,i2
88	ptab(ji,jj,jk) = e3v_n(ji,jj,jk)*vmask(ji,jj,jk)
89	END DO
90	END DO
91	END DO
92	ENDIF
93	ELSE
94	IF (scaleT) THEN
95	Allocate(update_scales_t(i1:i2,j1:j2,k1:k2))
96	update_scales_t(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
97	ELSEIF (scaleU) THEN
98	Allocate(update_scales_u(i1:i2,j1:j2,k1:k2))
99	update_scales_u(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
100	ELSEIF (scaleV) THEN
101	Allocate(update_scales_v(i1:i2,j1:j2,k1:k2))
102	update_scales_v(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
103	ENDIF
104	ENDIF
105
106	END SUBROUTINE updatescales
107	! VERTICAL REFINEMENT END
108
109	RECURSIVE SUBROUTINE Agrif_Update_Tra( )
110	!!---------------------------------------------
111	!! * ROUTINE Agrif_Update_Tra *
112	!!---------------------------------------------
113	!
114	IF (Agrif_Root()) RETURN
115	!
116	#if defined TWO_WAY
117	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update tracers from grid Number',Agrif_Fixed(), 'nbcline', nbcline
118
119	Agrif_UseSpecialValueInUpdate = .TRUE.
120	Agrif_SpecialValueFineGrid = 0.
121	!
122	IF (MOD(nbcline,nbclineupdate) == 0) THEN
123	# if ! defined DECAL_FEEDBACK
124	CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
125	# else
126	CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
127	# endif
128	ELSE
129	# if ! defined DECAL_FEEDBACK
130	CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
131	# else
132	CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
133	# endif
134	ENDIF
135	!
136	Agrif_UseSpecialValueInUpdate = .FALSE.
137	!
138	IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
139	CALL Agrif_ChildGrid_To_ParentGrid()
140	CALL Agrif_Update_Tra()
141	CALL Agrif_ParentGrid_To_ChildGrid()
142	ENDIF
143	!
144	#endif
145	!
146	END SUBROUTINE Agrif_Update_Tra
147
148
149	RECURSIVE SUBROUTINE Agrif_Update_Dyn( )
150	!!---------------------------------------------
151	!! * ROUTINE Agrif_Update_Dyn *
152	!!---------------------------------------------
153	!
154	IF (Agrif_Root()) RETURN
155	!
156	#if defined TWO_WAY
157	IF (lwp.AND.lk_agrif_debug) Write(,) 'Update momentum from grid Number',Agrif_Fixed(), 'nbcline', nbcline
158
159	Agrif_UseSpecialValueInUpdate = .FALSE.
160	Agrif_SpecialValueFineGrid = 0.
161	!
162	IF (mod(nbcline,nbclineupdate) == 0) THEN
163	# if ! defined DECAL_FEEDBACK
164	CALL Agrif_Update_Variable(un_update_id,procname = updateU)
165	CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
166	# else
167	CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
168	CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
169	# endif
170	ELSE
171	# if ! defined DECAL_FEEDBACK
172	CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
173	CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
174	# else
175	CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
176	CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
177	# endif
178	ENDIF
179
180	# if ! defined DECAL_FEEDBACK
181	CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
182	CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)
183	# else
184	CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
185	CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)
186	# endif
187
188	IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
189	! Update time integrated transports
190	IF (mod(nbcline,nbclineupdate) == 0) THEN
191	# if ! defined DECAL_FEEDBACK
192	CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
193	CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
194	# else
195	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
196	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
197	# endif
198	ELSE
199	# if ! defined DECAL_FEEDBACK
200	CALL Agrif_Update_Variable(ub2b_update_id,locupdate=(/0,1/),procname = updateub2b)
201	CALL Agrif_Update_Variable(vb2b_update_id,locupdate=(/0,1/),procname = updatevb2b)
202	# else
203	CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateub2b)
204	CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updatevb2b)
205	# endif
206	ENDIF
207	END IF
208	!
209	nbcline = nbcline + 1
210	!
211	Agrif_UseSpecialValueInUpdate = .TRUE.
212	Agrif_SpecialValueFineGrid = 0.
213	# if ! defined DECAL_FEEDBACK
214	CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
215	# else
216	CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
217	# endif
218	Agrif_UseSpecialValueInUpdate = .FALSE.
219	!
220	#endif
221	!
222	! Do recursive update:
223	IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
224	CALL Agrif_ChildGrid_To_ParentGrid()
225	CALL Agrif_Update_Dyn()
226	CALL Agrif_ParentGrid_To_ChildGrid()
227	ENDIF
228	!
229	END SUBROUTINE Agrif_Update_Dyn
230
231	# if defined key_zdftke
232
233	SUBROUTINE Agrif_Update_Tke( kt )
234	!!---------------------------------------------
235	!! * ROUTINE Agrif_Update_Tke *
236	!!---------------------------------------------
237	!!
238	INTEGER, INTENT(in) :: kt
239	!
240	return
241
242	IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
243	# if defined TWO_WAY
244
245	Agrif_UseSpecialValueInUpdate = .TRUE.
246	Agrif_SpecialValueFineGrid = 0.
247
248	CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN )
249	CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
250	CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
251
252	Agrif_UseSpecialValueInUpdate = .FALSE.
253
254	# endif
255
256	END SUBROUTINE Agrif_Update_Tke
257
258	# endif /* key_zdftke */
259
260	SUBROUTINE updateTS( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
261	!!---------------------------------------------
262	!! * ROUTINE updateT *
263	!!---------------------------------------------
264	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
265	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
266	LOGICAL, INTENT(in) :: before
267	!!
268	INTEGER :: ji,jj,jk,jn
269	! VERTICAL REFINEMENT BEGIN
270	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
271	REAL(wp) :: h_in(k1:k2)
272	REAL(wp) :: h_out(1:jpk)
273	INTEGER :: N_in, N_out
274	REAL(wp) :: h_diff
275	REAL(wp) :: zrho_xy
276	REAL(wp) :: tabin(k1:k2,n1:n2)
277	! VERTICAL REFINEMENT END
278	!!---------------------------------------------
279	!
280	IF (before) THEN
281	zrho_xy = Agrif_rhox() * Agrif_rhoy()
282	DO jn = n1,n2-1
283	DO jk=k1,k2
284	DO jj=j1,j2
285	DO ji=i1,i2
286	ptab(ji,jj,jk,jn) = zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
287	END DO
288	END DO
289	END DO
290	END DO
291	DO jk=k1,k2
292	DO jj=j1,j2
293	DO ji=i1,i2
294	ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * zrho_xy * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
295	END DO
296	END DO
297	END DO
298
299	ELSE
300	! VERTICAL REFINEMENT BEGIN
301	ptab_child(:,:,:,:) = 0.
302
303	DO jj=j1,j2
304	DO ji=i1,i2
305	N_in = 0
306	DO jk=k1,k2 !k2 = jpk of child grid
307	IF (ptab(ji,jj,jk,n2) == 0) EXIT
308	N_in = N_in + 1
309	tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2)
310	h_in(N_in) = ptab(ji,jj,jk,n2)/e1e2t(ji,jj)
311	ENDDO
312	N_out = 0
313	DO jk=1,jpk ! jpk of parent grid
314	IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF
315	N_out = N_out + 1
316	h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
317	ENDDO
318	IF (N_in > 0) THEN
319	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
320	! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
321	! IF abs(h_diff > 1.e-8) THEN
322	! N_in = N_in + 1
323	! h_in(N_in) = h_diff
324	! tabin(N_in,:) = tabin(N_in-1,:)
325	IF (h_diff < 0) THEN
326	print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
327	print *,'Nval = ',N_out,mbathy(ji,jj)
328	print *,'BATHY = ',gdepw_0(ji,jj,mbathy(ji,jj)+1),sum(e3t_0(ji,jj,1:mbathy(ji,jj)))
329	STOP
330	! N_out = N_out + 1
331	! h_out(N_out) = - h_diff
332	ENDIF
333	DO jn=n1,n2-1
334	CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),ptab_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
335	ENDDO
336	ENDIF
337	ENDDO
338	ENDDO
339
340	! WARNING :
341	! ptab replaced by ptab_child in the following
342	! k1 k2 replaced by 1 jpk
343	! VERTICAL REFINEMENT END
344
345	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
346	! Add asselin part
347	DO jn = n1,n2-1
348	DO jk=1,jpk
349	DO jj=j1,j2
350	DO ji=i1,i2
351	IF( ptab_child(ji,jj,jk,jn) .NE. 0. ) THEN
352	tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &
353	& + atfp * ( ptab_child(ji,jj,jk,jn) &
354	& - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
355	ENDIF
356	ENDDO
357	ENDDO
358	ENDDO
359	ENDDO
360	ENDIF
361	DO jn = n1,n2-1
362	DO jk=1,jpk
363	DO jj=j1,j2
364	DO ji=i1,i2
365	IF( ptab_child(ji,jj,jk,jn) .NE. 0. ) THEN
366	tsn(ji,jj,jk,jn) = ptab_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
367	END IF
368	END DO
369	END DO
370	END DO
371	END DO
372	ENDIF
373	!
374	END SUBROUTINE updateTS
375
376	SUBROUTINE updateu( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
377	!!---------------------------------------------
378	!! * ROUTINE updateu *
379	!!---------------------------------------------
380	!!
381	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2, n1, n2
382	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: ptab
383	LOGICAL , INTENT(in ) :: before
384	!
385	INTEGER :: ji, jj, jk
386	REAL(wp) :: zrhoy
387	! VERTICAL REFINEMENT BEGIN
388	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
389	REAL(wp) :: h_in(k1:k2)
390	REAL(wp) :: h_out(1:jpk)
391	INTEGER :: N_in, N_out
392	REAL(wp) :: h_diff
393	REAL(wp) :: tabin(k1:k2)
394	! VERTICAL REFINEMENT END
395	!!---------------------------------------------
396	!
397	IF( before ) THEN
398	zrhoy = Agrif_Rhoy()
399	DO jk=k1,k2
400	DO jj=j1,j2
401	DO ji=i1,i2
402	ptab(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)
403	ptab(ji,jj,jk,2) = umask(ji,jj,jk) * zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk)
404	END DO
405	END DO
406	END DO
407	ELSE
408	! VERTICAL REFINEMENT BEGIN
409	ptab_child(:,:,:) = 0.
410
411	DO jj=j1,j2
412	DO ji=i1,i2
413	N_in = 0
414	DO jk=k1,k2 !k2=jpk of child grid
415	IF (ptab(ji,jj,jk,2) == 0) EXIT
416	N_in = N_in + 1
417	tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
418	h_in(N_in) = ptab(ji,jj,jk,2)/e2u(ji,jj)
419	ENDDO
420	N_out = 0
421	DO jk=1,jpk
422	IF (umask(ji,jj,jk) == 0) EXIT
423	N_out = N_out + 1
424	h_out(N_out) = e3u_n(ji,jj,jk)
425	ENDDO
426	IF (N_in * N_out > 0) THEN
427	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
428	! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
429	if (h_diff < 0.) then
430	print *,'CHECK YOUR BATHY ...'
431	stop
432	! else ! Extends with 0
433	! N_in = N_in + 1
434	! tabin(N_in) = 0.
435	! h_in(N_in) = h_diff
436	endif
437	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
438	ENDIF
439	ENDDO
440	ENDDO
441
442	! WARNING :
443	! ptab replaced by ptab_child in the following
444	! k1 k2 replaced by 1 jpk
445	! remove division by fs e3u (already done)
446	! VERTICAL REFINEMENT END
447
448	DO jk=1,jpk
449	DO jj=j1,j2
450	DO ji=i1,i2
451	!Following line now replaced by division higher up I think
452	! ptab_child(ji,jj,jk) = ptab_child(ji,jj,jk) / e2u(ji,jj)
453	!
454	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
455	ub(ji,jj,jk) = ub(ji,jj,jk) &
456	& + atfp * ( ptab_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
457	ENDIF
458	!
459	un(ji,jj,jk) = ptab_child(ji,jj,jk) * umask(ji,jj,jk)
460	END DO
461	END DO
462	END DO
463	ENDIF
464	!
465	END SUBROUTINE updateu
466
467	SUBROUTINE updatev( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
468	!!---------------------------------------------
469	!! * ROUTINE updatev *
470	!!---------------------------------------------
471	!!
472	INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2
473	INTEGER :: ji,jj,jk
474	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: ptab
475	LOGICAL :: before
476	!!
477	REAL(wp) :: zrhox
478	! VERTICAL REFINEMENT BEGIN
479	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
480	REAL(wp) :: h_in(k1:k2)
481	REAL(wp) :: h_out(1:jpk)
482	INTEGER :: N_in, N_out
483	REAL(wp) :: h_diff
484	REAL(wp) :: tabin(k1:k2)
485	! VERTICAL REFINEMENT END
486	!!---------------------------------------------
487	!
488	IF (before) THEN
489	zrhox = Agrif_Rhox()
490	DO jk=k1,k2
491	DO jj=j1,j2
492	DO ji=i1,i2
493	ptab(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)
494	ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk)
495	END DO
496	END DO
497	END DO
498	ELSE
499	! VERTICAL REFINEMENT BEGIN
500	ptab_child(:,:,:) = 0.
501
502	DO jj=j1,j2
503	DO ji=i1,i2
504	N_in = 0
505	DO jk=k1,k2
506	IF (ptab(ji,jj,jk,2) == 0) EXIT
507	N_in = N_in + 1
508	tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
509	h_in(N_in) = ptab(ji,jj,jk,2)/e1v(ji,jj)
510	ENDDO
511	N_out = 0
512	DO jk=1,jpk
513	IF (vmask(ji,jj,jk) == 0) EXIT
514	N_out = N_out + 1
515	h_out(N_out) = e3v_n(ji,jj,jk)
516	ENDDO
517	IF (N_in * N_out > 0) THEN
518	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
519	! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
520	if (h_diff < 0.) then
521	print *,'CHECK YOUR BATHY ...'
522	stop
523	! else ! Extends with 0
524	! N_in = N_in + 1
525	! tabin(N_in) = 0.
526	! h_in(N_in) = h_diff
527	endif
528	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
529	ENDIF
530	ENDDO
531	ENDDO
532
533	! WARNING :
534	! ptab replaced by ptab_child in the following
535	! k1 k2 replaced by 1 jpk
536	! remove division by fs e3v (already done)
537	! VERTICAL REFINEMENT END
538
539	DO jk=k1,jpk
540	DO jj=j1,j2
541	DO ji=i1,i2
542	!Following line now replaced by division higher up I think
543	! ptab_child(ji,jj,jk) = ptab_child(ji,jj,jk) / e1v(ji,jj)
544	!
545	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
546	vb(ji,jj,jk) = vb(ji,jj,jk) &
547	& + atfp * ( ptab_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
548	ENDIF
549	!
550	vn(ji,jj,jk) = ptab_child(ji,jj,jk) * vmask(ji,jj,jk)
551	END DO
552	END DO
553	END DO
554	ENDIF
555	END SUBROUTINE updatev
556
557
558	SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
559	!!---------------------------------------------
560	!! * ROUTINE updateu2d *
561	!!---------------------------------------------
562	INTEGER, INTENT(in) :: i1, i2, j1, j2
563	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
564	LOGICAL, INTENT(in) :: before
565	!!
566	INTEGER :: ji, jj, jk
567	REAL(wp) :: zrhoy
568	REAL(wp) :: zcorr
569	!!---------------------------------------------
570	!
571	IF (before) THEN
572	zrhoy = Agrif_Rhoy()
573	DO jj=j1,j2
574	DO ji=i1,i2
575	tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
576	END DO
577	END DO
578	ELSE
579	DO jj=j1,j2
580	DO ji=i1,i2
581	tabres(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * r1_e2u(ji,jj)
582	!
583	! Update "now" 3d velocities:
584	spgu(ji,jj) = 0._wp
585	DO jk=1,jpkm1
586	spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
587	END DO
588	spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj)
589	!
590	zcorr = tabres(ji,jj) - spgu(ji,jj)
591	DO jk=1,jpkm1
592	un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)
593	END DO
594	!
595	! Update barotropic velocities:
596	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
597	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
598	zcorr = tabres(ji,jj) - un_b(ji,jj)
599	ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
600	END IF
601	ENDIF
602	un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1)
603	!
604	! Correct "before" velocities to hold correct bt component:
605	spgu(ji,jj) = 0.e0
606	DO jk=1,jpkm1
607	spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
608	END DO
609	spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj)
610	!
611	zcorr = ub_b(ji,jj) - spgu(ji,jj)
612	DO jk=1,jpkm1
613	ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)
614	END DO
615	!
616	END DO
617	END DO
618	ENDIF
619	!
620	END SUBROUTINE updateu2d
621
622
623	SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
624	!!---------------------------------------------
625	!! * ROUTINE updatev2d *
626	!!---------------------------------------------
627	INTEGER, INTENT(in) :: i1, i2, j1, j2
628	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
629	LOGICAL, INTENT(in) :: before
630	!!
631	INTEGER :: ji, jj, jk
632	REAL(wp) :: zrhox
633	REAL(wp) :: zcorr
634	!!---------------------------------------------
635	!
636	IF (before) THEN
637	zrhox = Agrif_Rhox()
638	DO jj=j1,j2
639	DO ji=i1,i2
640	tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj)
641	END DO
642	END DO
643	ELSE
644	DO jj=j1,j2
645	DO ji=i1,i2
646	tabres(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * r1_e1v(ji,jj)
647	!
648	! Update "now" 3d velocities:
649	spgv(ji,jj) = 0.e0
650	DO jk=1,jpkm1
651	spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
652	END DO
653	spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj)
654	!
655	zcorr = tabres(ji,jj) - spgv(ji,jj)
656	DO jk=1,jpkm1
657	vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)
658	END DO
659	!
660	! Update barotropic velocities:
661	IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
662	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
663	zcorr = tabres(ji,jj) - vn_b(ji,jj)
664	vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
665	END IF
666	ENDIF
667	vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1)
668	!
669	! Correct "before" velocities to hold correct bt component:
670	spgv(ji,jj) = 0.e0
671	DO jk=1,jpkm1
672	spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
673	END DO
674	spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj)
675	!
676	zcorr = vb_b(ji,jj) - spgv(ji,jj)
677	DO jk=1,jpkm1
678	vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)
679	END DO
680	!
681	END DO
682	END DO
683	ENDIF
684	!
685	END SUBROUTINE updatev2d
686
687
688	SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
689	!!---------------------------------------------
690	!! * ROUTINE updateSSH *
691	!!---------------------------------------------
692	INTEGER, INTENT(in) :: i1, i2, j1, j2
693	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
694	LOGICAL, INTENT(in) :: before
695	!!
696	INTEGER :: ji, jj
697	!!---------------------------------------------
698	!
699	IF (before) THEN
700	DO jj=j1,j2
701	DO ji=i1,i2
702	tabres(ji,jj) = sshn(ji,jj)
703	END DO
704	END DO
705	ELSE
706	IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN
707	IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
708	DO jj=j1,j2
709	DO ji=i1,i2
710	sshb(ji,jj) = sshb(ji,jj) &
711	& + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
712	END DO
713	END DO
714	ENDIF
715	ENDIF
716	!
717	DO jj=j1,j2
718	DO ji=i1,i2
719	sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
720	END DO
721	END DO
722	ENDIF
723	!
724	END SUBROUTINE updateSSH
725
726
727	SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
728	!!---------------------------------------------
729	!! * ROUTINE updateub2b *
730	!!---------------------------------------------
731	INTEGER, INTENT(in) :: i1, i2, j1, j2
732	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
733	LOGICAL, INTENT(in) :: before
734	!!
735	INTEGER :: ji, jj
736	REAL(wp) :: zrhoy
737	!!---------------------------------------------
738	!
739	IF (before) THEN
740	zrhoy = Agrif_Rhoy()
741	DO jj=j1,j2
742	DO ji=i1,i2
743	tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
744	END DO
745	END DO
746	tabres = zrhoy * tabres
747	ELSE
748	DO jj=j1,j2
749	DO ji=i1,i2
750	ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj)
751	END DO
752	END DO
753	ENDIF
754	!
755	END SUBROUTINE updateub2b
756
757
758	SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
759	!!---------------------------------------------
760	!! * ROUTINE updatevb2b *
761	!!---------------------------------------------
762	INTEGER, INTENT(in) :: i1, i2, j1, j2
763	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
764	LOGICAL, INTENT(in) :: before
765	!!
766	INTEGER :: ji, jj
767	REAL(wp) :: zrhox
768	!!---------------------------------------------
769	!
770	IF (before) THEN
771	zrhox = Agrif_Rhox()
772	DO jj=j1,j2
773	DO ji=i1,i2
774	tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)
775	END DO
776	END DO
777	tabres = zrhox * tabres
778	ELSE
779	DO jj=j1,j2
780	DO ji=i1,i2
781	vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj)
782	END DO
783	END DO
784	ENDIF
785	!
786	END SUBROUTINE updatevb2b
787
788
789	SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
790	! currently not used
791	!!---------------------------------------------
792	!! * ROUTINE updateT *
793	!!---------------------------------------------
794	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
795	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
796	LOGICAL, iNTENT(in) :: before
797	!
798	INTEGER :: ji,jj,jk
799	REAL(wp) :: ztemp
800	!!---------------------------------------------
801
802	IF (before) THEN
803	DO jk=k1,k2
804	DO jj=j1,j2
805	DO ji=i1,i2
806	tabres(ji,jj,jk,1) = e1t(ji,jj)e2t(ji,jj)tmask(ji,jj,jk)
807	tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
808	tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
809	END DO
810	END DO
811	END DO
812	tabres(:,:,:,1)=tabres(:,:,:,1)Agrif_Rhox()Agrif_Rhoy()
813	tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
814	tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
815	ELSE
816	DO jk=k1,k2
817	DO jj=j1,j2
818	DO ji=i1,i2
819	IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
820	print ,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)e2t(ji,jj)*tmask(ji,jj,jk)
821	print ,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)tmask(ji,jj,jk)
822	print ,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)tmask(ji,jj,jk)
823	ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
824	print *,'CORR = ',ztemp-1.
825	print ,'NEW VALS = ',tabres(ji,jj,jk,2)ztemp,tabres(ji,jj,jk,3)*ztemp, &
826	tabres(ji,jj,jk,2)ztemptabres(ji,jj,jk,3)*ztemp
827	e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
828	e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
829	END IF
830	END DO
831	END DO
832	END DO
833	ENDIF
834	!
835	END SUBROUTINE update_scales
836
837	# if defined key_zdftke
838
839	SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
840	!!---------------------------------------------
841	!! * ROUTINE updateen *
842	!!---------------------------------------------
843	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
844	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
845	LOGICAL, INTENT(in) :: before
846	!!---------------------------------------------
847	!
848	IF (before) THEN
849	ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
850	ELSE
851	en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
852	ENDIF
853	!
854	END SUBROUTINE updateEN
855
856
857	SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
858	!!---------------------------------------------
859	!! * ROUTINE updateavt *
860	!!---------------------------------------------
861	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
862	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
863	LOGICAL, INTENT(in) :: before
864	!!---------------------------------------------
865	!
866	IF (before) THEN
867	ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
868	ELSE
869	avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
870	ENDIF
871	!
872	END SUBROUTINE updateAVT
873
874
875	SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
876	!!---------------------------------------------
877	!! * ROUTINE updateavm *
878	!!---------------------------------------------
879	INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
880	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
881	LOGICAL, INTENT(in) :: before
882	!!---------------------------------------------
883	!
884	IF (before) THEN
885	ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
886	ELSE
887	avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
888	ENDIF
889	!
890	END SUBROUTINE updateAVM
891
892	# endif /* key_zdftke */
893
894	#else
895	CONTAINS
896	SUBROUTINE agrif_opa_update_empty
897	!!---------------------------------------------
898	!! * ROUTINE agrif_opa_update_empty *
899	!!---------------------------------------------
900	WRITE(,) 'agrif_opa_update : You should not have seen this print! error?'
901	END SUBROUTINE agrif_opa_update_empty
902	#endif
903	END MODULE agrif_opa_update
904

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

Download in other formats: