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/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 10246

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

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

Download in other formats: