Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sshwzv_tam.F90 @ 5226

Last change on this file since 5226 was 2578, checked in by rblod, 13 years ago
first import of NEMOTAM 3.2.2
File size: 41.6 KB

Line
1	MODULE sshwzv_tam
2	#if defined key_tam
3	!!==============================================================================
4	!! * MODULE sshwzv *
5	!! Ocean dynamics : sea surface height and vertical velocity
6	!!==============================================================================
7	!! History of the direct module:
8	!! 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code
9	!! History of the TAM module:
10	!! 3.2 ! 2010-04 (F. Vigilant) Original code
11	!!----------------------------------------------------------------------
12
13	!!----------------------------------------------------------------------
14	!! ssh_wzv : after ssh & now vertical velocity
15	!! ssh_nxt : filter ans swap the ssh arrays
16	!!----------------------------------------------------------------------
17	!! * Modules used
18	USE par_kind , ONLY: & ! Precision variables
19	& wp
20	USE par_oce , ONLY: & ! Ocean space and time domain variables
21	& jpi, &
22	& jpj, &
23	& jpk, &
24	& jpim1, &
25	& jpjm1, &
26	& jpkm1, &
27	& jpiglo
28	USE in_out_manager, ONLY: & ! I/O manager
29	& lwp, &
30	& numout, &
31	& nit000, &
32	& nitend, &
33	& ln_ctl
34	USE dom_oce , ONLY: & ! Ocean space and time domain
35	& n_cla, &
36	& e2u, &
37	& e1v, &
38	& e1t, &
39	& e1f, &
40	& e2t, &
41	& e2f, &
42	# if defined key_vvl
43	& e3t_1, &
44	# else
45	# if defined key_zco
46	& e3t_0, &
47	# else
48	& e3t, &
49	# endif
50	# endif
51	# if defined key_zco
52	! & e3u_0, & scale factor is identical to e3t_0
53	! & e3v_0, &
54	# else
55	& e3u, &
56	& e3v, &
57	& e3f, &
58	# endif
59	& lk_vvl, &
60	& rdt, &
61	& neuler, &
62	& atfp, &
63	& tmask, &
64	& mig, &
65	& mjg, &
66	& nldi, &
67	& nldj, &
68	& nlei, &
69	& nlej
70	USE prtctl , ONLY: & ! Print control
71	& prt_ctl
72	USE phycst , ONLY: & ! Physical constants
73	& rau0
74	# if defined key_obc
75	USE obc_par , ONLY: & ! Open boundary conditions
76	& lp_obc_east, &
77	& lp_obc_west, &
78	& lp_obc_north, &
79	& lp_obc_south
80	USE obc_oce , ONLY: & ! Open boundary conditions
81	& nie0p1, &
82	& nie1p1, &
83	& njn0p1, &
84	& njn1p1, &
85	& nje0, &
86	& nje1, &
87	& niw0, &
88	& niw1, &
89	& njw0, &
90	& njw1, &
91	& nin0, &
92	& nin1, &
93	& nis0, &
94	& nis1, &
95	& njs0, &
96	& njs1
97	# endif
98	USE lbclnk , ONLY: & ! Lateral boundary conditions
99	& lbc_lnk
100
101	USE lbclnk_tam , ONLY: & ! TAM lateral boundary conditions
102	& lbc_lnk_adj
103	USE divcur_tam , ONLY: & ! TAM horizontal divergence and relative
104	& div_cur_tan, & ! vorticity
105	& div_cur_adj ! vorticity
106	USE cla_div_tam , ONLY: &
107	& div_cla_tan, & !
108	& div_cla_adj !
109	USE oce_tam , ONLY: & ! TAM ocean dynamics and tracers variables
110	& un_tl, &
111	& un_ad, &
112	& vn_tl, &
113	& vn_ad, &
114	& wn_tl, &
115	& wn_ad, &
116	& wn_tl, &
117	& wn_ad, &
118	& hdivn_tl, &
119	& hdivn_ad, &
120	& hdivb_tl, &
121	& hdivb_ad, &
122	& rotn_tl, &
123	& rotn_ad, &
124	& rotb_tl, &
125	& rotb_ad, &
126	& sshb_tl, &
127	& sshb_ad, &
128	& sshn_tl, &
129	& sshn_ad, &
130	& ssha_tl, &
131	& ssha_ad
132	USE sbc_oce_tam , ONLY: & ! surface variables
133	& emp_tl, &
134	& emp_ad
135	USE gridrandom , ONLY: & ! Random Gaussian noise on grids
136	& grid_random
137	USE dotprodfld, ONLY: & ! Computes dot product for 3D and 2D fields
138	& dot_product
139	USE paresp , ONLY: & ! Normalized energy weights
140	& wesp_ssh
141	USE tstool_tam , ONLY: &
142	& prntst_adj, &
143	& stdssh, &
144	& stdu, &
145	& stdv
146
147	IMPLICIT NONE
148	PRIVATE
149
150	PUBLIC ssh_wzv_tan ! called by step.F90
151	PUBLIC ssh_nxt_tan ! called by step.F90
152	PUBLIC ssh_wzv_adj ! called by step.F90
153	PUBLIC ssh_nxt_adj ! called by step.F90
154	PUBLIC ssh_wzv_adj_tst ! called by tamtst.F90
155	PUBLIC ssh_nxt_adj_tst ! called by tamtst.F90
156
157	!! * Substitutions
158	# include "domzgr_substitute.h90"
159	# include "vectopt_loop_substitute.h90"
160
161	CONTAINS
162
163	SUBROUTINE ssh_wzv_tan( kt , kdum )
164	!!----------------------------------------------------------------------
165	!! * ROUTINE ssh_wzv_tan *
166	!!
167	!! ** Purpose of direct routine :
168	!! compute the after ssh (ssha), the now vertical velocity
169	!! and update the now vertical coordinate (lk_vvl=T).
170	!!
171	!! ** Method : -
172	!!
173	!! - Using the incompressibility hypothesis, the vertical
174	!! velocity is computed by integrating the horizontal divergence
175	!! from the bottom to the surface minus the scale factor evolution.
176	!! The boundary conditions are w=0 at the bottom (no flux) and.
177	!!
178	!! ** action : ssha : after sea surface height
179	!! wn : now vertical velocity
180	!! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height
181	!! at u-, v-, f-point s
182	!! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
183	!!----------------------------------------------------------------------
184	!!
185	INTEGER, INTENT(in) :: kt ! time step
186	!!
187	INTEGER :: jk ! dummy loop indices
188	REAL(wp) :: z2dt, zraur ! temporary scalars
189	REAL(wp), DIMENSION(jpi,jpj) :: zhdivtl ! 2D workspace
190	INTEGER, OPTIONAL :: kdum ! dummy argument to compute only vertical velocity
191	!!----------------------------------------------------------------------
192
193	IF( kt == nit000 ) THEN
194	IF(lwp) WRITE(numout,*)
195	IF(lwp) WRITE(numout,*) 'ssh_wzv_tan : after sea surface height and now vertical velocity '
196	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
197	!
198	wn_tl(:,:,jpk) = 0.0_wp ! bottom boundary condition: w=0 (set once for all)
199	!
200	IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)
201	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
202	CALL abort
203	ENDIF
204	!
205	ENDIF
206	! !------------------------------!
207	IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)
208	!------------------------------!
209	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
210	CALL abort
211	!
212	ENDIF
213
214	CALL div_cur_tan( kt ) ! Horizontal divergence & Relative vorticity
215
216	IF ( .NOT. PRESENT(kdum) ) THEN ! jump
217	IF( n_cla == 1 ) CALL div_cla_tan( kt ) ! Cross Land Advection (Update Hor. divergence)
218	ENDIF
219
220	! set time step size (Euler/Leapfrog)
221	z2dt = 2. * rdt
222	IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt
223
224	zraur = 1. / rau0
225
226	IF ( .NOT. PRESENT(kdum) ) THEN ! jump ssh computing
227	! !------------------------------!
228	! ! After Sea Surface Height !
229	! !------------------------------!
230	zhdivtl(:,:) = 0.0_wp
231	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
232	zhdivtl(:,:) = zhdivtl(:,:) + fse3t(:,:,jk) * hdivn_tl(:,:,jk)
233	END DO
234
235	! ! Sea surface elevation time stepping
236	ssha_tl(:,:) = ( sshb_tl(:,:) - z2dt * ( zraur * emp_tl(:,:) + zhdivtl(:,:) ) ) * tmask(:,:,1)
237
238	#if defined key_agrif
239	! CALL agrif_ssh_tan(kt)
240	IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
241	CALL abort
242	#endif
243	#if defined key_obc
244	IF ( Agrif_Root() ) THEN
245	ssha_tl(:,:) = ssha_tl(:,:) * obctmsk(:,:)
246	CALL lbc_lnk( ssha_tl, 'T', 1.0_wp ) ! absolutly compulsory !! (jmm)
247	ENDIF
248	#endif
249
250	! ! Sea Surface Height at u-,v- and f-points (vvl case only)
251	IF( lk_vvl ) THEN ! (required only in key_vvl case)
252	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
253	CALL abort
254	ENDIF
255
256	ENDIF
257	! !------------------------------!
258	! ! Now Vertical Velocity !
259	! !------------------------------!
260	! ! integrate from the bottom the hor. divergence
261	DO jk = jpkm1, 1, -1
262	wn_tl(:,:,jk) = wn_tl(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn_tl(:,:,jk)
263	END DO
264	!
265	END SUBROUTINE ssh_wzv_tan
266
267	SUBROUTINE ssh_nxt_tan( kt )
268	!!----------------------------------------------------------------------
269	!! * ROUTINE ssh_nxt_tan *
270	!!
271	!! ** Purpose of the direct :
272	!! achieve the sea surface height time stepping by
273	!! applying Asselin time filter and swapping the arrays
274	!! ssha already computed in ssh_wzv
275	!!
276	!! ** Method : - apply Asselin time fiter to now ssh and swap :
277	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
278	!! sshn = ssha
279	!!
280	!! ** action : - sshb, sshn : before & now sea surface height
281	!! ready for the next time step
282	!!----------------------------------------------------------------------
283	INTEGER, INTENT( in ) :: kt ! ocean time-step index
284	!!
285	INTEGER :: ji, jj ! dummy loop indices
286	!!----------------------------------------------------------------------
287
288	IF( kt == nit000 ) THEN
289	IF(lwp) WRITE(numout,*)
290	IF(lwp) WRITE(numout,*) 'ssh_nxt_tan : next sea surface height (Asselin time filter + swap)'
291	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
292	ENDIF
293
294	! Time filter and swap of the ssh
295	! -------------------------------
296
297	IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points
298	! ! ---------------------- !
299	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
300	CALL abort
301	!
302	ELSE ! fixed levels : ssh at t-point only
303	! ! ------------ !
304	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter
305	sshn_tl(:,:) = ssha_tl(:,:) ! now <-- after (before already = now)
306	!
307	ELSE ! Leap-Frog time-stepping: Asselin filter + swap
308	DO jj = 1, jpj
309	DO ji = 1, jpi ! before <-- now filtered
310	sshb_tl(ji,jj) = sshn_tl(ji,jj) + atfp * ( sshb_tl(ji,jj) - 2 * sshn_tl(ji,jj) + ssha_tl(ji,jj) )
311	sshn_tl(ji,jj) = ssha_tl(ji,jj) ! now <-- after
312	END DO
313	END DO
314	ENDIF
315	ENDIF
316	!
317	#if defined key_agrif
318	! Update velocity at AGRIF zoom boundaries
319	!IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn_tan( kt )
320	IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
321	CALL abort
322	#endif
323	!
324	END SUBROUTINE ssh_nxt_tan
325
326	SUBROUTINE ssh_wzv_adj( kt , kdum )
327	!!----------------------------------------------------------------------
328	!! * ROUTINE ssh_wzv_adj *
329	!!
330	!! ** Purpose of direct routine :
331	!! compute the after ssh (ssha), the now vertical velocity
332	!! and update the now vertical coordinate (lk_vvl=T).
333	!!
334	!! ** Method : -
335	!!
336	!! - Using the incompressibility hypothesis, the vertical
337	!! velocity is computed by integrating the horizontal divergence
338	!! from the bottom to the surface minus the scale factor evolution.
339	!! The boundary conditions are w=0 at the bottom (no flux) and.
340	!!
341	!! ** action : ssha : after sea surface height
342	!! wn : now vertical velocity
343	!! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height
344	!! at u-, v-, f-point s
345	!! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points
346	!!----------------------------------------------------------------------
347	!!
348	INTEGER, INTENT(in) :: kt ! time step
349	!!
350	INTEGER :: jk ! dummy loop indices
351	REAL(wp) :: z2dt, zraur ! temporary scalars
352	REAL(wp), DIMENSION(jpi,jpj) :: zhdivad ! 2D workspace
353	INTEGER, OPTIONAL :: kdum ! dummy argument to compute only vertical velocity
354	!!----------------------------------------------------------------------
355
356	! adjoint variable initialization
357	zhdivad = 0.0_wp
358
359	IF( kt == nitend ) THEN
360	IF(lwp) WRITE(numout,*)
361	IF(lwp) WRITE(numout,*) 'ssh_wzv_adj : after sea surface height and now vertical velocity '
362	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
363	ENDIF
364
365	! !------------------------------!
366	! ! Now Vertical Velocity !
367	! !------------------------------!
368	! ! integrate from the bottom the hor. divergence
369	DO jk = 1, jpkm1
370	hdivn_ad(:,:,jk ) = hdivn_ad(:,:,jk ) - fse3t_n(:,:,jk) * wn_ad(:,:,jk)
371	wn_ad( :,:,jk+1) = wn_ad(:,:,jk+1) + wn_ad(:,:,jk)
372	wn_ad( :,:,jk ) = 0.0_wp
373	END DO
374	!
375	! set time step size (Euler/Leapfrog)
376	z2dt = 2. * rdt
377	IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt
378
379	zraur = 1. / rau0
380
381	IF ( .NOT. PRESENT(kdum) ) THEN ! jump ssh computing
382	! !------------------------------!
383	! ! After Sea Surface Height !
384	! !------------------------------!
385	! ! Sea Surface Height at u-,v- and f-points (vvl case only)
386	IF( lk_vvl ) THEN ! (required only in key_vvl case)
387	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
388	CALL abort
389	ENDIF
390
391	#if defined key_obc
392	IF ( Agrif_Root() ) THEN
393	CALL lbc_lnk_ad( ssha_ad, 'T', 1.0_wp ) ! absolutly compulsory !! (jmm)
394	ssha_ad(:,:) = ssha_ad(:,:) * obctmsk(:,:)
395	ENDIF
396	#endif
397	#if defined key_agrif
398	! CALL agrif_ssh_adj(kt)
399	IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
400	CALL abort
401	#endif ! ! Sea surface elevation time stepping
402	sshb_ad(:,:) = sshb_ad(:,:) + ssha_ad(:,:)* tmask(:,:,1)
403	emp_ad( :,:) = emp_ad(:,:) - z2dt * zraur * ssha_ad(:,:) * tmask(:,:,1)
404	zhdivad(:,:) = zhdivad(:,:) - z2dt * tmask(:,:,1) * ssha_ad(:,:)
405	ssha_ad(:,:) = 0.0_wp
406
407	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
408	hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) + fse3t(:,:,jk) * zhdivad(:,:)
409	END DO
410
411	ENDIF
412
413	IF ( .NOT. PRESENT(kdum) ) THEN ! jump
414	IF( n_cla == 1 ) CALL div_cla_adj( kt ) ! Cross Land Advection (Update Hor. divergence)
415	ENDIF
416
417	CALL div_cur_adj( kt ) ! Horizontal divergence & Relative vorticity
418	!
419	! !------------------------------!
420	IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)
421	!------------------------------!
422	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
423	CALL abort
424	!
425	ENDIF
426	!
427	IF( kt == nit000 ) THEN
428	!
429	IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)
430	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
431	CALL abort
432	ENDIF
433	!
434	wn_ad(:,:,jpk) = 0.0_wp ! bottom boundary condition: w=0 (set once for all)
435	!
436	ENDIF
437	!
438	END SUBROUTINE ssh_wzv_adj
439
440	SUBROUTINE ssh_nxt_adj( kt )
441	!!----------------------------------------------------------------------
442	!! * ROUTINE ssh_nxt_adj *
443	!!
444	!! ** Purpose of the direct :
445	!! achieve the sea surface height time stepping by
446	!! applying Asselin time filter and swapping the arrays
447	!! ssha already computed in ssh_wzv
448	!!
449	!! ** Method : - apply Asselin time fiter to now ssh and swap :
450	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
451	!! sshn = ssha
452	!!
453	!! ** action : - sshb, sshn : before & now sea surface height
454	!! ready for the next time step
455	!!----------------------------------------------------------------------
456	INTEGER, INTENT( in ) :: kt ! ocean time-step index
457	!!
458	INTEGER :: ji, jj ! dummy loop indices
459	!!----------------------------------------------------------------------
460
461	IF( kt == nitend ) THEN
462	IF(lwp) WRITE(numout,*)
463	IF(lwp) WRITE(numout,*) 'ssh_nxt_adj : next sea surface height (Asselin time filter + swap)'
464	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
465	ENDIF
466
467	! Time filter and swap of the ssh
468	! -------------------------------
469	#if defined key_agrif
470	! Update velocity at AGRIF zoom boundaries
471	!IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn_adj( kt )
472	IF (lwp) WRITE(numout,*) 'key_agrif not available yet'
473	CALL abort
474	#endif
475	IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points
476	! ! ---------------------- !
477	IF (lwp) WRITE(numout,*) 'lk_vvl not available yet'
478	CALL abort
479	!
480	ELSE ! fixed levels : ssh at t-point only
481	! ! ------------ !
482	IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter
483	ssha_ad(:,:) = ssha_ad(:,:) + sshn_ad(:,:)
484	sshn_ad(:,:) = 0.0_wp
485	!
486	ELSE ! Leap-Frog time-stepping: Asselin filter + swap
487	DO jj = 1, jpj
488	DO ji = 1, jpi ! before <-- now filtered
489	ssha_ad(ji,jj) = ssha_ad(ji,jj) + sshn_ad(ji,jj)
490	sshn_ad(ji,jj) = (1.0_wp - 2.0 * atfp) * sshb_ad(ji,jj)
491	ssha_ad(ji,jj) = ssha_ad(ji,jj) + atfp * sshb_ad(ji,jj)
492	sshb_ad(ji,jj) = atfp * sshb_ad(ji,jj)
493	END DO
494	END DO
495	ENDIF
496	ENDIF
497
498	!
499	END SUBROUTINE ssh_nxt_adj
500
501	SUBROUTINE ssh_wzv_adj_tst( kumadt )
502	!!-----------------------------------------------------------------------
503	!!
504	!! * ROUTINE ssh_wzv_adj_tst : TEST OF wzv_adj *
505	!!
506	!! ** Purpose : Test the adjoint routine.
507	!!
508	!! ** Method : Verify the scalar product
509	!!
510	!! ( L dx )^T W dy = dx^T L^T W dy
511	!!
512	!! where L = tangent routine
513	!! L^T = adjoint routine
514	!! W = diagonal matrix of scale factors
515	!! dx = input perturbation (random field)
516	!! dy = L dx
517	!!
518	!! ** Action : Separate tests are applied for the following dx and dy:
519	!!
520	!! dx = ( un_tl, vn_tl, hdivn_tl, rotn_tl, emp_tl, sshb_tl ) and
521	!! dy = ( hdivn_tl, hdivb_tl, rotn_tl, rotb_tl, wn_tl, ssha_tl )
522	!!
523	!! History :
524	!! ! 2010-04 (F. Vigilant)
525	!!-----------------------------------------------------------------------
526
527	!! * Modules used
528	!! * Arguments
529	INTEGER, INTENT(IN) :: &
530	& kumadt ! Output unit
531
532	!! * Local declarations
533	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
534	& zun_tlin, & ! Tangent input: now u-velocity
535	& zvn_tlin, & ! Tangent input: now v-velocity
536	& zhdivn_tlin, & ! Tangent input: now horizontal divergence
537	& zrotn_tlin, & ! Tangent input: now horizontal divergence
538	& zhdivn_tlout, & ! Tangent output: now horizontal divergence
539	& zrotn_tlout, & ! Tangent output: now horizontal divergence
540	& zrotb_tlout, & ! Tangent output: now horizontal divergence
541	& zhdivb_tlout, & ! Tangent output: now horizontal divergence
542	& zwn_tlout, & ! Tangent output: now w-velocity
543	& zwn_adin, & ! Adjoint input: now w-velocity
544	& zhdivn_adout, & ! Adjoint output: now horizontal divergence
545	& zrotn_adin, & ! Adjoint input: now horizontal divergence
546	& zrotn_adout, & ! Adjoint output: now horizontal divergence
547	& zrotb_adin, & ! Adjoint input: now horizontal divergence
548	& zhdivn_adin, & ! Adjoint input: now horizontal divergence
549	& zhdivb_adin, & ! Adjoint output: now horizontal divergence
550	& zun_adout, & ! Adjoint output: now horizontal divergence
551	& zvn_adout, & ! Adjoint output: now horizontal divergence
552	& znu, & ! 3D random field for u
553	& znv ! 3D random field for v
554	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
555	& zsshb_tlin, & ! Tangent input: before SSH
556	& zssha_tlout, & ! Tangent input: before SSH
557	& zsshb_adout, & ! Adjoint output: before SSH
558	& zssha_adin, & ! Adjoint output: before SSH
559	& zemp_tlin, & ! Tangent input: EmP
560	& zemp_adout, & ! Adjoint output: EmP
561	& znssh, & ! 2D random field for SSH
562	& znemp ! 2D random field for EmP
563
564	INTEGER :: &
565	& ji, & ! dummy loop indices
566	& jj, &
567	& jk
568	INTEGER, DIMENSION(jpi,jpj) :: &
569	& iseed_2d ! 2D seed for the random number generator
570	REAL(KIND=wp) :: &
571	! random field standard deviation for:
572	& zstdssh, & ! SSH
573	& zstdemp, & ! EMP
574	& zsp1, & ! scalar product involving the tangent routine
575	& zsp2, & ! scalar product involving the adjoint routine
576	& zsp2_1, & ! scalar product components
577	& zsp2_2, &
578	& zsp2_3, &
579	& zsp2_4, &
580	& zsp2_5, &
581	& zsp2_6, &
582	& z2dt, & ! temporary scalars
583	& zraur
584	CHARACTER (LEN=14) :: &
585	& cl_name
586
587	! Allocate memory
588
589	ALLOCATE( &
590	& zhdivn_tlin(jpi,jpj,jpk), &
591	& zhdivb_tlout(jpi,jpj,jpk), &
592	& zhdivn_tlout(jpi,jpj,jpk), &
593	& zrotn_tlin(jpi,jpj,jpk), &
594	& zrotn_tlout(jpi,jpj,jpk), &
595	& zrotb_tlout(jpi,jpj,jpk), &
596	& zwn_tlout(jpi,jpj,jpk), &
597	& zwn_adin(jpi,jpj,jpk), &
598	& zhdivn_adout(jpi,jpj,jpk), &
599	& zhdivb_adin(jpi,jpj,jpk), &
600	& zrotn_adin(jpi,jpj,jpk), &
601	& zrotn_adout(jpi,jpj,jpk), &
602	& zrotb_adin(jpi,jpj,jpk), &
603	& zhdivn_adin(jpi,jpj,jpk), &
604	& zun_tlin(jpi,jpj,jpk), &
605	& zvn_tlin(jpi,jpj,jpk), &
606	& zun_adout(jpi,jpj,jpk), &
607	& zvn_adout(jpi,jpj,jpk), &
608	& znu(jpi,jpj,jpk), &
609	& znv(jpi,jpj,jpk) &
610	& )
611	ALLOCATE( &
612	& zsshb_tlin(jpi,jpj), &
613	& zsshb_adout(jpi,jpj), &
614	& zssha_tlout(jpi,jpj), &
615	& zssha_adin(jpi,jpj), &
616	& zemp_tlin(jpi,jpj), &
617	& zemp_adout(jpi,jpj), &
618	& znssh(jpi,jpj), &
619	& znemp(jpi,jpj) &
620	& )
621
622
623	! Initialize constants
624
625	z2dt = 2.0_wp * rdt ! time step: leap-frog
626	zraur = 1.0_wp / rau0 ! inverse density of pure water (m3/kg)
627
628	zhdivn_tlin(:,:,:) = 0.0_wp
629	zrotn_tlin(:,:,:) = 0.0_wp
630	zemp_tlin(:,:) = 0.0_wp
631	zsshb_tlin(:,:) = 0.0_wp
632	zun_tlin (:,:,:) = 0.0_wp
633	zvn_tlin (:,:,:) = 0.0_wp
634
635	zhdivn_tlout(:,:,:) = 0.0_wp
636	zhdivb_tlout(:,:,:) = 0.0_wp
637	zrotn_tlout(:,:,:) = 0.0_wp
638	zrotb_tlout(:,:,:) = 0.0_wp
639	zwn_tlout(:,:,:) = 0.0_wp
640	zssha_tlout(:,:) = 0.0_wp
641
642	zhdivn_adin(:,:,:) = 0.0_wp
643	zhdivb_adin(:,:,:) = 0.0_wp
644	zrotn_adin(:,:,:) = 0.0_wp
645	zrotb_adin(:,:,:) = 0.0_wp
646	zwn_adin(:,:,:) = 0.0_wp
647	zssha_adin(:,:) = 0.0_wp
648
649	zhdivn_adout(:,:,:) = 0.0_wp
650	zrotn_adout(:,:,:) = 0.0_wp
651	zemp_adout(:,:) = 0.0_wp
652	zsshb_adout(:,:) = 0.0_wp
653	zun_adout (:,:,:) = 0.0_wp
654	zvn_adout (:,:,:) = 0.0_wp
655
656	un_tl (:,:,:) = 0.0_wp
657	vn_tl (:,:,:) = 0.0_wp
658	hdivn_tl(:,:,:) = 0.0_wp
659	hdivb_tl(:,:,:) = 0.0_wp
660	rotn_tl (:,:,:) = 0.0_wp
661	rotb_tl (:,:,:) = 0.0_wp
662	wn_tl(:,:,:) = 0.0_wp
663	ssha_tl(:,:) = 0.0_wp
664	sshb_tl(:,:) = 0.0_wp
665	emp_tl(:,:) = 0.0_wp
666
667	un_ad (:,:,:) = 0.0_wp
668	vn_ad (:,:,:) = 0.0_wp
669	hdivn_ad(:,:,:) = 0.0_wp
670	hdivb_ad(:,:,:) = 0.0_wp
671	rotn_ad (:,:,:) = 0.0_wp
672	rotb_ad (:,:,:) = 0.0_wp
673	wn_ad(:,:,:) = 0.0_wp
674	sshb_ad(:,:) = 0.0_wp
675	ssha_ad(:,:) = 0.0_wp
676	emp_ad(:,:) = 0.0_wp
677
678	!=============================================================
679	! 1) dx = ( un_tl, vn_tl, emp_tl, sshb_tl ) and dy = ( wn_tl )
680	! - or -
681	! 2) dx = ( hdivn_tl ) and dy = ( wn_tl )
682	!=============================================================
683
684	!--------------------------------------------------------------------
685	! Initialize the tangent input with random noise: dx
686	!--------------------------------------------------------------------
687
688	DO jj = 1, jpj
689	DO ji = 1, jpi
690	iseed_2d(ji,jj) = - ( 596035 + &
691	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
692	END DO
693	END DO
694	CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu )
695
696	DO jj = 1, jpj
697	DO ji = 1, jpi
698	iseed_2d(ji,jj) = - ( 523432 + &
699	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
700	END DO
701	END DO
702	CALL grid_random( iseed_2d, znv, 'V', 0.0_wp, stdv )
703
704	DO jk = 1, jpk
705	DO jj = nldj, nlej
706	DO ji = nldi, nlei
707	zun_tlin(ji,jj,jk) = znu(ji,jj,jk)
708	zvn_tlin(ji,jj,jk) = znv(ji,jj,jk)
709	END DO
710	END DO
711	END DO
712
713	DO jj = 1, jpj
714	DO ji = 1, jpi
715	iseed_2d(ji,jj) = - ( 785483 + &
716	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
717	END DO
718	END DO
719	CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh )
720
721	DO jj = 1, jpj
722	DO ji = 1, jpi
723	iseed_2d(ji,jj) = - ( 358606 + &
724	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
725	END DO
726	END DO
727	CALL grid_random( iseed_2d, znemp, 'T', 0.0_wp, stdssh )
728
729	DO jj = nldj, nlej
730	DO ji = nldi, nlei
731	zsshb_tlin(ji,jj) = znssh(ji,jj)
732	zemp_tlin (ji,jj) = znemp(ji,jj) / ( z2dt * zraur )
733	END DO
734	END DO
735
736	un_tl(:,:,:) = zun_tlin(:,:,:)
737	vn_tl(:,:,:) = zvn_tlin(:,:,:)
738	CALL div_cur_tan( nit000 ) ! Generate noise hdiv/rot fields
739
740	DO jk = 1, jpk
741	DO jj = nldj, nlej
742	DO ji = nldi, nlei
743	zhdivn_tlin(ji,jj,jk) = 0.5_wp * hdivn_tl(ji,jj,jk)
744	zrotn_tlin (ji,jj,jk) = 0.5_wp * rotn_tl (ji,jj,jk)
745	END DO
746	END DO
747	END DO
748
749	! re-initialization to zero
750	un_tl (:,:,:) = 0.0_wp
751	vn_tl (:,:,:) = 0.0_wp
752	hdivb_tl(:,:,:) = 0.0_wp
753	hdivn_tl(:,:,:) = 0.0_wp
754	rotb_tl (:,:,:) = 0.0_wp
755	rotn_tl (:,:,:) = 0.0_wp
756
757	!--------------------------------------------------------------------
758	! Call the tangent routine: dy = L dx
759	!--------------------------------------------------------------------
760
761	hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:)
762	rotn_tl(:,:,:) = zrotn_tlin(:,:,:)
763	sshb_tl(:,:) = zsshb_tlin(:,:)
764	emp_tl (:,:) = zemp_tlin (:,:)
765	un_tl(:,:,:) = zun_tlin(:,:,:)
766	vn_tl(:,:,:) = zvn_tlin(:,:,:)
767
768	CALL ssh_wzv_tan( nit000+1 )
769
770	zwn_tlout(:,:,:) = wn_tl(:,:,:)
771	zssha_tlout(:,: ) = ssha_tl(:,:)
772	zhdivb_tlout(:,:,:) = hdivb_tl(:,:,:)
773	zhdivn_tlout(:,:,:) = hdivn_tl(:,:,:)
774	zrotb_tlout(:,:,:) = rotb_tl(:,:,:)
775	zrotn_tlout(:,:,:) = rotn_tl(:,:,:)
776	!--------------------------------------------------------------------
777	! Initialize the adjoint variables: dy^* = W dy
778	!--------------------------------------------------------------------
779
780	DO jk = 1, jpk
781	DO jj = nldj, nlej
782	DO ji = nldi, nlei
783	zwn_adin(ji,jj,jk) = zwn_tlout(ji,jj,jk) &
784	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
785	& * tmask(ji,jj,jk)
786	END DO
787	END DO
788	END DO
789	DO jj = nldj, nlej
790	DO ji = nldi, nlei
791	zssha_adin(ji,jj) = zssha_tlout(ji,jj) &
792	& * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
793	& * tmask(ji,jj,1)
794	END DO
795	END DO
796	DO jk = 1, jpk
797	DO jj = nldj, nlej
798	DO ji = nldi, nlei
799	zhdivb_adin(ji,jj,jk) = zhdivb_tlout(ji,jj,jk) &
800	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
801	& * tmask(ji,jj,jk)
802	zhdivn_adin(ji,jj,jk) = zhdivn_tlout(ji,jj,jk) &
803	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
804	& * tmask(ji,jj,jk)
805	END DO
806	END DO
807	END DO
808	DO jk = 1, jpk
809	DO jj = nldj, nlej
810	DO ji = nldi, nlei
811	zrotb_adin(ji,jj,jk) = zrotb_tlout(ji,jj,jk) &
812	& * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
813	zrotn_adin(ji,jj,jk) = zrotn_tlout(ji,jj,jk) &
814	& * e1f(ji,jj) * e2f(ji,jj) * fse3f(ji,jj,jk)
815	END DO
816	END DO
817	END DO
818
819	!--------------------------------------------------------------------
820	! Compute the scalar product: ( L dx )^T W dy
821	!--------------------------------------------------------------------
822
823
824	zsp1 = DOT_PRODUCT( zwn_tlout, zwn_adin ) + DOT_PRODUCT( zssha_tlout, zssha_adin ) &
825	& + DOT_PRODUCT( zhdivb_tlout, zhdivb_adin ) + DOT_PRODUCT( zhdivn_tlout, zhdivn_adin ) &
826	& + DOT_PRODUCT( zrotb_tlout, zrotb_adin ) + DOT_PRODUCT( zrotn_tlout, zrotn_adin )
827	!--------------------------------------------------------------------
828	! Call the adjoint routine: dx^* = L^T dy^*
829	!--------------------------------------------------------------------
830
831	wn_ad(:,:,:) = zwn_adin(:,:,:)
832	ssha_ad(:,:) = zssha_adin(:,:)
833	hdivb_ad(:,:,:) = zhdivb_adin(:,:,:)
834	hdivn_ad(:,:,:) = zhdivn_adin(:,:,:)
835	rotb_ad(:,:,:) = zrotb_adin(:,:,:)
836	rotn_ad(:,:,:) = zrotn_adin(:,:,:)
837
838	CALL ssh_wzv_adj( nit000+1 )
839
840	zrotn_adout(:,:,:) = rotn_ad(:,:,:)
841	zhdivn_adout(:,:,:) = hdivn_ad(:,:,:)
842	zsshb_adout(:,:) = sshb_ad(:,:)
843	zemp_adout (:,:) = emp_ad (:,:)
844	zun_adout(:,:,:) = un_ad(:,:,:)
845	zvn_adout(:,:,:) = vn_ad(:,:,:)
846
847	!--------------------------------------------------------------------
848	! Compute the scalar product: dx^T L^T W dy
849	!--------------------------------------------------------------------
850
851	zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout )
852	zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout )
853	zsp2_3 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout )
854	zsp2_4 = DOT_PRODUCT( zemp_tlin, zemp_adout )
855	zsp2_5 = DOT_PRODUCT( zsshb_tlin, zsshb_adout )
856	zsp2_6 = DOT_PRODUCT( zrotn_tlin, zrotn_adout )
857
858	zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 + zsp2_6
859
860	! Compare the scalar products
861	! 14 char:'12345678901234'
862	cl_name = 'sshwzv_adj '
863	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
864
865	END SUBROUTINE ssh_wzv_adj_tst
866
867	SUBROUTINE ssh_nxt_adj_tst( kumadt )
868	!!-----------------------------------------------------------------------
869	!!
870	!! * ROUTINE ssh_nxt_adj_tst : TEST OF nxt_adj *
871	!!
872	!! ** Purpose : Test the adjoint routine.
873	!!
874	!! ** Method : Verify the scalar product
875	!!
876	!! ( L dx )^T W dy = dx^T L^T W dy
877	!!
878	!! where L = tangent routine
879	!! L^T = adjoint routine
880	!! W = diagonal matrix of scale factors
881	!! dx = input perturbation (random field)
882	!! dy = L dx
883	!!
884	!! ** Action : Separate tests are applied for the following dx and dy:
885	!!
886	!! dx = ( sshb_tl, sshn_tl, ssha_tl ) and
887	!! dy = ( ssb_tl, sshn_tl )
888	!!
889	!! History :
890	!! ! 2010-05 (F. Vigilant)
891	!!-----------------------------------------------------------------------
892
893	!! * Modules used
894	!! * Arguments
895	INTEGER, INTENT(IN) :: &
896	& kumadt ! Output unit
897
898	!! * Local declarations
899	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
900	& zsshb_tlin, & ! Tangent input: before SSH
901	& zsshn_tlin, & ! Tangent input: before SSH
902	& zssha_tlin, & ! Tangent input: before SSH
903	& zsshb_tlout, & ! Tangent output: before SSH
904	& zsshn_tlout, & ! Tangent output: before SSH
905	& zsshb_adin, & ! Adjoint input: before SSH
906	& zsshn_adin, & ! Adjoint input: before SSH
907	& zsshb_adout, & ! Adjoint output: before SSH
908	& zsshn_adout, & ! Adjoint output: before SSH
909	& zssha_adout, & ! Adjoint output: before SSH
910	& znssh ! 2D random field for EmP
911
912	INTEGER :: &
913	& ji, & ! dummy loop indices
914	& jj, &
915	& jk
916	INTEGER, DIMENSION(jpi,jpj) :: &
917	& iseed_2d ! 2D seed for the random number generator
918	REAL(KIND=wp) :: &
919	! random field standard deviation for:
920	& zstdssh, & ! SSH
921	& zsp1, & ! scalar product involving the tangent routine
922	& zsp2, & ! scalar product involving the adjoint routine
923	& zsp1_1, & ! scalar product components
924	& zsp1_2, &
925	& zsp2_1, & ! scalar product components
926	& zsp2_2, &
927	& zsp2_3, &
928	& zsp2_4
929	CHARACTER (LEN=14) :: &
930	& cl_name
931
932	! Allocate memory
933
934	ALLOCATE( &
935	& zsshb_tlin(jpi,jpj), &
936	& zsshn_tlin(jpi,jpj), &
937	& zssha_tlin(jpi,jpj), &
938	& zsshb_tlout(jpi,jpj), &
939	& zsshn_tlout(jpi,jpj), &
940	& zsshb_adin(jpi,jpj), &
941	& zsshn_adin(jpi,jpj), &
942	& zsshb_adout(jpi,jpj), &
943	& zsshn_adout(jpi,jpj), &
944	& zssha_adout(jpi,jpj), &
945	& znssh(jpi,jpj) &
946	& )
947
948
949	! Initialize constants
950
951	zsshb_tlin(:,:) = 0.0_wp
952	zsshn_tlin(:,:) = 0.0_wp
953	zssha_tlin(:,:) = 0.0_wp
954
955	zsshb_tlout(:,:) = 0.0_wp
956	zsshn_tlout(:,:) = 0.0_wp
957
958	zsshb_adout(:,:) = 0.0_wp
959	zsshn_adout(:,:) = 0.0_wp
960	zssha_adout(:,:) = 0.0_wp
961
962	zsshb_adin(:,:) = 0.0_wp
963	zsshn_adin(:,:) = 0.0_wp
964
965	sshb_tl(:,:) = 0.0_wp
966	sshn_tl(:,:) = 0.0_wp
967	ssha_tl(:,:) = 0.0_wp
968
969	sshb_ad(:,:) = 0.0_wp
970	sshn_ad(:,:) = 0.0_wp
971	ssha_ad(:,:) = 0.0_wp
972
973	!=============================================================
974	! dx = ( sshb_tl, sshn_tl, ssha_tl ) and dy = ( ssb_tl, sshn_tl )
975	!=============================================================
976
977	!--------------------------------------------------------------------
978	! Initialize the tangent input with random noise: dx
979	!--------------------------------------------------------------------
980
981	DO jj = 1, jpj
982	DO ji = 1, jpi
983	iseed_2d(ji,jj) = - ( 785483 + &
984	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
985	END DO
986	END DO
987	CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh )
988
989	DO jj = nldj, nlej
990	DO ji = nldi, nlei
991	zsshb_tlin(ji,jj) = znssh(ji,jj)
992	END DO
993	END DO
994
995	DO jj = 1, jpj
996	DO ji = 1, jpi
997	iseed_2d(ji,jj) = - ( 358606 + &
998	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
999	END DO
1000	END DO
1001	CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh )
1002
1003	DO jj = nldj, nlej
1004	DO ji = nldi, nlei
1005	zsshn_tlin(ji,jj) = znssh(ji,jj)
1006	END DO
1007	END DO
1008
1009	DO jj = 1, jpj
1010	DO ji = 1, jpi
1011	iseed_2d(ji,jj) = - ( 596035 + &
1012	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
1013	END DO
1014	END DO
1015	CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh )
1016
1017	DO jj = nldj, nlej
1018	DO ji = nldi, nlei
1019	zssha_tlin(ji,jj) = znssh(ji,jj)
1020	END DO
1021	END DO
1022
1023	!--------------------------------------------------------------------
1024	! Call the tangent routine: dy = L dx
1025	!--------------------------------------------------------------------
1026
1027	sshb_tl(:,:) = zsshb_tlin(:,:)
1028	sshn_tl(:,:) = zsshn_tlin(:,:)
1029	ssha_tl(:,:) = zssha_tlin(:,:)
1030
1031	CALL ssh_nxt_tan( nit000+1 )
1032
1033	zsshb_tlout(:,: ) = sshb_tl(:,:)
1034	zsshn_tlout(:,: ) = sshn_tl(:,:)
1035	!--------------------------------------------------------------------
1036	! Initialize the adjoint variables: dy^* = W dy
1037	!--------------------------------------------------------------------
1038
1039	DO jj = nldj, nlej
1040	DO ji = nldi, nlei
1041	zsshb_adin(ji,jj) = zsshb_tlout(ji,jj) &
1042	& * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
1043	& * tmask(ji,jj,1)
1044	zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
1045	& * e1t(ji,jj) * e2t(ji,jj) * wesp_ssh &
1046	& * tmask(ji,jj,1)
1047	END DO
1048	END DO
1049
1050	!--------------------------------------------------------------------
1051	! Compute the scalar product: ( L dx )^T W dy
1052	!--------------------------------------------------------------------
1053	zsp1_1 = DOT_PRODUCT( zsshb_tlout, zsshb_adin )
1054	zsp1_2 = DOT_PRODUCT( zsshn_tlout, zsshn_adin )
1055
1056	zsp1 = zsp1_1 + zsp1_2
1057	!--------------------------------------------------------------------
1058	! Call the adjoint routine: dx^* = L^T dy^*
1059	!--------------------------------------------------------------------
1060
1061	sshb_ad(:,:) = zsshb_adin(:,:)
1062	sshn_ad(:,:) = zsshn_adin(:,:)
1063
1064	CALL ssh_nxt_adj( nit000+1 )
1065
1066	zsshb_adout(:,:) = sshb_ad(:,:)
1067	zsshn_adout(:,:) = sshn_ad(:,:)
1068	zssha_adout(:,:) = ssha_ad(:,:)
1069
1070	!--------------------------------------------------------------------
1071	! Compute the scalar product: dx^T L^T W dy
1072	!--------------------------------------------------------------------
1073
1074	zsp2_1 = DOT_PRODUCT( zsshb_tlin, zsshb_adout )
1075	zsp2_2 = DOT_PRODUCT( zsshn_tlin, zsshn_adout )
1076	zsp2_3 = DOT_PRODUCT( zssha_tlin, zssha_adout )
1077
1078	zsp2 = zsp2_1 + zsp2_2 + zsp2_3
1079
1080	! Compare the scalar products
1081	! 14 char:'12345678901234'
1082	cl_name = 'sshnxt_adj '
1083	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1084
1085	END SUBROUTINE ssh_nxt_adj_tst
1086
1087	!!======================================================================
1088	#endif
1089
1090	END MODULE sshwzv_tam

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

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

Context Navigation

source: branches/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/DYN/sshwzv_tam.F90 @ 5226

Download in other formats: