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/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 8822

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

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

Download in other formats: