1 | MODULE zpshde_tam |
---|
2 | #ifdef key_tam |
---|
3 | !!============================================================================== |
---|
4 | !! *** MODULE zpshde_tam *** |
---|
5 | !! z-coordinate - partial step : Horizontal Derivative |
---|
6 | !! Tangent and Adjoint Module |
---|
7 | !!============================================================================== |
---|
8 | |
---|
9 | !!---------------------------------------------------------------------- |
---|
10 | !! zps_hde : Horizontal DErivative of T, S and rd at the last |
---|
11 | !! ocean level (Z-coord. with Partial Steps) |
---|
12 | !!---------------------------------------------------------------------- |
---|
13 | !! * Modules used |
---|
14 | USE par_kind , ONLY: & ! Precision variables |
---|
15 | & wp |
---|
16 | USE par_oce , ONLY: & ! Ocean space and time domain variables |
---|
17 | & jpi, & |
---|
18 | & jpj, & |
---|
19 | & jpk, & |
---|
20 | & jpim1, & |
---|
21 | & jpjm1, & |
---|
22 | & jpiglo |
---|
23 | USE oce , ONLY: & ! ocean variables |
---|
24 | & tn, & |
---|
25 | & sn |
---|
26 | USE dom_oce , ONLY: & ! ocean space domain variables |
---|
27 | & umask, & |
---|
28 | & vmask, & |
---|
29 | & e2u, & |
---|
30 | & e1u, & |
---|
31 | & e2v, & |
---|
32 | & e1v, & |
---|
33 | # if defined key_zco |
---|
34 | & gdept_0, & |
---|
35 | & e3t_0, & |
---|
36 | & e3w_0, & |
---|
37 | # else |
---|
38 | & gdept, & |
---|
39 | & e3u, & |
---|
40 | & e3v, & |
---|
41 | & e3w, & |
---|
42 | # endif |
---|
43 | & mbathy, & |
---|
44 | & mig, & |
---|
45 | & mjg, & |
---|
46 | & nldi, & |
---|
47 | & nldj, & |
---|
48 | & nlei, & |
---|
49 | & nlej |
---|
50 | USE in_out_manager, ONLY: & ! I/O manager |
---|
51 | & lwp, & |
---|
52 | & numout, & |
---|
53 | & nit000, & |
---|
54 | & nitend |
---|
55 | USE eosbn2_tam , ONLY: & ! ocean equation of state |
---|
56 | & eos_tan, & |
---|
57 | & eos_adj |
---|
58 | USE lbclnk , ONLY: & ! lateral boundary conditions (or mpp link) |
---|
59 | & lbc_lnk |
---|
60 | USE lbclnk_tam , ONLY: & ! lateral boundary conditions (or mpp link) |
---|
61 | & lbc_lnk_adj |
---|
62 | USE gridrandom , ONLY: & ! Random Gaussian noise on grids |
---|
63 | & grid_random |
---|
64 | USE dotprodfld, ONLY: & ! Computes dot product for 3D and 2D fields |
---|
65 | & dot_product |
---|
66 | USE tstool_tam , ONLY: & |
---|
67 | & prntst_adj, & ! |
---|
68 | & prntst_tlm, & ! |
---|
69 | & stdt, & ! stdev for temperature |
---|
70 | & stds, & ! salinity |
---|
71 | & stdr ! |
---|
72 | |
---|
73 | IMPLICIT NONE |
---|
74 | PRIVATE |
---|
75 | |
---|
76 | !! * Routine accessibility |
---|
77 | PUBLIC zps_hde_tan ! routine called by step_tam.F90 |
---|
78 | PUBLIC zps_hde_adj ! routine called by step_tam.F90 |
---|
79 | PUBLIC zps_hde_adj_tst ! routine called by tst.F90 |
---|
80 | PUBLIC zps_hde_tlm_tst ! routine called by tamtst.F90 |
---|
81 | |
---|
82 | !! * module variables |
---|
83 | INTEGER, DIMENSION(jpi,jpj) :: & |
---|
84 | & mbatu, mbatv ! bottom ocean level index at U- and V-points |
---|
85 | |
---|
86 | !! * Substitutions |
---|
87 | # include "domzgr_substitute.h90" |
---|
88 | # include "vectopt_loop_substitute.h90" |
---|
89 | |
---|
90 | CONTAINS |
---|
91 | |
---|
92 | SUBROUTINE zps_hde_tan ( kt, ptem, psal, & |
---|
93 | & ptem_tl, psal_tl, prd_tl , & |
---|
94 | pgtu_tl, pgsu_tl, pgru_tl, & |
---|
95 | pgtv_tl, pgsv_tl, pgrv_tl ) |
---|
96 | !!---------------------------------------------------------------------- |
---|
97 | !! *** ROUTINE zps_hde_tan *** |
---|
98 | !! |
---|
99 | !! ** Purpose of the direct routine: |
---|
100 | !! Compute the horizontal derivative of T, S and rd |
---|
101 | !! at u- and v-points with a linear interpolation for z-coordinate |
---|
102 | !! with partial steps. |
---|
103 | !! |
---|
104 | !! ** Method of the direct routine: |
---|
105 | !! In z-coord with partial steps, scale factors on last |
---|
106 | !! levels are different for each grid point, so that T, S and rd |
---|
107 | !! points are not at the same depth as in z-coord. To have horizontal |
---|
108 | !! gradients again, we interpolate T and S at the good depth : |
---|
109 | !! Linear interpolation of T, S |
---|
110 | !! Computation of di(tb) and dj(tb) by vertical interpolation: |
---|
111 | !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ |
---|
112 | !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ |
---|
113 | !! This formulation computes the two cases: |
---|
114 | !! CASE 1 CASE 2 |
---|
115 | !! k-1 ___ ___________ k-1 ___ ___________ |
---|
116 | !! Ti T~ T~ Ti+1 |
---|
117 | !! _____ _____ |
---|
118 | !! k | |Ti+1 k Ti | | |
---|
119 | !! | |____ ____| | |
---|
120 | !! ___ | | | ___ | | | |
---|
121 | !! |
---|
122 | !! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then |
---|
123 | !! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) |
---|
124 | !! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) ) |
---|
125 | !! or |
---|
126 | !! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then |
---|
127 | !! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) |
---|
128 | !! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) |
---|
129 | !! Idem for di(s) and dj(s) |
---|
130 | !! |
---|
131 | !! For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at |
---|
132 | !! the good depth zh from interpolated T and S for the different |
---|
133 | !! formulation of the equation of state (eos). |
---|
134 | !! Gradient formulation for rho : |
---|
135 | !! di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~ |
---|
136 | !! |
---|
137 | !! ** Action : - pgtu, pgsu, pgru: horizontal gradient of T, S |
---|
138 | !! and rd at U-points |
---|
139 | !! - pgtv, pgsv, pgrv: horizontal gradient of T, S |
---|
140 | !! and rd at V-points |
---|
141 | !! |
---|
142 | !! History of the direct routine: |
---|
143 | !! 8.5 ! 02-04 (A. Bozec) Original code |
---|
144 | !! 8.5 ! 02-08 (G. Madec E. Durand) Optimization and Free form |
---|
145 | !! History of the TAM routine: |
---|
146 | !! 9.0 ! 08-06 (A. Vidard) Skeleton |
---|
147 | !! ! 08-06 (A. Vidard) tangent of the 02-08 version |
---|
148 | !!---------------------------------------------------------------------- |
---|
149 | !! * Arguments |
---|
150 | INTEGER, INTENT( in ) :: kt ! ocean time-step index |
---|
151 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & |
---|
152 | ptem, psal ! 3D T, S and rd direct fields |
---|
153 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & |
---|
154 | ptem_tl, psal_tl, prd_tl ! 3D T, S and rd tangent fields |
---|
155 | REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) :: & |
---|
156 | pgtu_tl, pgsu_tl, pgru_tl, & ! horizontal grad. of T, S and rd at u- |
---|
157 | pgtv_tl, pgsv_tl, pgrv_tl ! and v-points of the partial step level |
---|
158 | |
---|
159 | |
---|
160 | !! * Local declarations |
---|
161 | INTEGER :: ji, jj, & ! Dummy loop indices |
---|
162 | iku,ikv ! partial step level at u- and v-points |
---|
163 | REAL(wp), DIMENSION(jpi,jpj) :: & |
---|
164 | ztmpi, ztmpj, & |
---|
165 | zti, ztj, zsi, zsj, & ! interpolated value of T, S |
---|
166 | zri, zrj, & ! and rd |
---|
167 | ztitl, ztjtl, & ! interpolated value of Ttl |
---|
168 | zsitl, zsjtl, & !, Stl |
---|
169 | zritl, zrjtl, & ! and rdtl |
---|
170 | zhgi, zhgj ! depth of interpolation for eos2d |
---|
171 | REAL(wp) :: & |
---|
172 | ze3wu, ze3wv, & ! temporary scalars |
---|
173 | zmaxu1, zmaxu2, & ! " " |
---|
174 | zmaxv1, zmaxv2 ! " " |
---|
175 | |
---|
176 | ! Initialization (first time-step only): compute mbatu and mbatv |
---|
177 | IF( kt == nit000 ) THEN |
---|
178 | mbatu(:,:) = 0 |
---|
179 | mbatv(:,:) = 0 |
---|
180 | DO jj = 1, jpjm1 |
---|
181 | DO ji = 1, fs_jpim1 ! vector opt. |
---|
182 | mbatu(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1, 2 ) |
---|
183 | mbatv(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1, 2 ) |
---|
184 | END DO |
---|
185 | END DO |
---|
186 | ztmpi(:,:) = FLOAT( mbatu(:,:) ) |
---|
187 | ztmpj(:,:) = FLOAT( mbatv(:,:) ) |
---|
188 | ! lateral boundary conditions: T-point, sign unchanged |
---|
189 | CALL lbc_lnk( ztmpi , 'U', 1. ) |
---|
190 | CALL lbc_lnk( ztmpj , 'V', 1. ) |
---|
191 | mbatu(:,:) = MAX( INT( ztmpi(:,:) ), 2 ) |
---|
192 | mbatv(:,:) = MAX( INT( ztmpj(:,:) ), 2 ) |
---|
193 | ENDIF |
---|
194 | |
---|
195 | |
---|
196 | ! Interpolation of T and S at the last ocean level |
---|
197 | # if defined key_vectopt_loop |
---|
198 | jj = 1 |
---|
199 | DO ji = 1, jpij-jpi ! vector opt. (forced unrolled) |
---|
200 | # else |
---|
201 | DO jj = 1, jpjm1 |
---|
202 | DO ji = 1, jpim1 |
---|
203 | # endif |
---|
204 | ! last level |
---|
205 | iku = mbatu(ji,jj) |
---|
206 | ikv = mbatv(ji,jj) |
---|
207 | |
---|
208 | ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) |
---|
209 | ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) |
---|
210 | zmaxu1 = ze3wu / fse3w(ji+1,jj ,iku) |
---|
211 | zmaxu2 = -ze3wu / fse3w(ji ,jj ,iku) |
---|
212 | zmaxv1 = ze3wv / fse3w(ji ,jj+1,ikv) |
---|
213 | zmaxv2 = -ze3wv / fse3w(ji ,jj ,ikv) |
---|
214 | |
---|
215 | ! i- direction |
---|
216 | |
---|
217 | IF( ze3wu >= 0. ) THEN ! case 1 |
---|
218 | ! interpolated values of T and S |
---|
219 | zti(ji,jj) = ptem(ji+1,jj,iku) & |
---|
220 | & + zmaxu1 * ( ptem(ji+1,jj,iku-1) - ptem(ji+1,jj,iku) ) |
---|
221 | zsi(ji,jj) = psal(ji+1,jj,iku) & |
---|
222 | & + zmaxu1 * ( psal(ji+1,jj,iku-1) - psal(ji+1,jj,iku) ) |
---|
223 | ztitl(ji,jj) = ptem_tl(ji+1,jj,iku) & |
---|
224 | & + zmaxu1 * ( ptem_tl(ji+1,jj,iku-1) - ptem_tl(ji+1,jj,iku) ) |
---|
225 | zsitl(ji,jj) = psal_tl(ji+1,jj,iku) & |
---|
226 | & + zmaxu1 * ( psal_tl(ji+1,jj,iku-1) - psal_tl(ji+1,jj,iku) ) |
---|
227 | ! depth of the partial step level |
---|
228 | zhgi(ji,jj) = fsdept(ji,jj,iku) |
---|
229 | ! gradient of T and S |
---|
230 | pgtu_tl(ji,jj) = umask(ji,jj,1) * ( ztitl(ji,jj) - ptem_tl(ji,jj,iku) ) |
---|
231 | pgsu_tl(ji,jj) = umask(ji,jj,1) * ( zsitl(ji,jj) - psal_tl(ji,jj,iku) ) |
---|
232 | ELSE ! case 2 |
---|
233 | ! interpolated values of T and S |
---|
234 | zti(ji,jj) = ptem(ji,jj,iku) & |
---|
235 | & + zmaxu2 * ( ptem(ji,jj,iku-1) - ptem(ji,jj,iku) ) |
---|
236 | zsi(ji,jj) = psal(ji,jj,iku) & |
---|
237 | & + zmaxu2 * ( psal(ji,jj,iku-1) - psal(ji,jj,iku) ) |
---|
238 | ! interpolated values of T and S |
---|
239 | ztitl(ji,jj) = ptem_tl(ji,jj,iku) & |
---|
240 | & + zmaxu2 * ( ptem_tl(ji,jj,iku-1) - ptem_tl(ji,jj,iku) ) |
---|
241 | zsitl(ji,jj) = psal_tl(ji,jj,iku) & |
---|
242 | & + zmaxu2 * ( psal_tl(ji,jj,iku-1) - psal_tl(ji,jj,iku) ) |
---|
243 | ! depth of the partial step level |
---|
244 | zhgi(ji,jj) = fsdept(ji+1,jj,iku) |
---|
245 | ! gradient of T and S |
---|
246 | pgtu_tl(ji,jj) = umask(ji,jj,1) * ( ptem_tl(ji+1,jj,iku) - ztitl (ji,jj) ) |
---|
247 | pgsu_tl(ji,jj) = umask(ji,jj,1) * ( psal_tl(ji+1,jj,iku) - zsitl (ji,jj) ) |
---|
248 | ENDIF |
---|
249 | |
---|
250 | ! j- direction |
---|
251 | |
---|
252 | IF( ze3wv >= 0. ) THEN ! case 1 |
---|
253 | ! interpolated values of direct T and S |
---|
254 | ztj(ji,jj) = ptem(ji,jj+1,ikv) & |
---|
255 | & + zmaxv1 * ( ptem(ji,jj+1,ikv-1) - ptem(ji,jj+1,ikv) ) |
---|
256 | zsj(ji,jj) = psal(ji,jj+1,ikv) & |
---|
257 | & + zmaxv1 * ( psal(ji,jj+1,ikv-1) - psal(ji,jj+1,ikv) ) |
---|
258 | ! interpolated values of tangent T and S |
---|
259 | ztjtl(ji,jj) = ptem_tl(ji,jj+1,ikv) & |
---|
260 | & + zmaxv1 * ( ptem_tl(ji,jj+1,ikv-1) - ptem_tl(ji,jj+1,ikv) ) |
---|
261 | zsjtl(ji,jj) = psal_tl(ji,jj+1,ikv) & |
---|
262 | & + zmaxv1 * ( psal_tl(ji,jj+1,ikv-1) - psal_tl(ji,jj+1,ikv) ) |
---|
263 | ! depth of the partial step level |
---|
264 | zhgj(ji,jj) = fsdept(ji,jj,ikv) |
---|
265 | ! gradient of T and S |
---|
266 | pgtv_tl(ji,jj) = vmask(ji,jj,1) * ( ztjtl(ji,jj) - ptem_tl(ji,jj,ikv) ) |
---|
267 | pgsv_tl(ji,jj) = vmask(ji,jj,1) * ( zsjtl(ji,jj) - psal_tl(ji,jj,ikv) ) |
---|
268 | |
---|
269 | ELSE ! case 2 |
---|
270 | ! interpolated values of T and S |
---|
271 | ztj(ji,jj) = ptem(ji,jj,ikv) & |
---|
272 | & + zmaxv2 * ( ptem(ji,jj,ikv-1) - ptem(ji,jj,ikv) ) |
---|
273 | zsj(ji,jj) = psal(ji,jj,ikv) & |
---|
274 | & + zmaxv2 * ( psal(ji,jj,ikv-1) - psal(ji,jj,ikv) ) |
---|
275 | ! interpolated values of T and S |
---|
276 | ztjtl(ji,jj) = ptem_tl(ji,jj,ikv) & |
---|
277 | & + zmaxv2 * ( ptem_tl(ji,jj,ikv-1) - ptem_tl(ji,jj,ikv) ) |
---|
278 | zsjtl(ji,jj) = psal_tl(ji,jj,ikv) & |
---|
279 | & + zmaxv2 * ( psal_tl(ji,jj,ikv-1) - psal_tl(ji,jj,ikv) ) |
---|
280 | ! depth of the partial step level |
---|
281 | zhgj(ji,jj) = fsdept(ji,jj+1,ikv) |
---|
282 | ! gradient of T and S |
---|
283 | pgtv_tl(ji,jj) = vmask(ji,jj,1) * ( ptem_tl(ji,jj+1,ikv) - ztjtl(ji,jj) ) |
---|
284 | pgsv_tl(ji,jj) = vmask(ji,jj,1) * ( psal_tl(ji,jj+1,ikv) - zsjtl(ji,jj) ) |
---|
285 | ENDIF |
---|
286 | # if ! defined key_vectopt_loop |
---|
287 | END DO |
---|
288 | # endif |
---|
289 | END DO |
---|
290 | |
---|
291 | ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial |
---|
292 | ! step and store it in zri, zrj for each case |
---|
293 | CALL eos_tan( zti, zsi, zhgi, ztitl, zsitl, zritl ) |
---|
294 | CALL eos_tan( ztj, zsj, zhgj, ztjtl, zsjtl, zrjtl ) |
---|
295 | |
---|
296 | |
---|
297 | ! Gradient of density at the last level |
---|
298 | # if defined key_vectopt_loop |
---|
299 | jj = 1 |
---|
300 | DO ji = 1, jpij-jpi ! vector opt. (forced unrolled) |
---|
301 | # else |
---|
302 | DO jj = 1, jpjm1 |
---|
303 | DO ji = 1, jpim1 |
---|
304 | # endif |
---|
305 | iku = mbatu(ji,jj) |
---|
306 | ikv = mbatv(ji,jj) |
---|
307 | ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) |
---|
308 | ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) |
---|
309 | IF( ze3wu >= 0. ) THEN ! i-direction: case 1 |
---|
310 | pgru_tl(ji,jj) = umask(ji,jj,1) * ( zritl(ji,jj) - prd_tl(ji,jj,iku) ) |
---|
311 | ELSE ! i-direction: case 2 |
---|
312 | pgru_tl(ji,jj) = umask(ji,jj,1) * ( prd_tl(ji+1,jj,iku) - zritl(ji,jj) ) |
---|
313 | ENDIF |
---|
314 | IF( ze3wv >= 0. ) THEN ! j-direction: case 1 |
---|
315 | pgrv_tl(ji,jj) = vmask(ji,jj,1) * ( zrjtl(ji,jj) - prd_tl(ji,jj,ikv) ) |
---|
316 | ELSE ! j-direction: case 2 |
---|
317 | pgrv_tl(ji,jj) = vmask(ji,jj,1) * ( prd_tl(ji,jj+1,ikv) - zrjtl(ji,jj) ) |
---|
318 | ENDIF |
---|
319 | # if ! defined key_vectopt_loop |
---|
320 | END DO |
---|
321 | # endif |
---|
322 | END DO |
---|
323 | |
---|
324 | ! Lateral boundary conditions on each gradient |
---|
325 | CALL lbc_lnk( pgtu_tl , 'U', -1. ) ; CALL lbc_lnk( pgtv_tl , 'V', -1. ) |
---|
326 | CALL lbc_lnk( pgsu_tl , 'U', -1. ) ; CALL lbc_lnk( pgsv_tl , 'V', -1. ) |
---|
327 | CALL lbc_lnk( pgru_tl , 'U', -1. ) ; CALL lbc_lnk( pgrv_tl , 'V', -1. ) |
---|
328 | |
---|
329 | END SUBROUTINE zps_hde_tan |
---|
330 | |
---|
331 | SUBROUTINE zps_hde_adj ( kt, ptem, psal, & |
---|
332 | & ptem_ad, psal_ad, prd_ad , & |
---|
333 | pgtu_ad, pgsu_ad, pgru_ad, & |
---|
334 | pgtv_ad, pgsv_ad, pgrv_ad ) |
---|
335 | !!---------------------------------------------------------------------- |
---|
336 | !! *** ROUTINE zps_hde_adj *** |
---|
337 | !! |
---|
338 | !! ** Purpose of the direct routine: |
---|
339 | !! Compute the horizontal derivative of T, S and rd |
---|
340 | !! at u- and v-points with a linear interpolation for z-coordinate |
---|
341 | !! with partial steps. |
---|
342 | !! |
---|
343 | !! ** Method of the direct routine: |
---|
344 | !! In z-coord with partial steps, scale factors on last |
---|
345 | !! levels are different for each grid point, so that T, S and rd |
---|
346 | !! points are not at the same depth as in z-coord. To have horizontal |
---|
347 | !! gradients again, we interpolate T and S at the good depth : |
---|
348 | !! Linear interpolation of T, S |
---|
349 | !! Computation of di(tb) and dj(tb) by vertical interpolation: |
---|
350 | !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ |
---|
351 | !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ |
---|
352 | !! This formulation computes the two cases: |
---|
353 | !! CASE 1 CASE 2 |
---|
354 | !! k-1 ___ ___________ k-1 ___ ___________ |
---|
355 | !! Ti T~ T~ Ti+1 |
---|
356 | !! _____ _____ |
---|
357 | !! k | |Ti+1 k Ti | | |
---|
358 | !! | |____ ____| | |
---|
359 | !! ___ | | | ___ | | | |
---|
360 | !! |
---|
361 | !! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then |
---|
362 | !! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) |
---|
363 | !! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) ) |
---|
364 | !! or |
---|
365 | !! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then |
---|
366 | !! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) |
---|
367 | !! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) |
---|
368 | !! Idem for di(s) and dj(s) |
---|
369 | !! |
---|
370 | !! For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at |
---|
371 | !! the good depth zh from interpolated T and S for the different |
---|
372 | !! formulation of the equation of state (eos). |
---|
373 | !! Gradient formulation for rho : |
---|
374 | !! di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~ |
---|
375 | !! |
---|
376 | !! ** Action : - pgtu, pgsu, pgru: horizontal gradient of T, S |
---|
377 | !! and rd at U-points |
---|
378 | !! - pgtv, pgsv, pgrv: horizontal gradient of T, S |
---|
379 | !! and rd at V-points |
---|
380 | !! |
---|
381 | !! History of the direct routine: |
---|
382 | !! 8.5 ! 02-04 (A. Bozec) Original code |
---|
383 | !! 8.5 ! 02-08 (G. Madec E. Durand) Optimization and Free form |
---|
384 | !! History of the TAM routine: |
---|
385 | !! 9.0 ! 08-06 (A. Vidard) Skeleton |
---|
386 | !! ! 08-08 (A. Vidard) adjoint of the 02-08 version |
---|
387 | !!---------------------------------------------------------------------- |
---|
388 | !! * Arguments |
---|
389 | INTEGER, INTENT( in ) :: kt ! ocean time-step index |
---|
390 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & |
---|
391 | ptem, psal ! 3D T, S and rd direct fields |
---|
392 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & |
---|
393 | ptem_ad, psal_ad, prd_ad ! 3D T, S and rd tangent fields |
---|
394 | REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: & |
---|
395 | pgtu_ad, pgsu_ad, pgru_ad, & ! horizontal grad. of T, S and rd at u- |
---|
396 | pgtv_ad, pgsv_ad, pgrv_ad ! and v-points of the partial step level |
---|
397 | |
---|
398 | |
---|
399 | !! * Local declarations |
---|
400 | INTEGER :: ji, jj, & ! Dummy loop indices |
---|
401 | iku,ikv ! partial step level at u- and v-points |
---|
402 | REAL(wp), DIMENSION(jpi,jpj) :: & |
---|
403 | ztmpi, ztmpj, & |
---|
404 | zti, ztj, zsi, zsj, & ! interpolated value of T, S |
---|
405 | zri, zrj, & ! and rd |
---|
406 | ztiad, ztjad, zsiad, zsjad, & ! interpolated value of Tad, Sad |
---|
407 | zriad, zrjad, & ! and rdad |
---|
408 | zhgi, zhgj ! depth of interpolation for eos2d |
---|
409 | REAL(wp) :: & |
---|
410 | ze3wu, ze3wv, & ! temporary scalars |
---|
411 | zmaxu1, zmaxu2, & ! " " |
---|
412 | zmaxv1, zmaxv2 ! " " |
---|
413 | |
---|
414 | ! Initialization (first time-step only): compute mbatu and mbatv |
---|
415 | IF( kt == nitend ) THEN |
---|
416 | mbatu(:,:) = 0 |
---|
417 | mbatv(:,:) = 0 |
---|
418 | DO jj = 1, jpjm1 |
---|
419 | DO ji = 1, fs_jpim1 ! vector opt. |
---|
420 | mbatu(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1, 2 ) |
---|
421 | mbatv(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1, 2 ) |
---|
422 | END DO |
---|
423 | END DO |
---|
424 | ztmpi(:,:) = FLOAT( mbatu(:,:) ) |
---|
425 | ztmpj(:,:) = FLOAT( mbatv(:,:) ) |
---|
426 | ! lateral boundary conditions: T-point, sign unchanged |
---|
427 | CALL lbc_lnk( ztmpi , 'U', 1. ) |
---|
428 | CALL lbc_lnk( ztmpj , 'V', 1. ) |
---|
429 | mbatu(:,:) = MAX( INT( ztmpi(:,:) ), 2 ) |
---|
430 | mbatv(:,:) = MAX( INT( ztmpj(:,:) ), 2 ) |
---|
431 | ENDIF |
---|
432 | ! 1. Direct model recomputation |
---|
433 | ! Interpolation of T and S at the last ocean level |
---|
434 | # if defined key_vectopt_loop |
---|
435 | jj = 1 |
---|
436 | DO ji = 1, jpij-jpi ! vector opt. (forced unrolled) |
---|
437 | # else |
---|
438 | DO jj = 1, jpjm1 |
---|
439 | DO ji = 1, jpim1 |
---|
440 | # endif |
---|
441 | ! last level |
---|
442 | iku = mbatu(ji,jj) |
---|
443 | ikv = mbatv(ji,jj) |
---|
444 | |
---|
445 | ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) |
---|
446 | ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) |
---|
447 | zmaxu1 = ze3wu / fse3w(ji+1,jj ,iku) |
---|
448 | zmaxu2 = -ze3wu / fse3w(ji ,jj ,iku) |
---|
449 | zmaxv1 = ze3wv / fse3w(ji ,jj+1,ikv) |
---|
450 | zmaxv2 = -ze3wv / fse3w(ji ,jj ,ikv) |
---|
451 | |
---|
452 | ! i- direction |
---|
453 | |
---|
454 | IF( ze3wu >= 0. ) THEN ! case 1 |
---|
455 | ! interpolated values of T and S |
---|
456 | zti(ji,jj) = ptem(ji+1,jj,iku) & |
---|
457 | & + zmaxu1 * ( ptem(ji+1,jj,iku-1) - ptem(ji+1,jj,iku) ) |
---|
458 | zsi(ji,jj) = psal(ji+1,jj,iku) & |
---|
459 | & + zmaxu1 * ( psal(ji+1,jj,iku-1) - psal(ji+1,jj,iku) ) |
---|
460 | ! depth of the partial step level |
---|
461 | zhgi(ji,jj) = fsdept(ji,jj,iku) |
---|
462 | ELSE ! case 2 |
---|
463 | ! interpolated values of T and S |
---|
464 | zti(ji,jj) = ptem(ji,jj,iku) & |
---|
465 | & + zmaxu2 * ( ptem(ji,jj,iku-1) - ptem(ji,jj,iku) ) |
---|
466 | zsi(ji,jj) = psal(ji,jj,iku) + zmaxu2 * ( psal(ji,jj,iku-1) - psal(ji,jj,iku) ) |
---|
467 | ! depth of the partial step level |
---|
468 | zhgi(ji,jj) = fsdept(ji+1,jj,iku) |
---|
469 | ENDIF |
---|
470 | |
---|
471 | ! j- direction |
---|
472 | |
---|
473 | IF( ze3wv >= 0. ) THEN ! case 1 |
---|
474 | ! interpolated values of T and S |
---|
475 | ztj(ji,jj) = ptem(ji,jj+1,ikv) & |
---|
476 | & + zmaxv1 * ( ptem(ji,jj+1,ikv-1) - ptem(ji,jj+1,ikv) ) |
---|
477 | zsj(ji,jj) = psal(ji,jj+1,ikv) & |
---|
478 | & + zmaxv1 * ( psal(ji,jj+1,ikv-1) - psal(ji,jj+1,ikv) ) |
---|
479 | ! depth of the partial step level |
---|
480 | zhgj(ji,jj) = fsdept(ji,jj,ikv) |
---|
481 | ELSE ! case 2 |
---|
482 | ! interpolated values of T and S |
---|
483 | ztj(ji,jj) = ptem(ji,jj,ikv) & |
---|
484 | & + zmaxv2 * ( ptem(ji,jj,ikv-1) - ptem(ji,jj,ikv) ) |
---|
485 | zsj(ji,jj) = psal(ji,jj,ikv) & |
---|
486 | & + zmaxv2 * ( psal(ji,jj,ikv-1) - psal(ji,jj,ikv) ) |
---|
487 | ! depth of the partial step level |
---|
488 | zhgj(ji,jj) = fsdept(ji,jj+1,ikv) |
---|
489 | ENDIF |
---|
490 | # if ! defined key_vectopt_loop |
---|
491 | END DO |
---|
492 | # endif |
---|
493 | END DO |
---|
494 | |
---|
495 | !2. Adjoint model counterpart |
---|
496 | ztiad = 0.0_wp ; ztjad = 0.0_wp ; zsiad = 0.0_wp |
---|
497 | zsjad = 0.0_wp ; zriad = 0.0_wp ; zrjad = 0.0_wp |
---|
498 | |
---|
499 | ! Lateral boundary conditions on each gradient |
---|
500 | CALL lbc_lnk_adj( pgtu_ad , 'U', -1. ) ; CALL lbc_lnk_adj( pgtv_ad , 'V', -1. ) |
---|
501 | CALL lbc_lnk_adj( pgsu_ad , 'U', -1. ) ; CALL lbc_lnk_adj( pgsv_ad , 'V', -1. ) |
---|
502 | CALL lbc_lnk_adj( pgru_ad , 'U', -1. ) ; CALL lbc_lnk_adj( pgrv_ad , 'V', -1. ) |
---|
503 | |
---|
504 | ! Gradient of density at the last level |
---|
505 | # if defined key_vectopt_loop |
---|
506 | jj = 1 |
---|
507 | DO ji = jpij-jpi, -1 ! vector opt. (forced unrolled) |
---|
508 | # else |
---|
509 | DO jj = jpjm1, 1, -1 |
---|
510 | DO ji = jpim1, 1, -1 |
---|
511 | # endif |
---|
512 | iku = mbatu(ji,jj) |
---|
513 | ikv = mbatv(ji,jj) |
---|
514 | ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) |
---|
515 | ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) |
---|
516 | IF( ze3wv >= 0. ) THEN ! j-direction: case 1 |
---|
517 | zrjad(ji,jj) = zrjad(ji,jj) & |
---|
518 | & + pgrv_ad(ji,jj) * vmask(ji,jj,1) |
---|
519 | prd_ad(ji,jj,ikv) = prd_ad(ji,jj,ikv) & |
---|
520 | & - pgrv_ad(ji,jj) * vmask(ji,jj,1) |
---|
521 | pgrv_ad(ji,jj) = 0.0_wp |
---|
522 | ELSE ! j-direction: case 2 |
---|
523 | prd_ad(ji,jj+1,ikv) = prd_ad(ji,jj+1,ikv) & |
---|
524 | & + pgrv_ad(ji,jj) * vmask(ji,jj,1) |
---|
525 | zrjad(ji,jj) = zrjad(ji,jj) & |
---|
526 | & - pgrv_ad(ji,jj) * vmask(ji,jj,1) |
---|
527 | pgrv_ad(ji,jj) = 0.0_wp |
---|
528 | ENDIF |
---|
529 | IF( ze3wu >= 0. ) THEN ! i-direction: case 1 |
---|
530 | zriad(ji,jj) = zriad(ji,jj) & |
---|
531 | & + pgru_ad(ji,jj) * umask(ji,jj,1) |
---|
532 | prd_ad(ji,jj,iku) = prd_ad(ji,jj,iku) & |
---|
533 | & - pgru_ad(ji,jj) * umask(ji,jj,1) |
---|
534 | pgru_ad(ji,jj) = 0.0_wp |
---|
535 | ELSE ! i-direction: case 2 |
---|
536 | prd_ad(ji+1,jj,iku) = prd_ad(ji+1,jj,iku) & |
---|
537 | & + pgru_ad(ji,jj) * umask(ji,jj,1) |
---|
538 | zriad(ji,jj) = zriad(ji,jj) & |
---|
539 | & - pgru_ad(ji,jj) * umask(ji,jj,1) |
---|
540 | pgru_ad(ji,jj) = 0.0_wp |
---|
541 | ENDIF |
---|
542 | |
---|
543 | # if ! defined key_vectopt_loop |
---|
544 | END DO |
---|
545 | # endif |
---|
546 | END DO |
---|
547 | |
---|
548 | |
---|
549 | ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial |
---|
550 | ! step and store it in zri, zrj for each case |
---|
551 | CALL eos_adj( ztj, zsj, zhgj, ztjad, zsjad, zrjad ) |
---|
552 | CALL eos_adj( zti, zsi, zhgi, ztiad, zsiad, zriad ) |
---|
553 | |
---|
554 | |
---|
555 | # if defined key_vectopt_loop |
---|
556 | jj = 1 |
---|
557 | DO ji = jpij-jpi, 1, -1 ! vector opt. (forced unrolled) |
---|
558 | # else |
---|
559 | DO jj = jpjm1, 1, -1 |
---|
560 | DO ji = jpim1, 1, -1 |
---|
561 | # endif |
---|
562 | ! last level |
---|
563 | iku = mbatu(ji,jj) |
---|
564 | ikv = mbatv(ji,jj) |
---|
565 | |
---|
566 | ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) |
---|
567 | ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) |
---|
568 | zmaxu1 = ze3wu / fse3w(ji+1,jj ,iku) |
---|
569 | zmaxu2 = -ze3wu / fse3w(ji ,jj ,iku) |
---|
570 | zmaxv1 = ze3wv / fse3w(ji ,jj+1,ikv) |
---|
571 | zmaxv2 = -ze3wv / fse3w(ji ,jj ,ikv) |
---|
572 | |
---|
573 | ! j- direction |
---|
574 | |
---|
575 | IF( ze3wv >= 0. ) THEN ! case 1 |
---|
576 | ! gradient of T and S |
---|
577 | ztjad(ji,jj) = ztjad(ji,jj) & |
---|
578 | & + pgtv_ad(ji,jj) * vmask(ji,jj,1) |
---|
579 | ptem_ad(ji,jj,ikv) = ptem_ad(ji,jj,ikv) & |
---|
580 | & - pgtv_ad(ji,jj) * vmask(ji,jj,1) |
---|
581 | pgtv_ad(ji,jj) = 0.0_wp |
---|
582 | |
---|
583 | zsjad(ji,jj) = zsjad(ji,jj) & |
---|
584 | & + pgsv_ad(ji,jj) * vmask(ji,jj,1) |
---|
585 | psal_ad(ji,jj,ikv) = psal_ad(ji,jj,ikv) & |
---|
586 | & - pgsv_ad(ji,jj) * vmask(ji,jj,1) |
---|
587 | pgsv_ad(ji,jj) = 0.0_wp |
---|
588 | ! interpolated values of T and S |
---|
589 | ptem_ad(ji,jj+1,ikv) = ptem_ad(ji,jj+1,ikv) & |
---|
590 | & + ztjad(ji,jj) * (1 - zmaxv1) |
---|
591 | ptem_ad(ji,jj+1,ikv-1) = ptem_ad(ji,jj+1,ikv-1) & |
---|
592 | & + ztjad(ji,jj)* zmaxv1 |
---|
593 | ztjad(ji,jj) = 0.0_wp |
---|
594 | |
---|
595 | psal_ad(ji,jj+1,ikv) = psal_ad(ji,jj+1,ikv) & |
---|
596 | & + zsjad(ji,jj) * (1 - zmaxv1) |
---|
597 | psal_ad(ji,jj+1,ikv-1) = psal_ad(ji,jj+1,ikv-1) & |
---|
598 | & + zsjad(ji,jj) * zmaxv1 |
---|
599 | zsjad(ji,jj) = 0.0_wp |
---|
600 | ELSE ! case 2 |
---|
601 | ! gradient of T and S |
---|
602 | ptem_ad(ji,jj+1,ikv) = ptem_ad(ji,jj+1,ikv) & |
---|
603 | & + pgtv_ad(ji,jj) * vmask(ji,jj,1) |
---|
604 | ztjad(ji,jj) = ztjad(ji,jj) & |
---|
605 | & - pgtv_ad(ji,jj) * vmask(ji,jj,1) |
---|
606 | pgtv_ad(ji,jj) = 0.0_wp |
---|
607 | |
---|
608 | psal_ad(ji,jj+1,ikv) = psal_ad(ji,jj+1,ikv) & |
---|
609 | & + pgsv_ad(ji,jj) * vmask(ji,jj,1) |
---|
610 | zsjad(ji,jj) = zsjad(ji,jj) & |
---|
611 | & - pgsv_ad(ji,jj) * vmask(ji,jj,1) |
---|
612 | pgsv_ad(ji,jj) = 0.0_wp |
---|
613 | ! interpolated values of T and S |
---|
614 | ptem_ad(ji,jj,ikv) = ptem_ad(ji,jj,ikv) & |
---|
615 | & + ztjad(ji,jj) * (1 - zmaxv2) |
---|
616 | ptem_ad(ji,jj,ikv-1) = ptem_ad(ji,jj,ikv-1) & |
---|
617 | & + ztjad(ji,jj) * zmaxv2 |
---|
618 | ztjad(ji,jj) = 0.0_wp |
---|
619 | |
---|
620 | psal_ad(ji,jj,ikv) = psal_ad(ji,jj,ikv) & |
---|
621 | & + zsjad(ji,jj) * (1 - zmaxv2) |
---|
622 | psal_ad(ji,jj,ikv-1) = psal_ad(ji,jj,ikv-1) & |
---|
623 | & + zsjad(ji,jj) * zmaxv2 |
---|
624 | zsjad(ji,jj) = 0.0_wp |
---|
625 | ENDIF |
---|
626 | |
---|
627 | ! i- direction |
---|
628 | |
---|
629 | IF( ze3wu >= 0. ) THEN ! case 1 |
---|
630 | ! gradient of T and S |
---|
631 | ztiad(ji,jj) = ztiad(ji,jj) & |
---|
632 | & + pgtu_ad(ji,jj) * umask(ji,jj,1) |
---|
633 | ptem_ad(ji,jj,iku) = ptem_ad(ji,jj,iku) & |
---|
634 | & - pgtu_ad(ji,jj) * umask(ji,jj,1) |
---|
635 | pgtu_ad(ji,jj) = 0.0_wp |
---|
636 | |
---|
637 | zsiad(ji,jj) = zsiad(ji,jj) & |
---|
638 | & + pgsu_ad(ji,jj) * umask(ji,jj,1) |
---|
639 | psal_ad(ji,jj,iku) = psal_ad(ji,jj,iku) & |
---|
640 | & - pgsu_ad(ji,jj) * umask(ji,jj,1) |
---|
641 | pgsu_ad(ji,jj) = 0.0_wp |
---|
642 | ! interpolated values of T and S |
---|
643 | ptem_ad(ji+1,jj,iku) = ptem_ad(ji+1,jj,iku) & |
---|
644 | & + ztiad(ji,jj) * (1 - zmaxu1) |
---|
645 | ptem_ad(ji+1,jj,iku-1) = ptem_ad(ji+1,jj,iku-1) & |
---|
646 | & + ztiad(ji,jj) * zmaxu1 |
---|
647 | ztiad(ji,jj) = 0.0_wp |
---|
648 | |
---|
649 | psal_ad(ji+1,jj,iku) = psal_ad(ji+1,jj,iku) & |
---|
650 | & + zsiad(ji,jj) * (1 - zmaxu1) |
---|
651 | psal_ad(ji+1,jj,iku-1) = psal_ad(ji+1,jj,iku-1) & |
---|
652 | & + zsiad(ji,jj) * zmaxu1 |
---|
653 | zsiad(ji,jj) = 0.0_wp |
---|
654 | |
---|
655 | ELSE ! case 2 |
---|
656 | ! gradient of T and S |
---|
657 | ptem_ad(ji+1,jj,iku) = ptem_ad(ji+1,jj,iku) & |
---|
658 | & + pgtu_ad(ji,jj) * umask(ji,jj,1) |
---|
659 | ztiad (ji,jj) = ztiad (ji,jj) & |
---|
660 | & - pgtu_ad(ji,jj) * umask(ji,jj,1) |
---|
661 | pgtu_ad(ji,jj) = 0.0_wp |
---|
662 | |
---|
663 | psal_ad(ji+1,jj,iku) = psal_ad(ji+1,jj,iku) & |
---|
664 | & + pgsu_ad(ji,jj) * umask(ji,jj,1) |
---|
665 | zsiad (ji,jj) = zsiad (ji,jj) & |
---|
666 | & - pgsu_ad(ji,jj) * umask(ji,jj,1) |
---|
667 | pgsu_ad(ji,jj) = 0.0_wp |
---|
668 | ! interpolated values of T and S |
---|
669 | ptem_ad(ji,jj,iku) = ptem_ad(ji,jj,iku) & |
---|
670 | & + ztiad(ji,jj) * (1 - zmaxu2) |
---|
671 | ptem_ad(ji,jj,iku-1) = ptem_ad(ji,jj,iku-1) & |
---|
672 | & + ztiad(ji,jj) * zmaxu2 |
---|
673 | ztiad(ji,jj) = 0.0_wp |
---|
674 | |
---|
675 | psal_ad(ji,jj,iku) = psal_ad(ji,jj,iku) & |
---|
676 | & + zsiad(ji,jj) * (1 - zmaxu2) |
---|
677 | psal_ad(ji,jj,iku-1) = psal_ad(ji,jj,iku-1) & |
---|
678 | & + zsiad(ji,jj) * zmaxu2 |
---|
679 | zsiad(ji,jj) = 0.0_wp |
---|
680 | ENDIF |
---|
681 | |
---|
682 | # if ! defined key_vectopt_loop |
---|
683 | END DO |
---|
684 | # endif |
---|
685 | END DO |
---|
686 | |
---|
687 | END SUBROUTINE zps_hde_adj |
---|
688 | SUBROUTINE zps_hde_adj_tst( kumadt ) |
---|
689 | !!----------------------------------------------------------------------- |
---|
690 | !! |
---|
691 | !! *** ROUTINE dyn_adv_adj_tst *** |
---|
692 | !! |
---|
693 | !! ** Purpose : Test the adjoint routine. |
---|
694 | !! |
---|
695 | !! ** Method : Verify the scalar product |
---|
696 | !! |
---|
697 | !! ( L dx )^T W dy = dx^T L^T W dy |
---|
698 | !! |
---|
699 | !! where L = tangent routine |
---|
700 | !! L^T = adjoint routine |
---|
701 | !! W = diagonal matrix of scale factors |
---|
702 | !! dx = input perturbation (random field) |
---|
703 | !! dy = L dx |
---|
704 | !! |
---|
705 | !! |
---|
706 | !! History : |
---|
707 | !! ! 08-08 (A. Vidard) |
---|
708 | !!----------------------------------------------------------------------- |
---|
709 | !! * Modules used |
---|
710 | |
---|
711 | !! * Arguments |
---|
712 | INTEGER, INTENT(IN) :: & |
---|
713 | & kumadt ! Output unit |
---|
714 | |
---|
715 | INTEGER :: & |
---|
716 | & ji, & ! dummy loop indices |
---|
717 | & jj, & |
---|
718 | & jk, & |
---|
719 | & kt |
---|
720 | INTEGER, DIMENSION(jpi,jpj) :: & |
---|
721 | & iseed_2d ! 2D seed for the random number generator |
---|
722 | |
---|
723 | !! * Local declarations |
---|
724 | REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
725 | & ztem, & ! Direct field : temperature |
---|
726 | & zsal, & ! Direct field : salinity |
---|
727 | & ztem_tlin, & ! Tangent input: temperature |
---|
728 | & zsal_tlin, & ! Tangent input: salinity |
---|
729 | & zrd_tlin, & ! Tangent input: |
---|
730 | & ztem_adout, & ! Adjoint output: temperature |
---|
731 | & zsal_adout, & ! Adjoint output: salinity |
---|
732 | & zrd_adout, & ! Adjoint output: |
---|
733 | & zar, & ! 3D random field for rd |
---|
734 | & zat, & ! 3D random field for t |
---|
735 | & zas ! 3D random field for s |
---|
736 | |
---|
737 | REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
738 | & zgtu_tlout, & ! Tangent output: horizontal gradient |
---|
739 | & zgtv_tlout, & ! Tangent output: horizontal gradient |
---|
740 | & zgsu_tlout, & ! Tangent output: horizontal gradient |
---|
741 | & zgsv_tlout, & ! Tangent output: horizontal gradient |
---|
742 | & zgru_tlout, & ! Tangent output: horizontal gradient |
---|
743 | & zgrv_tlout, & ! Tangent output: horizontal gradient |
---|
744 | & zgtu_adin, & ! Adjoint input : horizontal gradient |
---|
745 | & zgtv_adin, & ! Adjoint input : horizontal gradient |
---|
746 | & zgsu_adin, & ! Adjoint input : horizontal gradient |
---|
747 | & zgsv_adin, & ! Adjoint input : horizontal gradient |
---|
748 | & zgru_adin, & ! Adjoint input : horizontal gradient |
---|
749 | & zgrv_adin ! Adjoint input : horizontal gradient |
---|
750 | |
---|
751 | REAL(KIND=wp) :: & |
---|
752 | ! random field standard deviation for: |
---|
753 | & zsp1, & ! scalar product involving the tangent routine |
---|
754 | & zsp1_1, & ! scalar product components |
---|
755 | & zsp1_2, & |
---|
756 | & zsp1_3, & ! scalar product components |
---|
757 | & zsp1_4, & |
---|
758 | & zsp1_5, & ! scalar product components |
---|
759 | & zsp1_6, & |
---|
760 | & zsp2, & ! scalar product involving the adjoint routine |
---|
761 | & zsp2_1, & ! scalar product components |
---|
762 | & zsp2_2, & |
---|
763 | & zsp2_3 |
---|
764 | CHARACTER (LEN=14) :: & |
---|
765 | & cl_name |
---|
766 | |
---|
767 | kt = nit000 |
---|
768 | ! Allocate memory |
---|
769 | ALLOCATE( & |
---|
770 | & ztem(jpi,jpj,jpk), & |
---|
771 | & zsal(jpi,jpj,jpk), & |
---|
772 | & ztem_tlin(jpi,jpj,jpk), & |
---|
773 | & zsal_tlin(jpi,jpj,jpk), & |
---|
774 | & zrd_tlin(jpi,jpj,jpk), & |
---|
775 | & ztem_adout(jpi,jpj,jpk), & |
---|
776 | & zsal_adout(jpi,jpj,jpk), & |
---|
777 | & zrd_adout(jpi,jpj,jpk), & |
---|
778 | & zar(jpi,jpj,jpk), & |
---|
779 | & zat(jpi,jpj,jpk), & |
---|
780 | & zas(jpi,jpj,jpk), & |
---|
781 | & zgtu_tlout(jpi,jpj), & |
---|
782 | & zgtv_tlout(jpi,jpj), & |
---|
783 | & zgsu_tlout(jpi,jpj), & |
---|
784 | & zgsv_tlout(jpi,jpj), & |
---|
785 | & zgru_tlout(jpi,jpj), & |
---|
786 | & zgrv_tlout(jpi,jpj), & |
---|
787 | & zgtu_adin(jpi,jpj), & |
---|
788 | & zgtv_adin(jpi,jpj), & |
---|
789 | & zgsu_adin(jpi,jpj), & |
---|
790 | & zgsv_adin(jpi,jpj), & |
---|
791 | & zgru_adin(jpi,jpj), & |
---|
792 | & zgrv_adin(jpi,jpj) & |
---|
793 | & ) |
---|
794 | ! Initialize random field standard deviationsthe reference state |
---|
795 | ztem = tn |
---|
796 | zsal = sn |
---|
797 | |
---|
798 | !============================================================= |
---|
799 | ! 1) dx = ( T ) and dy = ( T ) |
---|
800 | !============================================================= |
---|
801 | |
---|
802 | !-------------------------------------------------------------------- |
---|
803 | ! Reset the tangent and adjoint variables |
---|
804 | !-------------------------------------------------------------------- |
---|
805 | ztem_tlin(:,:,:) = 0.0_wp |
---|
806 | zsal_tlin(:,:,:) = 0.0_wp |
---|
807 | zrd_tlin(:,:,:) = 0.0_wp |
---|
808 | ztem_adout(:,:,:) = 0.0_wp |
---|
809 | zsal_adout(:,:,:) = 0.0_wp |
---|
810 | zrd_adout(:,:,:) = 0.0_wp |
---|
811 | zgtu_tlout(:,:) = 0.0_wp |
---|
812 | zgtv_tlout(:,:) = 0.0_wp |
---|
813 | zgsu_tlout(:,:) = 0.0_wp |
---|
814 | zgsv_tlout(:,:) = 0.0_wp |
---|
815 | zgru_tlout(:,:) = 0.0_wp |
---|
816 | zgrv_tlout(:,:) = 0.0_wp |
---|
817 | zgtu_adin(:,:) = 0.0_wp |
---|
818 | zgtv_adin(:,:) = 0.0_wp |
---|
819 | zgsu_adin(:,:) = 0.0_wp |
---|
820 | zgsv_adin(:,:) = 0.0_wp |
---|
821 | zgru_adin(:,:) = 0.0_wp |
---|
822 | zgrv_adin(:,:) = 0.0_wp |
---|
823 | |
---|
824 | !-------------------------------------------------------------------- |
---|
825 | ! Initialize the tangent input with random noise: dx |
---|
826 | !-------------------------------------------------------------------- |
---|
827 | DO jj = 1, jpj |
---|
828 | DO ji = 1, jpi |
---|
829 | iseed_2d(ji,jj) = - ( 456953 + & |
---|
830 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
831 | END DO |
---|
832 | END DO |
---|
833 | CALL grid_random( iseed_2d, zat, 'T', 0.0_wp, stdt ) |
---|
834 | DO jj = 1, jpj |
---|
835 | DO ji = 1, jpi |
---|
836 | iseed_2d(ji,jj) = - ( 526791 + & |
---|
837 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
838 | END DO |
---|
839 | END DO |
---|
840 | CALL grid_random( iseed_2d, zas, 'T', 0.0_wp, stds ) |
---|
841 | |
---|
842 | DO jj = 1, jpj |
---|
843 | DO ji = 1, jpi |
---|
844 | iseed_2d(ji,jj) = - ( 395703 + & |
---|
845 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
846 | END DO |
---|
847 | END DO |
---|
848 | CALL grid_random( iseed_2d, zar, 'T', 0.0_wp, stdr ) |
---|
849 | |
---|
850 | DO jk = 1, jpk |
---|
851 | DO jj = nldj, nlej |
---|
852 | DO ji = nldi, nlei |
---|
853 | ztem_tlin(ji,jj,jk) = zat(ji,jj,jk) |
---|
854 | zsal_tlin(ji,jj,jk) = zas(ji,jj,jk) |
---|
855 | zrd_tlin( ji,jj,jk) = zar(ji,jj,jk) |
---|
856 | END DO |
---|
857 | END DO |
---|
858 | END DO |
---|
859 | |
---|
860 | CALL zps_hde_tan ( kt, ztem , zsal , & |
---|
861 | & ztem_tlin , zsal_tlin , zrd_tlin , & |
---|
862 | & zgtu_tlout, zgsu_tlout, zgru_tlout, & |
---|
863 | & zgtv_tlout, zgsv_tlout, zgrv_tlout ) |
---|
864 | DO jj = nldj, nlej |
---|
865 | DO ji = nldi, nlei |
---|
866 | jk = mbatu(ji,jj) |
---|
867 | zgtu_adin(ji,jj) = zgtu_tlout(ji,jj) & |
---|
868 | & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & |
---|
869 | & * umask(ji,jj,jk) |
---|
870 | zgsu_adin(ji,jj) = zgsu_tlout(ji,jj) & |
---|
871 | & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & |
---|
872 | & * umask(ji,jj,jk) |
---|
873 | zgru_adin(ji,jj) = zgru_tlout(ji,jj) & |
---|
874 | & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & |
---|
875 | & * umask(ji,jj,jk) |
---|
876 | jk = mbatv(ji,jj) |
---|
877 | zgtv_adin(ji,jj) = zgtv_tlout(ji,jj) & |
---|
878 | & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & |
---|
879 | & * vmask(ji,jj,jk) |
---|
880 | zgsv_adin(ji,jj) = zgsv_tlout(ji,jj) & |
---|
881 | & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & |
---|
882 | & * vmask(ji,jj,jk) |
---|
883 | zgrv_adin(ji,jj) = zgrv_tlout(ji,jj) & |
---|
884 | & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & |
---|
885 | & * vmask(ji,jj,jk) |
---|
886 | END DO |
---|
887 | END DO |
---|
888 | !-------------------------------------------------------------------- |
---|
889 | ! Compute the scalar product: ( L dx )^T W dy |
---|
890 | !-------------------------------------------------------------------- |
---|
891 | |
---|
892 | zsp1_1 = DOT_PRODUCT( zgtu_adin, zgtu_tlout ) |
---|
893 | zsp1_2 = DOT_PRODUCT( zgsu_adin, zgsu_tlout ) |
---|
894 | zsp1_3 = DOT_PRODUCT( zgru_adin, zgru_tlout ) |
---|
895 | zsp1_4 = DOT_PRODUCT( zgtv_adin, zgtv_tlout ) |
---|
896 | zsp1_5 = DOT_PRODUCT( zgsv_adin, zgsv_tlout ) |
---|
897 | zsp1_6 = DOT_PRODUCT( zgrv_adin, zgrv_tlout ) |
---|
898 | zsp1 = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5 + zsp1_6 |
---|
899 | |
---|
900 | |
---|
901 | !-------------------------------------------------------------------- |
---|
902 | ! Call the adjoint routine: dx^* = L^T dy^* |
---|
903 | !-------------------------------------------------------------------- |
---|
904 | CALL zps_hde_adj ( kt, ztem , zsal , & |
---|
905 | & ztem_adout, zsal_adout, zrd_adout , & |
---|
906 | & zgtu_adin , zgsu_adin , zgru_adin , & |
---|
907 | & zgtv_adin , zgsv_adin , zgrv_adin ) |
---|
908 | zsp2_1 = DOT_PRODUCT( ztem_tlin, ztem_adout ) |
---|
909 | zsp2_2 = DOT_PRODUCT( zsal_tlin, zsal_adout ) |
---|
910 | zsp2_3 = DOT_PRODUCT( zrd_tlin , zrd_adout ) |
---|
911 | zsp2 = zsp2_1 + zsp2_2 + zsp2_3 |
---|
912 | |
---|
913 | ! Compare the scalar products |
---|
914 | |
---|
915 | cl_name = 'zps_hde_adj ' |
---|
916 | CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) |
---|
917 | |
---|
918 | ! Deallocate memory |
---|
919 | DEALLOCATE( & |
---|
920 | & ztem, & |
---|
921 | & zsal, & |
---|
922 | & ztem_tlin, & |
---|
923 | & zsal_tlin, & |
---|
924 | & zrd_tlin, & |
---|
925 | & ztem_adout, & |
---|
926 | & zsal_adout, & |
---|
927 | & zrd_adout, & |
---|
928 | & zar, & |
---|
929 | & zat, & |
---|
930 | & zas, & |
---|
931 | & zgtu_tlout, & |
---|
932 | & zgtv_tlout, & |
---|
933 | & zgsu_tlout, & |
---|
934 | & zgsv_tlout, & |
---|
935 | & zgru_tlout, & |
---|
936 | & zgrv_tlout, & |
---|
937 | & zgtu_adin, & |
---|
938 | & zgtv_adin, & |
---|
939 | & zgsu_adin, & |
---|
940 | & zgsv_adin, & |
---|
941 | & zgru_adin, & |
---|
942 | & zgrv_adin & |
---|
943 | & ) |
---|
944 | |
---|
945 | |
---|
946 | |
---|
947 | END SUBROUTINE zps_hde_adj_tst |
---|
948 | |
---|
949 | SUBROUTINE zps_hde_tlm_tst( kumadt ) |
---|
950 | !!----------------------------------------------------------------------- |
---|
951 | !! |
---|
952 | !! *** ROUTINE dyn_adv_tlm_tst *** |
---|
953 | !! |
---|
954 | !! ** Purpose : Test the adjoint routine. |
---|
955 | !! |
---|
956 | !! ** Method : Verify the tangent with Taylor expansion |
---|
957 | !! |
---|
958 | !! M(x+hdx) = M(x) + L(hdx) + O(h^2) |
---|
959 | !! |
---|
960 | !! where L = tangent routine |
---|
961 | !! M = direct routine |
---|
962 | !! dx = input perturbation (random field) |
---|
963 | !! h = ration on perturbation |
---|
964 | !! |
---|
965 | !! In the tangent test we verify that: |
---|
966 | !! M(x+h*dx) - M(x) |
---|
967 | !! g(h) = ------------------ ---> 1 as h ---> 0 |
---|
968 | !! L(h*dx) |
---|
969 | !! and |
---|
970 | !! g(h) - 1 |
---|
971 | !! f(h) = ---------- ---> k (costant) as h ---> 0 |
---|
972 | !! p |
---|
973 | !! |
---|
974 | !! History : |
---|
975 | !! ! 09-08 (A. Vigilant) |
---|
976 | !!----------------------------------------------------------------------- |
---|
977 | !! * Modules used |
---|
978 | USE zpshde |
---|
979 | USE oce_tam, ONLY: & |
---|
980 | & gtu_tl, gtv_tl, & |
---|
981 | & gsu_tl, gsv_tl, & |
---|
982 | & gru_tl, grv_tl, & |
---|
983 | & tn_tl, sn_tl, rhd_tl |
---|
984 | USE tamtrj ! writing out state trajectory |
---|
985 | USE par_tlm, ONLY: & |
---|
986 | & cur_loop, & |
---|
987 | & h_ratio |
---|
988 | USE istate_mod |
---|
989 | USE divcur ! horizontal divergence and relative vorticity |
---|
990 | USE wzvmod ! vertical velocity |
---|
991 | USE gridrandom, ONLY: & |
---|
992 | & grid_rd_sd |
---|
993 | USE trj_tam |
---|
994 | USE oce , ONLY: & ! ocean dynamics and tracers variables |
---|
995 | & tn, sn, rhd, & |
---|
996 | & gtu, gtv, gsu, gsv, & |
---|
997 | & gru, grv |
---|
998 | USE opatam_tst_ini, ONLY: & |
---|
999 | & tlm_namrd |
---|
1000 | USE tamctl, ONLY: & ! Control parameters |
---|
1001 | & numtan, numtan_sc |
---|
1002 | !! * Arguments |
---|
1003 | INTEGER, INTENT(IN) :: & |
---|
1004 | & kumadt ! Output unit |
---|
1005 | |
---|
1006 | !! * Local declarations |
---|
1007 | INTEGER :: & |
---|
1008 | & ji, & ! dummy loop indices |
---|
1009 | & jj, & |
---|
1010 | & jk |
---|
1011 | !! * Local declarations |
---|
1012 | REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
1013 | & ztem, & ! Direct field : temperature |
---|
1014 | & zsal, & ! Direct field : salinity |
---|
1015 | & ztn_tlin, & ! Tangent input: temperature |
---|
1016 | & zsn_tlin, & ! Tangent input: salinity |
---|
1017 | & zrd_tlin, & |
---|
1018 | & z3r ! 3D field |
---|
1019 | |
---|
1020 | REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
1021 | & zgtu_out, zgtu_wop, & ! Tangent output: horizontal gradient |
---|
1022 | & zgtv_out, zgtv_wop, & ! Tangent output: horizontal gradient |
---|
1023 | & zgsu_out, zgsu_wop, & ! Tangent output: horizontal gradient |
---|
1024 | & zgsv_out, zgsv_wop, & ! Tangent output: horizontal gradient |
---|
1025 | & zgru_out, zgru_wop, & ! Tangent output: horizontal gradient |
---|
1026 | & zgrv_out, zgrv_wop |
---|
1027 | |
---|
1028 | REAL(KIND=wp) :: & |
---|
1029 | & zsp1, & ! scalar product involving the tangent routine |
---|
1030 | & zsp1_1, zsp1_2, & |
---|
1031 | & zsp1_3, zsp1_4, & |
---|
1032 | & zsp1_5, zsp1_6, & |
---|
1033 | & zsp2, & |
---|
1034 | & zsp2_1, zsp2_2, & |
---|
1035 | & zsp2_3, zsp2_4, & |
---|
1036 | & zsp2_5, zsp2_6, & |
---|
1037 | & zsp3, & |
---|
1038 | & zsp3_1, zsp3_2, & |
---|
1039 | & zsp3_3, zsp3_4, & |
---|
1040 | & zsp3_5, zsp3_6, & |
---|
1041 | & zzsp, & |
---|
1042 | & zzsp_1, zzsp_2, & |
---|
1043 | & zzsp_3, zzsp_4, & |
---|
1044 | & zzsp_5, zzsp_6, & |
---|
1045 | & gamma, & |
---|
1046 | & zgsp1, zgsp2, & |
---|
1047 | & zgsp3, zgsp4, & |
---|
1048 | & zgsp5, zgsp6, & |
---|
1049 | & zgsp7 |
---|
1050 | CHARACTER(LEN=14) :: cl_name |
---|
1051 | CHARACTER (LEN=128) :: file_out, file_wop |
---|
1052 | CHARACTER (LEN=90) :: FMT |
---|
1053 | REAL(KIND=wp), DIMENSION(100):: & |
---|
1054 | & zscgtu, zscgtv, zscgsu, & |
---|
1055 | & zscgsv, zscgru, zscgrv, & |
---|
1056 | & zscerrgtu, zscerrgtv, & |
---|
1057 | & zscerrgsu, zscerrgsv, & |
---|
1058 | & zscerrgru, zscerrgrv |
---|
1059 | INTEGER, DIMENSION(100):: & |
---|
1060 | & iiposgtu, iiposgtv, iiposgsu, & |
---|
1061 | & iiposgsv, iiposgru, iiposgrv, & |
---|
1062 | & ijposgtu, ijposgtv, ijposgsu, & |
---|
1063 | & ijposgsv, ijposgru, ijposgrv |
---|
1064 | INTEGER:: & |
---|
1065 | & ii, & |
---|
1066 | & isamp=40, & |
---|
1067 | & jsamp=40, & |
---|
1068 | & numsctlm |
---|
1069 | REAL(KIND=wp), DIMENSION(jpi,jpj) :: & |
---|
1070 | & zerrgtu, zerrgtv, zerrgsu, & |
---|
1071 | & zerrgsv, zerrgru, zerrgrv |
---|
1072 | ! Allocate memory |
---|
1073 | |
---|
1074 | ALLOCATE( & |
---|
1075 | & ztem(jpi,jpj,jpk), & |
---|
1076 | & zsal(jpi,jpj,jpk), & |
---|
1077 | & ztn_tlin(jpi,jpj,jpk), & |
---|
1078 | & zsn_tlin(jpi,jpj,jpk), & |
---|
1079 | & zrd_tlin(jpi,jpj,jpk), & |
---|
1080 | & z3r (jpi,jpj,jpk), & |
---|
1081 | & zgtu_out(jpi,jpj), & |
---|
1082 | & zgtv_out(jpi,jpj), & |
---|
1083 | & zgsu_out(jpi,jpj), & |
---|
1084 | & zgsv_out(jpi,jpj), & |
---|
1085 | & zgru_out(jpi,jpj), & |
---|
1086 | & zgrv_out(jpi,jpj), & |
---|
1087 | & zgtu_wop(jpi,jpj), & |
---|
1088 | & zgtv_wop(jpi,jpj), & |
---|
1089 | & zgsu_wop(jpi,jpj), & |
---|
1090 | & zgsv_wop(jpi,jpj), & |
---|
1091 | & zgru_wop(jpi,jpj), & |
---|
1092 | & zgrv_wop(jpi,jpj) & |
---|
1093 | & ) |
---|
1094 | !-------------------------------------------------------------------- |
---|
1095 | ! Reset the tangent and adjoint variables |
---|
1096 | !-------------------------------------------------------------------- |
---|
1097 | ztn_tlin(:,:,:) = 0.0_wp |
---|
1098 | zsn_tlin(:,:,:) = 0.0_wp |
---|
1099 | zrd_tlin(:,:,:) = 0.0_wp |
---|
1100 | zgtu_out(:,:) = 0.0_wp |
---|
1101 | zgtv_out(:,:) = 0.0_wp |
---|
1102 | zgsu_out(:,:) = 0.0_wp |
---|
1103 | zgsv_out(:,:) = 0.0_wp |
---|
1104 | zgru_out(:,:) = 0.0_wp |
---|
1105 | zgrv_out(:,:) = 0.0_wp |
---|
1106 | gtu_tl(:,:) = 0.0_wp |
---|
1107 | gtv_tl(:,:) = 0.0_wp |
---|
1108 | gsu_tl(:,:) = 0.0_wp |
---|
1109 | gsv_tl(:,:) = 0.0_wp |
---|
1110 | gru_tl(:,:) = 0.0_wp |
---|
1111 | grv_tl(:,:) = 0.0_wp |
---|
1112 | |
---|
1113 | zscgtu(:) = 0.0_wp |
---|
1114 | zscgtv(:) = 0.0_wp |
---|
1115 | zscgsu(:) = 0.0_wp |
---|
1116 | zscgsv(:) = 0.0_wp |
---|
1117 | zscgru(:) = 0.0_wp |
---|
1118 | zscgrv(:) = 0.0_wp |
---|
1119 | zscerrgtu(:) = 0.0_wp |
---|
1120 | zscerrgtv(:) = 0.0_wp |
---|
1121 | zscerrgsu(:) = 0.0_wp |
---|
1122 | zscerrgsv(:) = 0.0_wp |
---|
1123 | zscerrgru(:) = 0.0_wp |
---|
1124 | zscerrgrv(:) = 0.0_wp |
---|
1125 | zerrgtu(:,:) = 0.0_wp |
---|
1126 | zerrgtv(:,:) = 0.0_wp |
---|
1127 | zerrgsu(:,:) = 0.0_wp |
---|
1128 | zerrgsv(:,:) = 0.0_wp |
---|
1129 | zerrgru(:,:) = 0.0_wp |
---|
1130 | zerrgrv(:,:) = 0.0_wp |
---|
1131 | !-------------------------------------------------------------------- |
---|
1132 | ! Output filename Xn=F(X0) |
---|
1133 | !-------------------------------------------------------------------- |
---|
1134 | file_wop='trj_wop_zps' |
---|
1135 | CALL tlm_namrd |
---|
1136 | gamma = h_ratio |
---|
1137 | !-------------------------------------------------------------------- |
---|
1138 | ! Initialize the tangent input with random noise: dx |
---|
1139 | !-------------------------------------------------------------------- |
---|
1140 | IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN |
---|
1141 | CALL grid_rd_sd( 596035, z3r, 'T', 0.0_wp, stdt) |
---|
1142 | DO jk = 1, jpk |
---|
1143 | DO jj = nldj, nlej |
---|
1144 | DO ji = nldi, nlei |
---|
1145 | ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1146 | END DO |
---|
1147 | END DO |
---|
1148 | END DO |
---|
1149 | CALL grid_rd_sd( 523432, z3r, 'T', 0.0_wp, stds) |
---|
1150 | DO jk = 1, jpk |
---|
1151 | DO jj = nldj, nlej |
---|
1152 | DO ji = nldi, nlei |
---|
1153 | zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1154 | END DO |
---|
1155 | END DO |
---|
1156 | END DO |
---|
1157 | CALL grid_rd_sd( 432545, z3r, 'T', 0.0_wp, stdr) |
---|
1158 | DO jk = 1, jpk |
---|
1159 | DO jj = nldj, nlej |
---|
1160 | DO ji = nldi, nlei |
---|
1161 | zrd_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1162 | END DO |
---|
1163 | END DO |
---|
1164 | END DO |
---|
1165 | ENDIF |
---|
1166 | !-------------------------------------------------------------------- |
---|
1167 | ! Complete Init for Direct |
---|
1168 | !------------------------------------------------------------------- |
---|
1169 | CALL istate_p |
---|
1170 | |
---|
1171 | ! *** initialize the reference trajectory |
---|
1172 | ! ------------ |
---|
1173 | CALL trj_rea( nit000-1, 1 ) |
---|
1174 | CALL trj_rea( nit000, 1 ) |
---|
1175 | |
---|
1176 | IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN |
---|
1177 | ztn_tlin(:,:,:) = gamma * ztn_tlin(:,:,:) |
---|
1178 | tn(:,:,:) = tn(:,:,:) + ztn_tlin(:,:,:) |
---|
1179 | |
---|
1180 | zsn_tlin(:,:,:) = gamma * zsn_tlin(:,:,:) |
---|
1181 | sn(:,:,:) = sn(:,:,:) + zsn_tlin(:,:,:) |
---|
1182 | |
---|
1183 | zrd_tlin(:,:,:) = gamma * zrd_tlin(:,:,:) |
---|
1184 | rhd(:,:,:) = rhd(:,:,:) + zrd_tlin(:,:,:) |
---|
1185 | ENDIF |
---|
1186 | !-------------------------------------------------------------------- |
---|
1187 | ! Compute the direct model F(X0,t=n) = Xn |
---|
1188 | !-------------------------------------------------------------------- |
---|
1189 | CALL zps_hde(nit000, tn, sn, rhd, gtu, gsu, gru, gtv, gsv, grv) |
---|
1190 | |
---|
1191 | IF ( cur_loop .EQ. 0) CALL trj_wri_spl(file_wop) |
---|
1192 | |
---|
1193 | !-------------------------------------------------------------------- |
---|
1194 | ! Compute the Tangent |
---|
1195 | !-------------------------------------------------------------------- |
---|
1196 | IF ( cur_loop .NE. 0) THEN |
---|
1197 | !-------------------------------------------------------------------- |
---|
1198 | ! Storing data |
---|
1199 | !-------------------------------------------------------------------- |
---|
1200 | zgtu_out (:,:) = gtu (:,:) |
---|
1201 | zgtv_out (:,:) = gtv (:,:) |
---|
1202 | zgsu_out (:,:) = gsu (:,:) |
---|
1203 | zgsv_out (:,:) = gsv (:,:) |
---|
1204 | zgru_out (:,:) = gru (:,:) |
---|
1205 | zgrv_out (:,:) = grv (:,:) |
---|
1206 | |
---|
1207 | !-------------------------------------------------------------------- |
---|
1208 | ! Initialize the tangent variables: |
---|
1209 | !-------------------------------------------------------------------- |
---|
1210 | CALL trj_rea( nit000-1, 1 ) |
---|
1211 | CALL trj_rea( nit000, 1 ) |
---|
1212 | tn_tl (:,:,:) = ztn_tlin (:,:,:) |
---|
1213 | sn_tl (:,:,:) = zsn_tlin (:,:,:) |
---|
1214 | rhd_tl (:,:,:) = zrd_tlin (:,:,:) |
---|
1215 | CALL zps_hde_tan(nit000, tn , sn , & |
---|
1216 | & ztn_tlin , zsn_tlin , zrd_tlin, & |
---|
1217 | & gtu_tl, gsu_tl, gru_tl, & |
---|
1218 | & gtv_tl, gsv_tl, grv_tl ) |
---|
1219 | |
---|
1220 | !-------------------------------------------------------------------- |
---|
1221 | ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) |
---|
1222 | !-------------------------------------------------------------------- |
---|
1223 | |
---|
1224 | zsp2_1 = DOT_PRODUCT( gtu_tl, gtu_tl ) |
---|
1225 | zsp2_2 = DOT_PRODUCT( gtv_tl, gtv_tl ) |
---|
1226 | zsp2_3 = DOT_PRODUCT( gsu_tl, gsu_tl ) |
---|
1227 | zsp2_4 = DOT_PRODUCT( gsv_tl, gsv_tl ) |
---|
1228 | zsp2_5 = DOT_PRODUCT( gru_tl, gru_tl ) |
---|
1229 | zsp2_6 = DOT_PRODUCT( grv_tl, grv_tl ) |
---|
1230 | |
---|
1231 | zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 + zsp2_6 |
---|
1232 | !-------------------------------------------------------------------- |
---|
1233 | ! Storing data |
---|
1234 | !-------------------------------------------------------------------- |
---|
1235 | |
---|
1236 | CALL trj_rd_spl(file_wop) |
---|
1237 | zgtu_wop (:,:) = gtu (:,:) |
---|
1238 | zgtv_wop (:,:) = gtv (:,:) |
---|
1239 | zgsu_wop (:,:) = gsu (:,:) |
---|
1240 | zgsv_wop (:,:) = gsv (:,:) |
---|
1241 | zgru_wop (:,:) = gru (:,:) |
---|
1242 | zgrv_wop (:,:) = grv (:,:) |
---|
1243 | !-------------------------------------------------------------------- |
---|
1244 | ! Compute the Linearization Error |
---|
1245 | ! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn) |
---|
1246 | ! and |
---|
1247 | ! Compute the Linearization Error |
---|
1248 | ! En = Nn -TL(gamma.dX0, t0,tn) |
---|
1249 | !-------------------------------------------------------------------- |
---|
1250 | ! Warning: Here we re-use local variables z()_out and z()_wop |
---|
1251 | ii=0 |
---|
1252 | DO jj = 1, jpj |
---|
1253 | DO ji = 1, jpi |
---|
1254 | zgtu_out (ji,jj) = zgtu_out (ji,jj) - zgtu_wop (ji,jj) |
---|
1255 | zgtu_wop (ji,jj) = zgtu_out (ji,jj) - gtu_tl (ji,jj) |
---|
1256 | IF ( gtu_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1257 | & zerrgtu(ji,jj) = zgtu_out(ji,jj)/gtu_tl(ji,jj) |
---|
1258 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1259 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1260 | ii = ii+1 |
---|
1261 | iiposgtu(ii) = ji |
---|
1262 | ijposgtu(ii) = jj |
---|
1263 | IF ( INT(umask(ji,jj,1)) .NE. 0) THEN |
---|
1264 | zscgtu (ii) = zgtu_wop(ji,jj) |
---|
1265 | zscerrgtu (ii) = ( zerrgtu(ji,jj) - 1.0_wp ) /gamma |
---|
1266 | ENDIF |
---|
1267 | ENDIF |
---|
1268 | END DO |
---|
1269 | END DO |
---|
1270 | ii=0 |
---|
1271 | DO jj = 1, jpj |
---|
1272 | DO ji = 1, jpi |
---|
1273 | zgtv_out (ji,jj) = zgtv_out (ji,jj) - zgtv_wop (ji,jj) |
---|
1274 | zgtv_wop (ji,jj) = zgtv_out (ji,jj) - gtv_tl (ji,jj) |
---|
1275 | IF ( gtv_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1276 | & zerrgtv(ji,jj) = zgtv_out(ji,jj)/gtv_tl(ji,jj) |
---|
1277 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1278 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1279 | ii = ii+1 |
---|
1280 | iiposgtv(ii) = ji |
---|
1281 | ijposgtv(ii) = jj |
---|
1282 | IF ( INT(vmask(ji,jj,1)) .NE. 0) THEN |
---|
1283 | zscgtv (ii) = zgtv_wop(ji,jj) |
---|
1284 | zscerrgtv (ii) = ( zerrgtv(ji,jj) - 1.0_wp ) / gamma |
---|
1285 | ENDIF |
---|
1286 | ENDIF |
---|
1287 | END DO |
---|
1288 | END DO |
---|
1289 | ii=0 |
---|
1290 | DO jj = 1, jpj |
---|
1291 | DO ji = 1, jpi |
---|
1292 | zgsu_out (ji,jj) = zgsu_out (ji,jj) - zgsu_wop (ji,jj) |
---|
1293 | zgsu_wop (ji,jj) = zgsu_out (ji,jj) - gsu_tl (ji,jj) |
---|
1294 | IF ( gsu_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1295 | & zerrgsu(ji,jj) = zgsu_out(ji,jj)/gsu_tl(ji,jj) |
---|
1296 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1297 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1298 | ii = ii+1 |
---|
1299 | iiposgsu(ii) = ji |
---|
1300 | ijposgsu(ii) = jj |
---|
1301 | IF ( INT(umask(ji,jj,1)) .NE. 0) THEN |
---|
1302 | zscgsu (ii) = zgsu_wop(ji,jj) |
---|
1303 | zscerrgsu (ii) = ( zerrgsu(ji,jj) - 1.0_wp ) /gamma |
---|
1304 | ENDIF |
---|
1305 | ENDIF |
---|
1306 | END DO |
---|
1307 | END DO |
---|
1308 | ii=0 |
---|
1309 | DO jj = 1, jpj |
---|
1310 | DO ji = 1, jpi |
---|
1311 | zgsv_out (ji,jj) = zgsv_out (ji,jj) - zgsv_wop (ji,jj) |
---|
1312 | zgsv_wop (ji,jj) = zgsv_out (ji,jj) - gsv_tl (ji,jj) |
---|
1313 | IF ( gtv_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1314 | & zerrgsv(ji,jj) = zgsv_out(ji,jj)/gsv_tl(ji,jj) |
---|
1315 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1316 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1317 | ii = ii+1 |
---|
1318 | iiposgsv(ii) = ji |
---|
1319 | ijposgsv(ii) = jj |
---|
1320 | IF ( INT(vmask(ji,jj,1)) .NE. 0) THEN |
---|
1321 | zscgsv (ii) = zgsv_wop(ji,jj) |
---|
1322 | zscerrgsv (ii) = ( zerrgsv(ji,jj) - 1.0_wp ) /gamma |
---|
1323 | ENDIF |
---|
1324 | ENDIF |
---|
1325 | END DO |
---|
1326 | END DO |
---|
1327 | ii=0 |
---|
1328 | DO jj = 1, jpj |
---|
1329 | DO ji = 1, jpi |
---|
1330 | zgru_out (ji,jj) = zgru_out (ji,jj) - zgru_wop (ji,jj) |
---|
1331 | zgru_wop (ji,jj) = zgru_out (ji,jj) - gru_tl (ji,jj) |
---|
1332 | IF ( gru_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1333 | & zerrgru(ji,jj) = zgru_out(ji,jj)/gru_tl(ji,jj) |
---|
1334 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1335 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1336 | ii = ii+1 |
---|
1337 | iiposgru(ii) = ji |
---|
1338 | ijposgru(ii) = jj |
---|
1339 | IF ( INT(umask(ji,jj,1)) .NE. 0) THEN |
---|
1340 | zscgru (ii) = zgru_wop(ji,jj) |
---|
1341 | zscerrgru (ii) = ( zerrgru(ji,jj) - 1.0_wp ) /gamma |
---|
1342 | ENDIF |
---|
1343 | ENDIF |
---|
1344 | END DO |
---|
1345 | END DO |
---|
1346 | ii=0 |
---|
1347 | DO jj = 1, jpj |
---|
1348 | DO ji = 1, jpi |
---|
1349 | zgrv_out (ji,jj) = zgrv_out (ji,jj) - zgrv_wop (ji,jj) |
---|
1350 | zgrv_wop (ji,jj) = zgrv_out (ji,jj) - grv_tl (ji,jj) |
---|
1351 | IF ( grv_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1352 | & zerrgrv(ji,jj) = zgrv_out(ji,jj)/grv_tl(ji,jj) |
---|
1353 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1354 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1355 | ii = ii+1 |
---|
1356 | iiposgrv(ii) = ji |
---|
1357 | ijposgrv(ii) = jj |
---|
1358 | IF ( INT(vmask(ji,jj,1)) .NE. 0) THEN |
---|
1359 | zscgrv (ii) = zgrv_wop(ji,jj) |
---|
1360 | zscerrgrv (ii) = ( zerrgrv(ji,jj) - 1.0_wp ) / gamma |
---|
1361 | ENDIF |
---|
1362 | ENDIF |
---|
1363 | END DO |
---|
1364 | END DO |
---|
1365 | zsp1_1 = DOT_PRODUCT( zgtu_out, zgtu_out ) |
---|
1366 | zsp1_2 = DOT_PRODUCT( zgtv_out, zgtv_out ) |
---|
1367 | zsp1_3 = DOT_PRODUCT( zgsu_out, zgsu_out ) |
---|
1368 | zsp1_4 = DOT_PRODUCT( zgsv_out, zgsv_out ) |
---|
1369 | zsp1_5 = DOT_PRODUCT( zgru_out, zgru_out ) |
---|
1370 | zsp1_6 = DOT_PRODUCT( zgrv_out, zgrv_out ) |
---|
1371 | zsp1 = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5 + zsp1_6 |
---|
1372 | zsp3_1 = DOT_PRODUCT( zgtu_wop, zgtu_wop ) |
---|
1373 | zsp3_2 = DOT_PRODUCT( zgtv_wop, zgtv_wop ) |
---|
1374 | zsp3_3 = DOT_PRODUCT( zgsu_wop, zgsu_wop ) |
---|
1375 | zsp3_4 = DOT_PRODUCT( zgsv_wop, zgsv_wop ) |
---|
1376 | zsp3_5 = DOT_PRODUCT( zgru_wop, zgru_wop ) |
---|
1377 | zsp3_6 = DOT_PRODUCT( zgrv_wop, zgrv_wop ) |
---|
1378 | zsp3 = zsp3_1 + zsp3_2 + zsp3_3 + zsp3_4 + zsp3_5 + zsp3_6 |
---|
1379 | !-------------------------------------------------------------------- |
---|
1380 | ! Print the linearization error En - norme 2 |
---|
1381 | !-------------------------------------------------------------------- |
---|
1382 | ! 14 char:'12345678901234' |
---|
1383 | cl_name = 'zps_tam:En ' |
---|
1384 | zzsp = SQRT(zsp3) |
---|
1385 | zzsp_1 = SQRT(zsp3_1) |
---|
1386 | zzsp_2 = SQRT(zsp3_2) |
---|
1387 | zzsp_3 = SQRT(zsp3_3) |
---|
1388 | zzsp_4 = SQRT(zsp3_4) |
---|
1389 | zzsp_5 = SQRT(zsp3_5) |
---|
1390 | zzsp_6 = SQRT(zsp3_6) |
---|
1391 | zgsp5 = zzsp |
---|
1392 | CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) |
---|
1393 | !-------------------------------------------------------------------- |
---|
1394 | ! Compute TLM norm2 |
---|
1395 | !-------------------------------------------------------------------- |
---|
1396 | zzsp = SQRT(zsp2) |
---|
1397 | zzsp_1 = SQRT(zsp2_1) |
---|
1398 | zzsp_2 = SQRT(zsp2_2) |
---|
1399 | zzsp_3 = SQRT(zsp2_3) |
---|
1400 | zzsp_4 = SQRT(zsp2_4) |
---|
1401 | zzsp_5 = SQRT(zsp2_5) |
---|
1402 | zzsp_6 = SQRT(zsp2_6) |
---|
1403 | zgsp4 = zzsp |
---|
1404 | cl_name = 'zps_tam:Ln2' |
---|
1405 | CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) |
---|
1406 | !-------------------------------------------------------------------- |
---|
1407 | ! Print the linearization error Nn - norme 2 |
---|
1408 | !-------------------------------------------------------------------- |
---|
1409 | zzsp = SQRT(zsp1) |
---|
1410 | zzsp_1 = SQRT(zsp1_1) |
---|
1411 | zzsp_2 = SQRT(zsp1_2) |
---|
1412 | zzsp_3 = SQRT(zsp1_3) |
---|
1413 | zzsp_4 = SQRT(zsp1_4) |
---|
1414 | zzsp_5 = SQRT(zsp1_5) |
---|
1415 | zzsp_6 = SQRT(zsp1_6) |
---|
1416 | cl_name = 'zpshde:Mhdx-Mx' |
---|
1417 | CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) |
---|
1418 | |
---|
1419 | zgsp3 = SQRT( zsp3/zsp2 ) |
---|
1420 | zgsp7 = zgsp3/gamma |
---|
1421 | zgsp1 = zzsp |
---|
1422 | zgsp2 = zgsp1 / zgsp4 |
---|
1423 | zgsp6 = (zgsp2 - 1.0_wp)/gamma |
---|
1424 | |
---|
1425 | FMT = "(A8,2X,I4.4,2X,E6.1,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13)" |
---|
1426 | WRITE(numtan,FMT) 'zpshde ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7 |
---|
1427 | !-------------------------------------------------------------------- |
---|
1428 | ! Unitary calculus |
---|
1429 | !-------------------------------------------------------------------- |
---|
1430 | FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)" |
---|
1431 | cl_name = 'zpshde ' |
---|
1432 | IF (lwp) THEN |
---|
1433 | DO ii=1, 100, 1 |
---|
1434 | IF ( zscgtu(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgtu ', & |
---|
1435 | & cur_loop, h_ratio, ii, iiposgtu(ii), ijposgtu(ii), zscgtu(ii) |
---|
1436 | ENDDO |
---|
1437 | DO ii=1, 100, 1 |
---|
1438 | IF ( zscgtv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgtv ', & |
---|
1439 | & cur_loop, h_ratio, ii, iiposgtv(ii), ijposgtv(ii), zscgtv(ii) |
---|
1440 | ENDDO |
---|
1441 | DO ii=1, 100, 1 |
---|
1442 | IF ( zscgsu(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgsu ', & |
---|
1443 | & cur_loop, h_ratio, ii, iiposgsu(ii), ijposgsu(ii), zscgsu(ii) |
---|
1444 | ENDDO |
---|
1445 | DO ii=1, 100, 1 |
---|
1446 | IF ( zscgsv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgsv ', & |
---|
1447 | & cur_loop, h_ratio, ii, iiposgsv(ii), ijposgsv(ii), zscgsv(ii) |
---|
1448 | ENDDO |
---|
1449 | DO ii=1, 100, 1 |
---|
1450 | IF ( zscgru(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgru ', & |
---|
1451 | & cur_loop, h_ratio, ii, iiposgru(ii), ijposgru(ii), zscgru(ii) |
---|
1452 | ENDDO |
---|
1453 | DO ii=1, 100, 1 |
---|
1454 | IF ( zscgrv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscgrv ', & |
---|
1455 | & cur_loop, h_ratio, ii, iiposgrv(ii), ijposgrv(ii), zscgrv(ii) |
---|
1456 | ENDDO |
---|
1457 | ! write separator |
---|
1458 | WRITE(numtan_sc,"(A4)") '====' |
---|
1459 | ENDIF |
---|
1460 | ENDIF |
---|
1461 | DEALLOCATE( & |
---|
1462 | & ztem, zsal, & |
---|
1463 | & ztn_tlin, zsn_tlin, & |
---|
1464 | & zrd_tlin, zgtu_out, & |
---|
1465 | & zgtv_out, zgsu_out, & |
---|
1466 | & zgsv_out, zgru_out, & |
---|
1467 | & zgrv_out, & |
---|
1468 | & z3r & |
---|
1469 | & ) |
---|
1470 | END SUBROUTINE zps_hde_tlm_tst |
---|
1471 | !!====================================================================== |
---|
1472 | #endif |
---|
1473 | END MODULE zpshde_tam |
---|