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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

obs_oper.F90 in branches/devukmo2010/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/devukmo2010/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 2128

Last change on this file since 2128 was 2128, checked in by rfurner, 14 years ago
merged branches OBS, ASM, Rivers, BDY & mixed_dynldf ready for vn3.3 merge
File size: 47.8 KB

Line
1	MODULE obs_oper
2	!!======================================================================
3	!! * MODULE obs_oper *
4	!! Observation diagnostics: Observation operators for various observation
5	!! types
6	!!======================================================================
7
8	!!----------------------------------------------------------------------
9	!! obs_pro_opt : Compute the model counterpart of temperature and
10	!! salinity observations from profiles
11	!! obs_sla_opt : Compute the model counterpart of sea level anomaly
12	!! observations
13	!! obs_sst_opt : Compute the model counterpart of sea surface temperature
14	!! observations
15	!! obs_sss_opt : Compute the model counterpart of sea surface salinity
16	!! observations
17	!! obs_seaice_opt : Compute the model counterpart of sea ice concentration
18	!! observations
19	!!
20	!! obs_vel_opt : Compute the model counterpart of zonal and meridional
21	!! components of velocity from observations.
22	!!----------------------------------------------------------------------
23
24	!! * Modules used
25	USE par_kind, ONLY : & ! Precision variables
26	& wp
27	USE in_out_manager ! I/O manager
28	USE obs_inter_sup ! Interpolation support
29	USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the observation pt
30	& obs_int_h2d, &
31	& obs_int_h2d_init
32	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt
33	& obs_int_z1d, &
34	& obs_int_z1d_spl
35	USE obs_const, ONLY : &
36	& obfillflt ! Fillvalue
37	USE dom_oce, ONLY : &
38	& glamt, glamu, glamv, &
39	& gphit, gphiu, gphiv
40
41	IMPLICIT NONE
42
43	!! * Routine accessibility
44	PRIVATE
45
46	PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations
47	& obs_sla_opt, & ! Compute the model counterpart of SLA observations
48	& obs_sst_opt, & ! Compute the model counterpart of SST observations
49	& obs_sss_opt, & ! Compute the model counterpart of SSS observations
50	& obs_seaice_opt, &
51	& obs_vel_opt ! Compute the model counterpart of velocity profile data
52
53	INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
54
55	CONTAINS
56
57	SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
58	& ptn, psn, pgdept, ptmask, k1dint, k2dint, &
59	& kdailyavtypes )
60	!!-----------------------------------------------------------------------
61	!!
62	!! * ROUTINE obs_pro_opt *
63	!!
64	!! ** Purpose : Compute the model counterpart of profiles
65	!! data by interpolating from the model grid to the
66	!! observation point.
67	!!
68	!! ** Method : Linearly interpolate to each observation point using
69	!! the model values at the corners of the surrounding grid box.
70	!!
71	!! First, a vertical profile of horizontally interpolated model
72	!! now temperatures is computed at the obs (lon, lat) point.
73	!! Several horizontal interpolation schemes are available:
74	!! - distance-weighted (great circle) (k2dint = 0)
75	!! - distance-weighted (small angle) (k2dint = 1)
76	!! - bilinear (geographical grid) (k2dint = 2)
77	!! - bilinear (quadrilateral grid) (k2dint = 3)
78	!! - polynomial (quadrilateral grid) (k2dint = 4)
79	!!
80	!! Next, the vertical temperature profile is interpolated to the
81	!! data depth points. Two vertical interpolation schemes are
82	!! available:
83	!! - linear (k1dint = 0)
84	!! - Cubic spline (k1dint = 1)
85	!!
86	!! For the cubic spline the 2nd derivative of the interpolating
87	!! polynomial is computed before entering the vertical interpolation
88	!! routine.
89	!!
90	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
91	!! a daily mean model temperature field. So, we first compute
92	!! the mean, then interpolate only at the end of the day.
93	!!
94	!! Note: the in situ temperature observations must be converted
95	!! to potential temperature (the model variable) prior to
96	!! assimilation.
97	!!??????????????????????????????????????????????????????????????
98	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
99	!!??????????????????????????????????????????????????????????????
100	!!
101	!! ** Action :
102	!!
103	!! History :
104	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
105	!! ! 06-03 (G. Smith) NEMOVAR migration
106	!! ! 06-10 (A. Weaver) Cleanup
107	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
108	!! ! 07-03 (K. Mogensen) General handling of profiles
109	!!-----------------------------------------------------------------------
110
111	!! * Modules used
112	USE obs_profiles_def ! Definition of storage space for profile obs.
113
114	IMPLICIT NONE
115
116	!! * Arguments
117	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
118	INTEGER, INTENT(IN) :: kt ! Time step
119	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
120	INTEGER, INTENT(IN) :: kpj
121	INTEGER, INTENT(IN) :: kpk
122	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
123	! (kit000-1 = restart time)
124	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
125	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
126	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
127	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
128	& ptn, & ! Model temperature field
129	& psn, & ! Model salinity field
130	& ptmask ! Land-sea mask
131	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
132	& pgdept ! Model array of depth levels
133	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
134	& kdailyavtypes! Types for daily averages
135	!! * Local declarations
136	INTEGER :: ji
137	INTEGER :: jj
138	INTEGER :: jk
139	INTEGER :: jobs
140	INTEGER :: inrc
141	INTEGER :: ipro
142	INTEGER :: idayend
143	INTEGER :: ista
144	INTEGER :: iend
145	INTEGER :: iobs
146	INTEGER, DIMENSION(imaxavtypes) :: &
147	& idailyavtypes
148	REAL(KIND=wp) :: zlam
149	REAL(KIND=wp) :: zphi
150	REAL(KIND=wp) :: zdaystp
151	REAL(KIND=wp), DIMENSION(kpk) :: &
152	& zobsmask, &
153	& zobsk, &
154	& zobs2k
155	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
156	& zweig
157	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
158	& zmask, &
159	& zintt, &
160	& zints, &
161	& zinmt, &
162	& zinms
163	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
164	& zglam, &
165	& zgphi
166	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
167	& igrdi, &
168	& igrdj
169
170	!------------------------------------------------------------------------
171	! Local initialization
172	!------------------------------------------------------------------------
173	! ... Record and data counters
174	inrc = kt - kit000 + 2
175	ipro = prodatqc%npstp(inrc)
176
177	! Daily average types
178	IF ( PRESENT(kdailyavtypes) ) THEN
179	idailyavtypes(:) = kdailyavtypes(:)
180	ELSE
181	idailyavtypes(:) = -1
182	ENDIF
183
184	! Initialize daily mean for first timestep
185	idayend = MOD( kt - kit000 + 1, kdaystp )
186
187	! Added kt == 0 test to catch restart case
188	IF ( idayend == 1 .OR. kt == 0) THEN
189	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
190	DO jk = 1, jpk
191	DO jj = 1, jpj
192	DO ji = 1, jpi
193	prodatqc%vdmean(ji,jj,jk,1) = 0.0
194	prodatqc%vdmean(ji,jj,jk,2) = 0.0
195	END DO
196	END DO
197	END DO
198	ENDIF
199
200	DO jk = 1, jpk
201	DO jj = 1, jpj
202	DO ji = 1, jpi
203	! Increment the temperature field for computing daily mean
204	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
205	& + ptn(ji,jj,jk)
206	! Increment the salinity field for computing daily mean
207	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
208	& + psn(ji,jj,jk)
209	END DO
210	END DO
211	END DO
212
213	! Compute the daily mean at the end of day
214	zdaystp = 1.0 / REAL( kdaystp )
215	IF ( idayend == 0 ) THEN
216	DO jk = 1, jpk
217	DO jj = 1, jpj
218	DO ji = 1, jpi
219	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
220	& * zdaystp
221	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
222	& * zdaystp
223	END DO
224	END DO
225	END DO
226	ENDIF
227
228	! Get the data for interpolation
229	ALLOCATE( &
230	& igrdi(2,2,ipro), &
231	& igrdj(2,2,ipro), &
232	& zglam(2,2,ipro), &
233	& zgphi(2,2,ipro), &
234	& zmask(2,2,kpk,ipro), &
235	& zintt(2,2,kpk,ipro), &
236	& zints(2,2,kpk,ipro) &
237	& )
238
239	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
240	iobs = jobs - prodatqc%nprofup
241	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
242	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
243	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
244	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
245	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
246	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
247	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
248	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
249	END DO
250
251	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
252	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
253	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
254	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
255	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
256
257	! At the end of the day also get interpolated means
258	IF ( idayend == 0 ) THEN
259
260	ALLOCATE( &
261	& zinmt(2,2,kpk,ipro), &
262	& zinms(2,2,kpk,ipro) &
263	& )
264
265	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
266	& prodatqc%vdmean(:,:,:,1), zinmt )
267	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
268	& prodatqc%vdmean(:,:,:,2), zinms )
269
270	ENDIF
271
272	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
273
274	iobs = jobs - prodatqc%nprofup
275
276	IF ( kt /= prodatqc%mstp(jobs) ) THEN
277
278	IF(lwp) THEN
279	WRITE(numout,*)
280	WRITE(numout,*) ' E R R O R : Observation', &
281	& ' time step is not consistent with the', &
282	& ' model time step'
283	WRITE(numout,*) ' ========='
284	WRITE(numout,*)
285	WRITE(numout,*) ' Record = ', jobs, &
286	& ' kt = ', kt, &
287	& ' mstp = ', prodatqc%mstp(jobs), &
288	& ' ntyp = ', prodatqc%ntyp(jobs)
289	ENDIF
290	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
291	ENDIF
292
293	zlam = prodatqc%rlam(jobs)
294	zphi = prodatqc%rphi(jobs)
295
296	! Horizontal weights and vertical mask
297
298	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
299	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
300
301	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
302	& zglam(:,:,iobs), zgphi(:,:,iobs), &
303	& zmask(:,:,:,iobs), zweig, zobsmask )
304
305	ENDIF
306
307	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
308
309	zobsk(:) = obfillflt
310
311	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
312
313	IF ( idayend == 0 ) THEN
314
315	! Daily averaged moored buoy (MRB) data
316
317	CALL obs_int_h2d( kpk, kpk, &
318	& zweig, zinmt(:,:,:,iobs), zobsk )
319
320
321	ELSE
322
323	CALL ctl_stop( ' A nonzero' // &
324	& ' number of profile T BUOY data should' // &
325	& ' only occur at the end of a given day' )
326
327	ENDIF
328
329	ELSE
330
331	! Point data
332
333	CALL obs_int_h2d( kpk, kpk, &
334	& zweig, zintt(:,:,:,iobs), zobsk )
335
336	ENDIF
337
338	!-------------------------------------------------------------
339	! Compute vertical second-derivative of the interpolating
340	! polynomial at obs points
341	!-------------------------------------------------------------
342
343	IF ( k1dint == 1 ) THEN
344	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
345	& pgdept, zobsmask )
346	ENDIF
347
348	!-----------------------------------------------------------------
349	! Vertical interpolation to the observation point
350	!-----------------------------------------------------------------
351	ista = prodatqc%npvsta(jobs,1)
352	iend = prodatqc%npvend(jobs,1)
353	CALL obs_int_z1d( kpk, &
354	& prodatqc%var(1)%mvk(ista:iend), &
355	& k1dint, iend - ista + 1, &
356	& prodatqc%var(1)%vdep(ista:iend), &
357	& zobsk, zobs2k, &
358	& prodatqc%var(1)%vmod(ista:iend), &
359	& pgdept, zobsmask )
360
361	ENDIF
362
363	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
364
365	zobsk(:) = obfillflt
366
367	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
368
369	IF ( idayend == 0 ) THEN
370
371	! Daily averaged moored buoy (MRB) data
372
373	CALL obs_int_h2d( kpk, kpk, &
374	& zweig, zinms(:,:,:,iobs), zobsk )
375
376	ELSE
377
378	CALL ctl_stop( ' A nonzero' // &
379	& ' number of profile S BUOY data should' // &
380	& ' only occur at the end of a given day' )
381
382	ENDIF
383
384	ELSE
385
386	! Point data
387
388	CALL obs_int_h2d( kpk, kpk, &
389	& zweig, zints(:,:,:,iobs), zobsk )
390
391	ENDIF
392
393
394	!-------------------------------------------------------------
395	! Compute vertical second-derivative of the interpolating
396	! polynomial at obs points
397	!-------------------------------------------------------------
398
399	IF ( k1dint == 1 ) THEN
400	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
401	& pgdept, zobsmask )
402	ENDIF
403
404	!----------------------------------------------------------------
405	! Vertical interpolation to the observation point
406	!----------------------------------------------------------------
407	ista = prodatqc%npvsta(jobs,2)
408	iend = prodatqc%npvend(jobs,2)
409	CALL obs_int_z1d( kpk, &
410	& prodatqc%var(2)%mvk(ista:iend),&
411	& k1dint, iend - ista + 1, &
412	& prodatqc%var(2)%vdep(ista:iend),&
413	& zobsk, zobs2k, &
414	& prodatqc%var(2)%vmod(ista:iend),&
415	& pgdept, zobsmask )
416
417	ENDIF
418
419	END DO
420
421	! Deallocate the data for interpolation
422	DEALLOCATE( &
423	& igrdi, &
424	& igrdj, &
425	& zglam, &
426	& zgphi, &
427	& zmask, &
428	& zintt, &
429	& zints &
430	& )
431	! At the end of the day also get interpolated means
432	IF ( idayend == 0 ) THEN
433	DEALLOCATE( &
434	& zinmt, &
435	& zinms &
436	& )
437	ENDIF
438
439	prodatqc%nprofup = prodatqc%nprofup + ipro
440
441	END SUBROUTINE obs_pro_opt
442
443	SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
444	& psshn, psshmask, k2dint )
445	!!-----------------------------------------------------------------------
446	!!
447	!! * ROUTINE obs_sla_opt *
448	!!
449	!! ** Purpose : Compute the model counterpart of sea level anomaly
450	!! data by interpolating from the model grid to the
451	!! observation point.
452	!!
453	!! ** Method : Linearly interpolate to each observation point using
454	!! the model values at the corners of the surrounding grid box.
455	!!
456	!! The now model SSH is first computed at the obs (lon, lat) point.
457	!!
458	!! Several horizontal interpolation schemes are available:
459	!! - distance-weighted (great circle) (k2dint = 0)
460	!! - distance-weighted (small angle) (k2dint = 1)
461	!! - bilinear (geographical grid) (k2dint = 2)
462	!! - bilinear (quadrilateral grid) (k2dint = 3)
463	!! - polynomial (quadrilateral grid) (k2dint = 4)
464	!!
465	!! The sea level anomaly at the observation points is then computed
466	!! by removing a mean dynamic topography (defined at the obs. point).
467	!!
468	!! ** Action :
469	!!
470	!! History :
471	!! ! 07-03 (A. Weaver)
472	!!-----------------------------------------------------------------------
473
474	!! * Modules used
475	USE obs_surf_def ! Definition of storage space for surface observations
476
477	IMPLICIT NONE
478
479	!! * Arguments
480	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening
481	INTEGER, INTENT(IN) :: kt ! Time step
482	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
483	INTEGER, INTENT(IN) :: kpj
484	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
485	! (kit000-1 = restart time)
486	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
487	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
488	& psshn, & ! Model SSH field
489	& psshmask ! Land-sea mask
490
491	!! * Local declarations
492	INTEGER :: ji
493	INTEGER :: jj
494	INTEGER :: jobs
495	INTEGER :: inrc
496	INTEGER :: isla
497	INTEGER :: iobs
498	REAL(KIND=wp) :: zlam
499	REAL(KIND=wp) :: zphi
500	REAL(KIND=wp) :: zext(1), zobsmask(1)
501	REAL(kind=wp), DIMENSION(2,2,1) :: &
502	& zweig
503	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
504	& zmask, &
505	& zsshl, &
506	& zglam, &
507	& zgphi
508	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
509	& igrdi, &
510	& igrdj
511
512	!------------------------------------------------------------------------
513	! Local initialization
514	!------------------------------------------------------------------------
515	! ... Record and data counters
516	inrc = kt - kit000 + 2
517	isla = sladatqc%nsstp(inrc)
518
519	! Get the data for interpolation
520
521	ALLOCATE( &
522	& igrdi(2,2,isla), &
523	& igrdj(2,2,isla), &
524	& zglam(2,2,isla), &
525	& zgphi(2,2,isla), &
526	& zmask(2,2,isla), &
527	& zsshl(2,2,isla) &
528	& )
529
530	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
531	iobs = jobs - sladatqc%nsurfup
532	igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
533	igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
534	igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
535	igrdj(1,2,iobs) = sladatqc%mj(jobs)
536	igrdi(2,1,iobs) = sladatqc%mi(jobs)
537	igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
538	igrdi(2,2,iobs) = sladatqc%mi(jobs)
539	igrdj(2,2,iobs) = sladatqc%mj(jobs)
540	END DO
541
542	CALL obs_int_comm_2d( 2, 2, isla, &
543	& igrdi, igrdj, glamt, zglam )
544	CALL obs_int_comm_2d( 2, 2, isla, &
545	& igrdi, igrdj, gphit, zgphi )
546	CALL obs_int_comm_2d( 2, 2, isla, &
547	& igrdi, igrdj, psshmask, zmask )
548	CALL obs_int_comm_2d( 2, 2, isla, &
549	& igrdi, igrdj, psshn, zsshl )
550
551	! Loop over observations
552
553	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
554
555	iobs = jobs - sladatqc%nsurfup
556
557	IF ( kt /= sladatqc%mstp(jobs) ) THEN
558
559	IF(lwp) THEN
560	WRITE(numout,*)
561	WRITE(numout,*) ' E R R O R : Observation', &
562	& ' time step is not consistent with the', &
563	& ' model time step'
564	WRITE(numout,*) ' ========='
565	WRITE(numout,*)
566	WRITE(numout,*) ' Record = ', jobs, &
567	& ' kt = ', kt, &
568	& ' mstp = ', sladatqc%mstp(jobs), &
569	& ' ntyp = ', sladatqc%ntyp(jobs)
570	ENDIF
571	CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
572
573	ENDIF
574
575	zlam = sladatqc%rlam(jobs)
576	zphi = sladatqc%rphi(jobs)
577
578	! Get weights to interpolate the model SSH to the observation point
579	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
580	& zglam(:,:,iobs), zgphi(:,:,iobs), &
581	& zmask(:,:,iobs), zweig, zobsmask )
582
583
584	! Interpolate the model SSH to the observation point
585	CALL obs_int_h2d( 1, 1, &
586	& zweig, zsshl(:,:,iobs), zext )
587
588	sladatqc%rext(jobs,1) = zext(1)
589	! ... Remove the MDT at the observation point
590	sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
591
592	END DO
593
594	! Deallocate the data for interpolation
595	DEALLOCATE( &
596	& igrdi, &
597	& igrdj, &
598	& zglam, &
599	& zgphi, &
600	& zmask, &
601	& zsshl &
602	& )
603
604	sladatqc%nsurfup = sladatqc%nsurfup + isla
605
606	END SUBROUTINE obs_sla_opt
607
608	SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, &
609	& psstn, psstmask, k2dint )
610
611	!!-----------------------------------------------------------------------
612	!!
613	!! * ROUTINE obs_sst_opt *
614	!!
615	!! ** Purpose : Compute the model counterpart of surface temperature
616	!! data by interpolating from the model grid to the
617	!! observation point.
618	!!
619	!! ** Method : Linearly interpolate to each observation point using
620	!! the model values at the corners of the surrounding grid box.
621	!!
622	!! The now model SST is first computed at the obs (lon, lat) point.
623	!!
624	!! Several horizontal interpolation schemes are available:
625	!! - distance-weighted (great circle) (k2dint = 0)
626	!! - distance-weighted (small angle) (k2dint = 1)
627	!! - bilinear (geographical grid) (k2dint = 2)
628	!! - bilinear (quadrilateral grid) (k2dint = 3)
629	!! - polynomial (quadrilateral grid) (k2dint = 4)
630	!!
631	!!
632	!! ** Action :
633	!!
634	!! History :
635	!! ! 07-07 (S. Ricci ) : Original
636	!!
637	!!-----------------------------------------------------------------------
638
639	!! * Modules used
640	USE obs_surf_def ! Definition of storage space for surface observations
641
642	IMPLICIT NONE
643
644	!! * Arguments
645	TYPE(obs_surf), INTENT(INOUT) :: &
646	& sstdatqc ! Subset of surface data not failing screening
647	INTEGER, INTENT(IN) :: kt ! Time step
648	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
649	INTEGER, INTENT(IN) :: kpj
650	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
651	! (kit000-1 = restart time)
652	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
653	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
654	& psstn, & ! Model SST field
655	& psstmask ! Land-sea mask
656
657	!! * Local declarations
658	INTEGER :: ji
659	INTEGER :: jj
660	INTEGER :: jobs
661	INTEGER :: inrc
662	INTEGER :: isst
663	INTEGER :: iobs
664	REAL(KIND=wp) :: zlam
665	REAL(KIND=wp) :: zphi
666	REAL(KIND=wp) :: zext(1), zobsmask(1)
667	REAL(kind=wp), DIMENSION(2,2,1) :: &
668	& zweig
669	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
670	& zmask, &
671	& zsstl, &
672	& zglam, &
673	& zgphi
674	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
675	& igrdi, &
676	& igrdj
677
678	!-----------------------------------------------------------------------
679	! Local initialization
680	!-----------------------------------------------------------------------
681	! ... Record and data counters
682	inrc = kt - kit000 + 2
683	isst = sstdatqc%nsstp(inrc)
684
685	! Get the data for interpolation
686
687	ALLOCATE( &
688	& igrdi(2,2,isst), &
689	& igrdj(2,2,isst), &
690	& zglam(2,2,isst), &
691	& zgphi(2,2,isst), &
692	& zmask(2,2,isst), &
693	& zsstl(2,2,isst) &
694	& )
695
696	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
697	iobs = jobs - sstdatqc%nsurfup
698	igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
699	igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
700	igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
701	igrdj(1,2,iobs) = sstdatqc%mj(jobs)
702	igrdi(2,1,iobs) = sstdatqc%mi(jobs)
703	igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
704	igrdi(2,2,iobs) = sstdatqc%mi(jobs)
705	igrdj(2,2,iobs) = sstdatqc%mj(jobs)
706	END DO
707
708	CALL obs_int_comm_2d( 2, 2, isst, &
709	& igrdi, igrdj, glamt, zglam )
710	CALL obs_int_comm_2d( 2, 2, isst, &
711	& igrdi, igrdj, gphit, zgphi )
712	CALL obs_int_comm_2d( 2, 2, isst, &
713	& igrdi, igrdj, psstmask, zmask )
714	CALL obs_int_comm_2d( 2, 2, isst, &
715	& igrdi, igrdj, psstn, zsstl )
716
717	! Loop over observations
718
719	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
720
721	iobs = jobs - sstdatqc%nsurfup
722
723	IF ( kt /= sstdatqc%mstp(jobs) ) THEN
724
725	IF(lwp) THEN
726	WRITE(numout,*)
727	WRITE(numout,*) ' E R R O R : Observation', &
728	& ' time step is not consistent with the', &
729	& ' model time step'
730	WRITE(numout,*) ' ========='
731	WRITE(numout,*)
732	WRITE(numout,*) ' Record = ', jobs, &
733	& ' kt = ', kt, &
734	& ' mstp = ', sstdatqc%mstp(jobs), &
735	& ' ntyp = ', sstdatqc%ntyp(jobs)
736	ENDIF
737	CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
738
739	ENDIF
740
741	zlam = sstdatqc%rlam(jobs)
742	zphi = sstdatqc%rphi(jobs)
743
744	! Get weights to interpolate the model SST to the observation point
745	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
746	& zglam(:,:,iobs), zgphi(:,:,iobs), &
747	& zmask(:,:,iobs), zweig, zobsmask )
748
749	! Interpolate the model SST to the observation point
750	CALL obs_int_h2d( 1, 1, &
751	& zweig, zsstl(:,:,iobs), zext )
752
753	sstdatqc%rmod(jobs,1) = zext(1)
754
755	END DO
756
757	! Deallocate the data for interpolation
758	DEALLOCATE( &
759	& igrdi, &
760	& igrdj, &
761	& zglam, &
762	& zgphi, &
763	& zmask, &
764	& zsstl &
765	& )
766
767	sstdatqc%nsurfup = sstdatqc%nsurfup + isst
768
769	END SUBROUTINE obs_sst_opt
770
771	SUBROUTINE obs_sss_opt
772	!!-----------------------------------------------------------------------
773	!!
774	!! * ROUTINE obs_sss_opt *
775	!!
776	!! ** Purpose : Compute the model counterpart of sea surface salinity
777	!! data by interpolating from the model grid to the
778	!! observation point.
779	!!
780	!! ** Method :
781	!!
782	!! ** Action :
783	!!
784	!! History :
785	!! ! ??-??
786	!!-----------------------------------------------------------------------
787
788	IMPLICIT NONE
789
790	END SUBROUTINE obs_sss_opt
791
792	SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
793	& pseaicen, pseaicemask, k2dint )
794
795	!!-----------------------------------------------------------------------
796	!!
797	!! * ROUTINE obs_seaice_opt *
798	!!
799	!! ** Purpose : Compute the model counterpart of surface temperature
800	!! data by interpolating from the model grid to the
801	!! observation point.
802	!!
803	!! ** Method : Linearly interpolate to each observation point using
804	!! the model values at the corners of the surrounding grid box.
805	!!
806	!! The now model sea ice is first computed at the obs (lon, lat) point.
807	!!
808	!! Several horizontal interpolation schemes are available:
809	!! - distance-weighted (great circle) (k2dint = 0)
810	!! - distance-weighted (small angle) (k2dint = 1)
811	!! - bilinear (geographical grid) (k2dint = 2)
812	!! - bilinear (quadrilateral grid) (k2dint = 3)
813	!! - polynomial (quadrilateral grid) (k2dint = 4)
814	!!
815	!!
816	!! ** Action :
817	!!
818	!! History :
819	!! ! 07-07 (S. Ricci ) : Original
820	!!
821	!!-----------------------------------------------------------------------
822
823	!! * Modules used
824	USE obs_surf_def ! Definition of storage space for surface observations
825
826	IMPLICIT NONE
827
828	!! * Arguments
829	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening
830	INTEGER, INTENT(IN) :: kt ! Time step
831	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
832	INTEGER, INTENT(IN) :: kpj
833	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
834	! (kit000-1 = restart time)
835	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
836	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
837	& pseaicen, & ! Model sea ice field
838	& pseaicemask ! Land-sea mask
839
840	!! * Local declarations
841	INTEGER :: ji
842	INTEGER :: jj
843	INTEGER :: jobs
844	INTEGER :: inrc
845	INTEGER :: iseaice
846	INTEGER :: iobs
847
848	REAL(KIND=wp) :: zlam
849	REAL(KIND=wp) :: zphi
850	REAL(KIND=wp) :: zext(1), zobsmask(1)
851	REAL(kind=wp), DIMENSION(2,2,1) :: &
852	& zweig
853	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
854	& zmask, &
855	& zseaicel, &
856	& zglam, &
857	& zgphi
858	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
859	& igrdi, &
860	& igrdj
861
862	!------------------------------------------------------------------------
863	! Local initialization
864	!------------------------------------------------------------------------
865	! ... Record and data counters
866	inrc = kt - kit000 + 2
867	iseaice = seaicedatqc%nsstp(inrc)
868
869	! Get the data for interpolation
870
871	ALLOCATE( &
872	& igrdi(2,2,iseaice), &
873	& igrdj(2,2,iseaice), &
874	& zglam(2,2,iseaice), &
875	& zgphi(2,2,iseaice), &
876	& zmask(2,2,iseaice), &
877	& zseaicel(2,2,iseaice) &
878	& )
879
880	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
881	iobs = jobs - seaicedatqc%nsurfup
882	igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
883	igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
884	igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
885	igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
886	igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
887	igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
888	igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
889	igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
890	END DO
891
892	CALL obs_int_comm_2d( 2, 2, iseaice, &
893	& igrdi, igrdj, glamt, zglam )
894	CALL obs_int_comm_2d( 2, 2, iseaice, &
895	& igrdi, igrdj, gphit, zgphi )
896	CALL obs_int_comm_2d( 2, 2, iseaice, &
897	& igrdi, igrdj, pseaicemask, zmask )
898	CALL obs_int_comm_2d( 2, 2, iseaice, &
899	& igrdi, igrdj, pseaicen, zseaicel )
900
901	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
902
903	iobs = jobs - seaicedatqc%nsurfup
904
905	IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
906
907	IF(lwp) THEN
908	WRITE(numout,*)
909	WRITE(numout,*) ' E R R O R : Observation', &
910	& ' time step is not consistent with the', &
911	& ' model time step'
912	WRITE(numout,*) ' ========='
913	WRITE(numout,*)
914	WRITE(numout,*) ' Record = ', jobs, &
915	& ' kt = ', kt, &
916	& ' mstp = ', seaicedatqc%mstp(jobs), &
917	& ' ntyp = ', seaicedatqc%ntyp(jobs)
918	ENDIF
919	CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
920
921	ENDIF
922
923	zlam = seaicedatqc%rlam(jobs)
924	zphi = seaicedatqc%rphi(jobs)
925
926	! Get weights to interpolate the model sea ice to the observation point
927	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
928	& zglam(:,:,iobs), zgphi(:,:,iobs), &
929	& zmask(:,:,iobs), zweig, zobsmask )
930
931	! ... Interpolate the model sea ice to the observation point
932	CALL obs_int_h2d( 1, 1, &
933	& zweig, zseaicel(:,:,iobs), zext )
934
935	seaicedatqc%rmod(jobs,1) = zext(1)
936
937	END DO
938
939	! Deallocate the data for interpolation
940	DEALLOCATE( &
941	& igrdi, &
942	& igrdj, &
943	& zglam, &
944	& zgphi, &
945	& zmask, &
946	& zseaicel &
947	& )
948
949	seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
950
951	END SUBROUTINE obs_seaice_opt
952
953	SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
954	& pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
955	& ld_dailyav )
956	!!-----------------------------------------------------------------------
957	!!
958	!! * ROUTINE obs_vel_opt *
959	!!
960	!! ** Purpose : Compute the model counterpart of velocity profile
961	!! data by interpolating from the model grid to the
962	!! observation point.
963	!!
964	!! ** Method : Linearly interpolate zonal and meridional components of velocity
965	!! to each observation point using the model values at the corners of
966	!! the surrounding grid box. The model velocity components are on a
967	!! staggered C- grid.
968	!!
969	!! For velocity data from the TAO array, the model equivalent is
970	!! a daily mean velocity field. So, we first compute
971	!! the mean, then interpolate only at the end of the day.
972	!!
973	!! ** Action :
974	!!
975	!! History :
976	!! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles
977	!! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
978	!!-----------------------------------------------------------------------
979
980	!! * Modules used
981	USE obs_profiles_def ! Definition of storage space for profile obs.
982
983	IMPLICIT NONE
984
985	!! * Arguments
986	TYPE(obs_prof), INTENT(INOUT) :: &
987	& prodatqc ! Subset of profile data not failing screening
988	INTEGER, INTENT(IN) :: kt ! Time step
989	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
990	INTEGER, INTENT(IN) :: kpj
991	INTEGER, INTENT(IN) :: kpk
992	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
993	! (kit000-1 = restart time)
994	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
995	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
996	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
997	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
998	& pun, & ! Model zonal component of velocity
999	& pvn, & ! Model meridional component of velocity
1000	& pumask, & ! Land-sea mask
1001	& pvmask ! Land-sea mask
1002	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1003	& pgdept ! Model array of depth levels
1004	LOGICAL, INTENT(IN) :: ld_dailyav
1005
1006	!! * Local declarations
1007	INTEGER :: ji
1008	INTEGER :: jj
1009	INTEGER :: jk
1010	INTEGER :: jobs
1011	INTEGER :: inrc
1012	INTEGER :: ipro
1013	INTEGER :: idayend
1014	INTEGER :: ista
1015	INTEGER :: iend
1016	INTEGER :: iobs
1017	INTEGER, DIMENSION(imaxavtypes) :: &
1018	& idailyavtypes
1019	REAL(KIND=wp) :: zlam
1020	REAL(KIND=wp) :: zphi
1021	REAL(KIND=wp) :: zdaystp
1022	REAL(KIND=wp), DIMENSION(kpk) :: &
1023	& zobsmask, &
1024	& zobsk, &
1025	& zobs2k
1026	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1027	& zweigu,zweigv
1028	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1029	& zumask, zvmask, &
1030	& zintu, &
1031	& zintv, &
1032	& zinmu, &
1033	& zinmv
1034	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1035	& zglamu, zglamv, &
1036	& zgphiu, zgphiv
1037	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1038	& igrdiu, &
1039	& igrdju, &
1040	& igrdiv, &
1041	& igrdjv
1042
1043	!------------------------------------------------------------------------
1044	! Local initialization
1045	!------------------------------------------------------------------------
1046	! ... Record and data counters
1047	inrc = kt - kit000 + 2
1048	ipro = prodatqc%npstp(inrc)
1049
1050	! Initialize daily mean for first timestep
1051	idayend = MOD( kt - kit000 + 1, kdaystp )
1052
1053	! Added kt == 0 test to catch restart case
1054	IF ( idayend == 1 .OR. kt == 0) THEN
1055	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1056	prodatqc%vdmean(:,:,:,1) = 0.0
1057	prodatqc%vdmean(:,:,:,2) = 0.0
1058	ENDIF
1059
1060	! Increment the zonal velocity field for computing daily mean
1061	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1062	! Increment the meridional velocity field for computing daily mean
1063	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1064
1065	! Compute the daily mean at the end of day
1066	zdaystp = 1.0 / REAL( kdaystp )
1067	IF ( idayend == 0 ) THEN
1068	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1069	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1070	ENDIF
1071
1072	! Get the data for interpolation
1073	ALLOCATE( &
1074	& igrdiu(2,2,ipro), &
1075	& igrdju(2,2,ipro), &
1076	& igrdiv(2,2,ipro), &
1077	& igrdjv(2,2,ipro), &
1078	& zglamu(2,2,ipro), zglamv(2,2,ipro), &
1079	& zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1080	& zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1081	& zintu(2,2,kpk,ipro), &
1082	& zintv(2,2,kpk,ipro) &
1083	& )
1084
1085	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1086	iobs = jobs - prodatqc%nprofup
1087	igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1088	igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1089	igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1090	igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1091	igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1092	igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1093	igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1094	igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1095	igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1096	igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1097	igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1098	igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1099	igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1100	igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1101	igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1102	igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1103	END DO
1104
1105	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1106	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1107	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1108	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1109
1110	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1111	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1112	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1113	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1114
1115	! At the end of the day also get interpolated means
1116	IF ( idayend == 0 ) THEN
1117
1118	ALLOCATE( &
1119	& zinmu(2,2,kpk,ipro), &
1120	& zinmv(2,2,kpk,ipro) &
1121	& )
1122
1123	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1124	& prodatqc%vdmean(:,:,:,1), zinmu )
1125	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1126	& prodatqc%vdmean(:,:,:,2), zinmv )
1127
1128	ENDIF
1129
1130	! loop over observations
1131
1132	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1133
1134	iobs = jobs - prodatqc%nprofup
1135
1136	IF ( kt /= prodatqc%mstp(jobs) ) THEN
1137
1138	IF(lwp) THEN
1139	WRITE(numout,*)
1140	WRITE(numout,*) ' E R R O R : Observation', &
1141	& ' time step is not consistent with the', &
1142	& ' model time step'
1143	WRITE(numout,*) ' ========='
1144	WRITE(numout,*)
1145	WRITE(numout,*) ' Record = ', jobs, &
1146	& ' kt = ', kt, &
1147	& ' mstp = ', prodatqc%mstp(jobs), &
1148	& ' ntyp = ', prodatqc%ntyp(jobs)
1149	ENDIF
1150	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1151	ENDIF
1152
1153	zlam = prodatqc%rlam(jobs)
1154	zphi = prodatqc%rphi(jobs)
1155
1156	! Horizontal weights and vertical mask
1157
1158	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1159
1160	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1161	& zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1162	& zumask(:,:,:,iobs), zweigu, zobsmask )
1163
1164	ENDIF
1165
1166
1167	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1168
1169	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1170	& zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1171	& zvmask(:,:,:,iobs), zweigv, zobsmask )
1172
1173	ENDIF
1174
1175
1176	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1177
1178	zobsk(:) = obfillflt
1179
1180	IF ( ld_dailyav ) THEN
1181
1182	IF ( idayend == 0 ) THEN
1183
1184	! Daily averaged data
1185
1186	CALL obs_int_h2d( kpk, kpk, &
1187	& zweigu, zinmu(:,:,:,iobs), zobsk )
1188
1189
1190	ELSE
1191
1192	CALL ctl_stop( ' A nonzero' // &
1193	& ' number of U profile data should' // &
1194	& ' only occur at the end of a given day' )
1195
1196	ENDIF
1197
1198	ELSE
1199
1200	! Point data
1201
1202	CALL obs_int_h2d( kpk, kpk, &
1203	& zweigu, zintu(:,:,:,iobs), zobsk )
1204
1205	ENDIF
1206
1207	!-------------------------------------------------------------
1208	! Compute vertical second-derivative of the interpolating
1209	! polynomial at obs points
1210	!-------------------------------------------------------------
1211
1212	IF ( k1dint == 1 ) THEN
1213	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1214	& pgdept, zobsmask )
1215	ENDIF
1216
1217	!-----------------------------------------------------------------
1218	! Vertical interpolation to the observation point
1219	!-----------------------------------------------------------------
1220	ista = prodatqc%npvsta(jobs,1)
1221	iend = prodatqc%npvend(jobs,1)
1222	CALL obs_int_z1d( kpk, &
1223	& prodatqc%var(1)%mvk(ista:iend), &
1224	& k1dint, iend - ista + 1, &
1225	& prodatqc%var(1)%vdep(ista:iend), &
1226	& zobsk, zobs2k, &
1227	& prodatqc%var(1)%vmod(ista:iend), &
1228	& pgdept, zobsmask )
1229
1230	ENDIF
1231
1232	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1233
1234	zobsk(:) = obfillflt
1235
1236	IF ( ld_dailyav ) THEN
1237
1238	IF ( idayend == 0 ) THEN
1239
1240	! Daily averaged data
1241
1242	CALL obs_int_h2d( kpk, kpk, &
1243	& zweigv, zinmv(:,:,:,iobs), zobsk )
1244
1245	ELSE
1246
1247	CALL ctl_stop( ' A nonzero' // &
1248	& ' number of V profile data should' // &
1249	& ' only occur at the end of a given day' )
1250
1251	ENDIF
1252
1253	ELSE
1254
1255	! Point data
1256
1257	CALL obs_int_h2d( kpk, kpk, &
1258	& zweigv, zintv(:,:,:,iobs), zobsk )
1259
1260	ENDIF
1261
1262
1263	!-------------------------------------------------------------
1264	! Compute vertical second-derivative of the interpolating
1265	! polynomial at obs points
1266	!-------------------------------------------------------------
1267
1268	IF ( k1dint == 1 ) THEN
1269	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1270	& pgdept, zobsmask )
1271	ENDIF
1272
1273	!----------------------------------------------------------------
1274	! Vertical interpolation to the observation point
1275	!----------------------------------------------------------------
1276	ista = prodatqc%npvsta(jobs,2)
1277	iend = prodatqc%npvend(jobs,2)
1278	CALL obs_int_z1d( kpk, &
1279	& prodatqc%var(2)%mvk(ista:iend),&
1280	& k1dint, iend - ista + 1, &
1281	& prodatqc%var(2)%vdep(ista:iend),&
1282	& zobsk, zobs2k, &
1283	& prodatqc%var(2)%vmod(ista:iend),&
1284	& pgdept, zobsmask )
1285
1286	ENDIF
1287
1288	END DO
1289
1290	! Deallocate the data for interpolation
1291	DEALLOCATE( &
1292	& igrdiu, &
1293	& igrdju, &
1294	& igrdiv, &
1295	& igrdjv, &
1296	& zglamu, zglamv, &
1297	& zgphiu, zgphiv, &
1298	& zumask, zvmask, &
1299	& zintu, &
1300	& zintv &
1301	& )
1302	! At the end of the day also get interpolated means
1303	IF ( idayend == 0 ) THEN
1304	DEALLOCATE( &
1305	& zinmu, &
1306	& zinmv &
1307	& )
1308	ENDIF
1309
1310	prodatqc%nprofup = prodatqc%nprofup + ipro
1311
1312	END SUBROUTINE obs_vel_opt
1313
1314	END MODULE obs_oper
1315

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

Download in other formats: