Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

obs_oper.F90

Last change on this file was 7713, checked in by dford, 7 years ago
Add observation operator code for surface log10(chlorophyll), SPM, pCO2 and fCO2, for use with FABM-ERSEM, HadOCC and MEDUSA. See internal Met Office NEMO ticket 660.
File size: 99.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_pro_sco_opt: Compute the model counterpart of temperature and
12	!! salinity observations from profiles in generalised
13	!! vertical coordinates
14	!! obs_sla_opt : Compute the model counterpart of sea level anomaly
15	!! observations
16	!! obs_sst_opt : Compute the model counterpart of sea surface temperature
17	!! observations
18	!! obs_sss_opt : Compute the model counterpart of sea surface salinity
19	!! observations
20	!! obs_seaice_opt : Compute the model counterpart of sea ice concentration
21	!! observations
22	!!
23	!! obs_vel_opt : Compute the model counterpart of zonal and meridional
24	!! components of velocity from observations.
25	!! obs_logchl_opt : Compute the model counterpart of log10(chlorophyll)
26	!! observations
27	!! obs_spm_opt : Compute the model counterpart of spm
28	!! observations
29	!! obs_fco2_opt : Compute the model counterpart of fco2
30	!! observations
31	!! obs_pco2_opt : Compute the model counterpart of pco2
32	!! observations
33	!!----------------------------------------------------------------------
34
35	!! * Modules used
36	USE par_kind, ONLY : & ! Precision variables
37	& wp
38	USE in_out_manager ! I/O manager
39	USE obs_inter_sup ! Interpolation support
40	USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the observation pt
41	& obs_int_h2d, &
42	& obs_int_h2d_init
43	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt
44	& obs_int_z1d, &
45	& obs_int_z1d_spl
46	USE obs_const, ONLY : &
47	& obfillflt ! Fillvalue
48	USE dom_oce, ONLY : &
49	& glamt, glamu, glamv, &
50	& gphit, gphiu, gphiv, &
51	#if defined key_vvl
52	& gdept_n
53	#else
54	& gdept_0
55	#endif
56	USE lib_mpp, ONLY : &
57	& ctl_warn, ctl_stop
58	USE obs_grid, ONLY : &
59	& obs_level_search
60
61	IMPLICIT NONE
62
63	!! * Routine accessibility
64	PRIVATE
65
66	PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations
67	& obs_pro_sco_opt, & ! Compute the model counterpart of profile observations
68	! in generalised vertical coordinates
69	& obs_sla_opt, & ! Compute the model counterpart of SLA observations
70	& obs_sst_opt, & ! Compute the model counterpart of SST observations
71	& obs_sss_opt, & ! Compute the model counterpart of SSS observations
72	& obs_seaice_opt, &
73	& obs_vel_opt, & ! Compute the model counterpart of velocity profile data
74	& obs_logchl_opt, & ! Compute the model counterpart of logchl data
75	& obs_spm_opt, & ! Compute the model counterpart of spm data
76	& obs_fco2_opt, & ! Compute the model counterpart of fco2 data
77	& obs_pco2_opt ! Compute the model counterpart of pco2 data
78
79	INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
80
81	!!----------------------------------------------------------------------
82	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
83	!! $Id$
84	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
85	!!----------------------------------------------------------------------
86
87	!! * Substitutions
88	# include "domzgr_substitute.h90"
89	CONTAINS
90
91	SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
92	& ptn, psn, pgdept, ptmask, k1dint, k2dint, &
93	& kdailyavtypes )
94	!!-----------------------------------------------------------------------
95	!!
96	!! * ROUTINE obs_pro_opt *
97	!!
98	!! ** Purpose : Compute the model counterpart of profiles
99	!! data by interpolating from the model grid to the
100	!! observation point.
101	!!
102	!! ** Method : Linearly interpolate to each observation point using
103	!! the model values at the corners of the surrounding grid box.
104	!!
105	!! First, a vertical profile of horizontally interpolated model
106	!! now temperatures is computed at the obs (lon, lat) point.
107	!! Several horizontal interpolation schemes are available:
108	!! - distance-weighted (great circle) (k2dint = 0)
109	!! - distance-weighted (small angle) (k2dint = 1)
110	!! - bilinear (geographical grid) (k2dint = 2)
111	!! - bilinear (quadrilateral grid) (k2dint = 3)
112	!! - polynomial (quadrilateral grid) (k2dint = 4)
113	!!
114	!! Next, the vertical temperature profile is interpolated to the
115	!! data depth points. Two vertical interpolation schemes are
116	!! available:
117	!! - linear (k1dint = 0)
118	!! - Cubic spline (k1dint = 1)
119	!!
120	!! For the cubic spline the 2nd derivative of the interpolating
121	!! polynomial is computed before entering the vertical interpolation
122	!! routine.
123	!!
124	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
125	!! a daily mean model temperature field. So, we first compute
126	!! the mean, then interpolate only at the end of the day.
127	!!
128	!! Note: the in situ temperature observations must be converted
129	!! to potential temperature (the model variable) prior to
130	!! assimilation.
131	!!??????????????????????????????????????????????????????????????
132	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
133	!!??????????????????????????????????????????????????????????????
134	!!
135	!! ** Action :
136	!!
137	!! History :
138	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
139	!! ! 06-03 (G. Smith) NEMOVAR migration
140	!! ! 06-10 (A. Weaver) Cleanup
141	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
142	!! ! 07-03 (K. Mogensen) General handling of profiles
143	!!-----------------------------------------------------------------------
144
145	!! * Modules used
146	USE obs_profiles_def ! Definition of storage space for profile obs.
147
148	IMPLICIT NONE
149
150	!! * Arguments
151	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
152	INTEGER, INTENT(IN) :: kt ! Time step
153	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
154	INTEGER, INTENT(IN) :: kpj
155	INTEGER, INTENT(IN) :: kpk
156	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
157	! (kit000-1 = restart time)
158	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
159	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
160	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
161	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
162	& ptn, & ! Model temperature field
163	& psn, & ! Model salinity field
164	& ptmask ! Land-sea mask
165	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
166	& pgdept ! Model array of depth levels
167	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
168	& kdailyavtypes! Types for daily averages
169	!! * Local declarations
170	INTEGER :: ji
171	INTEGER :: jj
172	INTEGER :: jk
173	INTEGER :: jobs
174	INTEGER :: inrc
175	INTEGER :: ipro
176	INTEGER :: idayend
177	INTEGER :: ista
178	INTEGER :: iend
179	INTEGER :: iobs
180	INTEGER, DIMENSION(imaxavtypes) :: &
181	& idailyavtypes
182	REAL(KIND=wp) :: zlam
183	REAL(KIND=wp) :: zphi
184	REAL(KIND=wp) :: zdaystp
185	REAL(KIND=wp), DIMENSION(kpk) :: &
186	& zobsmask, &
187	& zobsk, &
188	& zobs2k
189	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
190	& zweig
191	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
192	& zmask, &
193	& zintt, &
194	& zints, &
195	& zinmt, &
196	& zinms
197	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
198	& zglam, &
199	& zgphi
200	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
201	& igrdi, &
202	& igrdj
203
204	!------------------------------------------------------------------------
205	! Local initialization
206	!------------------------------------------------------------------------
207	! ... Record and data counters
208	inrc = kt - kit000 + 2
209	ipro = prodatqc%npstp(inrc)
210
211	! Daily average types
212	IF ( PRESENT(kdailyavtypes) ) THEN
213	idailyavtypes(:) = kdailyavtypes(:)
214	ELSE
215	idailyavtypes(:) = -1
216	ENDIF
217
218	! Initialize daily mean for first timestep
219	idayend = MOD( kt - kit000 + 1, kdaystp )
220
221	! Added kt == 0 test to catch restart case
222	IF ( idayend == 1 .OR. kt == 0) THEN
223	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
224	DO jk = 1, jpk
225	DO jj = 1, jpj
226	DO ji = 1, jpi
227	prodatqc%vdmean(ji,jj,jk,1) = 0.0
228	prodatqc%vdmean(ji,jj,jk,2) = 0.0
229	END DO
230	END DO
231	END DO
232	ENDIF
233
234	DO jk = 1, jpk
235	DO jj = 1, jpj
236	DO ji = 1, jpi
237	! Increment the temperature field for computing daily mean
238	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
239	& + ptn(ji,jj,jk)
240	! Increment the salinity field for computing daily mean
241	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
242	& + psn(ji,jj,jk)
243	END DO
244	END DO
245	END DO
246
247	! Compute the daily mean at the end of day
248	zdaystp = 1.0 / REAL( kdaystp )
249	IF ( idayend == 0 ) THEN
250	DO jk = 1, jpk
251	DO jj = 1, jpj
252	DO ji = 1, jpi
253	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
254	& * zdaystp
255	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
256	& * zdaystp
257	END DO
258	END DO
259	END DO
260	ENDIF
261
262	! Get the data for interpolation
263	ALLOCATE( &
264	& igrdi(2,2,ipro), &
265	& igrdj(2,2,ipro), &
266	& zglam(2,2,ipro), &
267	& zgphi(2,2,ipro), &
268	& zmask(2,2,kpk,ipro), &
269	& zintt(2,2,kpk,ipro), &
270	& zints(2,2,kpk,ipro) &
271	& )
272
273	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
274	iobs = jobs - prodatqc%nprofup
275	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
276	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
277	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
278	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
279	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
280	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
281	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
282	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
283	END DO
284
285	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
286	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
287	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
288	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
289	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
290
291	! At the end of the day also get interpolated means
292	IF ( idayend == 0 ) THEN
293
294	ALLOCATE( &
295	& zinmt(2,2,kpk,ipro), &
296	& zinms(2,2,kpk,ipro) &
297	& )
298
299	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
300	& prodatqc%vdmean(:,:,:,1), zinmt )
301	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
302	& prodatqc%vdmean(:,:,:,2), zinms )
303
304	ENDIF
305
306	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
307
308	iobs = jobs - prodatqc%nprofup
309
310	IF ( kt /= prodatqc%mstp(jobs) ) THEN
311
312	IF(lwp) THEN
313	WRITE(numout,*)
314	WRITE(numout,*) ' E R R O R : Observation', &
315	& ' time step is not consistent with the', &
316	& ' model time step'
317	WRITE(numout,*) ' ========='
318	WRITE(numout,*)
319	WRITE(numout,*) ' Record = ', jobs, &
320	& ' kt = ', kt, &
321	& ' mstp = ', prodatqc%mstp(jobs), &
322	& ' ntyp = ', prodatqc%ntyp(jobs)
323	ENDIF
324	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
325	ENDIF
326
327	zlam = prodatqc%rlam(jobs)
328	zphi = prodatqc%rphi(jobs)
329
330	! Horizontal weights and vertical mask
331
332	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
333	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
334
335	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
336	& zglam(:,:,iobs), zgphi(:,:,iobs), &
337	& zmask(:,:,:,iobs), zweig, zobsmask )
338
339	ENDIF
340
341	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
342
343	zobsk(:) = obfillflt
344
345	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
346
347	IF ( idayend == 0 ) THEN
348
349	! Daily averaged moored buoy (MRB) data
350
351	CALL obs_int_h2d( kpk, kpk, &
352	& zweig, zinmt(:,:,:,iobs), zobsk )
353
354
355	ELSE
356
357	CALL ctl_stop( ' A nonzero' // &
358	& ' number of profile T BUOY data should' // &
359	& ' only occur at the end of a given day' )
360
361	ENDIF
362
363	ELSE
364
365	! Point data
366
367	CALL obs_int_h2d( kpk, kpk, &
368	& zweig, zintt(:,:,:,iobs), zobsk )
369
370	ENDIF
371
372	!-------------------------------------------------------------
373	! Compute vertical second-derivative of the interpolating
374	! polynomial at obs points
375	!-------------------------------------------------------------
376
377	IF ( k1dint == 1 ) THEN
378	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
379	& pgdept, zobsmask )
380	ENDIF
381
382	!-----------------------------------------------------------------
383	! Vertical interpolation to the observation point
384	!-----------------------------------------------------------------
385	ista = prodatqc%npvsta(jobs,1)
386	iend = prodatqc%npvend(jobs,1)
387	CALL obs_int_z1d( kpk, &
388	& prodatqc%var(1)%mvk(ista:iend), &
389	& k1dint, iend - ista + 1, &
390	& prodatqc%var(1)%vdep(ista:iend), &
391	& zobsk, zobs2k, &
392	& prodatqc%var(1)%vmod(ista:iend), &
393	& pgdept, zobsmask )
394
395	ENDIF
396
397	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
398
399	zobsk(:) = obfillflt
400
401	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
402
403	IF ( idayend == 0 ) THEN
404
405	! Daily averaged moored buoy (MRB) data
406
407	CALL obs_int_h2d( kpk, kpk, &
408	& zweig, zinms(:,:,:,iobs), zobsk )
409
410	ELSE
411
412	CALL ctl_stop( ' A nonzero' // &
413	& ' number of profile S BUOY data should' // &
414	& ' only occur at the end of a given day' )
415
416	ENDIF
417
418	ELSE
419
420	! Point data
421
422	CALL obs_int_h2d( kpk, kpk, &
423	& zweig, zints(:,:,:,iobs), zobsk )
424
425	ENDIF
426
427
428	!-------------------------------------------------------------
429	! Compute vertical second-derivative of the interpolating
430	! polynomial at obs points
431	!-------------------------------------------------------------
432
433	IF ( k1dint == 1 ) THEN
434	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
435	& pgdept, zobsmask )
436	ENDIF
437
438	!----------------------------------------------------------------
439	! Vertical interpolation to the observation point
440	!----------------------------------------------------------------
441	ista = prodatqc%npvsta(jobs,2)
442	iend = prodatqc%npvend(jobs,2)
443	CALL obs_int_z1d( kpk, &
444	& prodatqc%var(2)%mvk(ista:iend),&
445	& k1dint, iend - ista + 1, &
446	& prodatqc%var(2)%vdep(ista:iend),&
447	& zobsk, zobs2k, &
448	& prodatqc%var(2)%vmod(ista:iend),&
449	& pgdept, zobsmask )
450
451	ENDIF
452
453	END DO
454
455	! Deallocate the data for interpolation
456	DEALLOCATE( &
457	& igrdi, &
458	& igrdj, &
459	& zglam, &
460	& zgphi, &
461	& zmask, &
462	& zintt, &
463	& zints &
464	& )
465	! At the end of the day also get interpolated means
466	IF ( idayend == 0 ) THEN
467	DEALLOCATE( &
468	& zinmt, &
469	& zinms &
470	& )
471	ENDIF
472
473	prodatqc%nprofup = prodatqc%nprofup + ipro
474
475	END SUBROUTINE obs_pro_opt
476
477	SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
478	& ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, &
479	& kdailyavtypes )
480	!!-----------------------------------------------------------------------
481	!!
482	!! * ROUTINE obs_pro_opt *
483	!!
484	!! ** Purpose : Compute the model counterpart of profiles
485	!! data by interpolating from the model grid to the
486	!! observation point. Generalised vertical coordinate version
487	!!
488	!! ** Method : Linearly interpolate to each observation point using
489	!! the model values at the corners of the surrounding grid box.
490	!!
491	!! First, model values on the model grid are interpolated vertically to the
492	!! Depths of the profile observations. Two vertical interpolation schemes are
493	!! available:
494	!! - linear (k1dint = 0)
495	!! - Cubic spline (k1dint = 1)
496	!!
497	!!
498	!! Secondly the interpolated values are interpolated horizontally to the
499	!! obs (lon, lat) point.
500	!! Several horizontal interpolation schemes are available:
501	!! - distance-weighted (great circle) (k2dint = 0)
502	!! - distance-weighted (small angle) (k2dint = 1)
503	!! - bilinear (geographical grid) (k2dint = 2)
504	!! - bilinear (quadrilateral grid) (k2dint = 3)
505	!! - polynomial (quadrilateral grid) (k2dint = 4)
506	!!
507	!! For the cubic spline the 2nd derivative of the interpolating
508	!! polynomial is computed before entering the vertical interpolation
509	!! routine.
510	!!
511	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
512	!! a daily mean model temperature field. So, we first compute
513	!! the mean, then interpolate only at the end of the day.
514	!!
515	!! This is the procedure to be used with generalised vertical model
516	!! coordinates (ie s-coordinates. It is ~4x slower than the equivalent
517	!! horizontal then vertical interpolation algorithm, but can deal with situations
518	!! where the model levels are not flat.
519	!! ONLY PERFORMED if ln_sco=.TRUE.
520	!!
521	!! Note: the in situ temperature observations must be converted
522	!! to potential temperature (the model variable) prior to
523	!! assimilation.
524	!!??????????????????????????????????????????????????????????????
525	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
526	!!??????????????????????????????????????????????????????????????
527	!!
528	!! ** Action :
529	!!
530	!! History :
531	!! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised
532	!! vertical coordinates
533	!!-----------------------------------------------------------------------
534
535	!! * Modules used
536	USE obs_profiles_def ! Definition of storage space for profile obs.
537	USE dom_oce, ONLY : &
538	#if defined key_vvl
539	& gdepw_n
540	#else
541	& gdepw_0
542	#endif
543
544	IMPLICIT NONE
545
546	!! * Arguments
547	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
548	INTEGER, INTENT(IN) :: kt ! Time step
549	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
550	INTEGER, INTENT(IN) :: kpj
551	INTEGER, INTENT(IN) :: kpk
552	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
553	! (kit000-1 = restart time)
554	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
555	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
556	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
557	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
558	& ptn, & ! Model temperature field
559	& psn, & ! Model salinity field
560	& ptmask ! Land-sea mask
561	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
562	& pgdept, & ! Model array of depth T levels
563	& pgdepw ! Model array of depth W levels
564	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
565	& kdailyavtypes ! Types for daily averages
566
567	!! * Local declarations
568	INTEGER :: ji
569	INTEGER :: jj
570	INTEGER :: jk
571	INTEGER :: iico, ijco
572	INTEGER :: jobs
573	INTEGER :: inrc
574	INTEGER :: ipro
575	INTEGER :: idayend
576	INTEGER :: ista
577	INTEGER :: iend
578	INTEGER :: iobs
579	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
580	INTEGER, DIMENSION(imaxavtypes) :: &
581	& idailyavtypes
582	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
583	& igrdi, &
584	& igrdj
585	INTEGER :: &
586	& inum_obs
587	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
588	REAL(KIND=wp) :: zlam
589	REAL(KIND=wp) :: zphi
590	REAL(KIND=wp) :: zdaystp
591	REAL(KIND=wp), DIMENSION(kpk) :: &
592	& zobsmask, &
593	& zobsk, &
594	& zobs2k
595	REAL(KIND=wp), DIMENSION(2,2,1) :: &
596	& zweig, &
597	& l_zweig
598	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
599	& zmask, &
600	& zintt, &
601	& zints, &
602	& zinmt, &
603	& zgdept,&
604	& zgdepw,&
605	& zinms
606	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
607	& zglam, &
608	& zgphi
609	REAL(KIND=wp), DIMENSION(1) :: zmsk_1
610	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
611
612	!------------------------------------------------------------------------
613	! Local initialization
614	!------------------------------------------------------------------------
615	! ... Record and data counters
616	inrc = kt - kit000 + 2
617	ipro = prodatqc%npstp(inrc)
618
619	! Daily average types
620	IF ( PRESENT(kdailyavtypes) ) THEN
621	idailyavtypes(:) = kdailyavtypes(:)
622	ELSE
623	idailyavtypes(:) = -1
624	ENDIF
625
626	! Initialize daily mean for first time-step
627	idayend = MOD( kt - kit000 + 1, kdaystp )
628
629	! Added kt == 0 test to catch restart case
630	IF ( idayend == 1 .OR. kt == 0) THEN
631
632	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
633	DO jk = 1, jpk
634	DO jj = 1, jpj
635	DO ji = 1, jpi
636	prodatqc%vdmean(ji,jj,jk,1) = 0.0
637	prodatqc%vdmean(ji,jj,jk,2) = 0.0
638	END DO
639	END DO
640	END DO
641
642	ENDIF
643
644	DO jk = 1, jpk
645	DO jj = 1, jpj
646	DO ji = 1, jpi
647	! Increment the temperature field for computing daily mean
648	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
649	& + ptn(ji,jj,jk)
650	! Increment the salinity field for computing daily mean
651	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
652	& + psn(ji,jj,jk)
653	END DO
654	END DO
655	END DO
656
657	! Compute the daily mean at the end of day
658	zdaystp = 1.0 / REAL( kdaystp )
659	IF ( idayend == 0 ) THEN
660	DO jk = 1, jpk
661	DO jj = 1, jpj
662	DO ji = 1, jpi
663	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
664	& * zdaystp
665	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
666	& * zdaystp
667	END DO
668	END DO
669	END DO
670	ENDIF
671
672	! Get the data for interpolation
673	ALLOCATE( &
674	& igrdi(2,2,ipro), &
675	& igrdj(2,2,ipro), &
676	& zglam(2,2,ipro), &
677	& zgphi(2,2,ipro), &
678	& zmask(2,2,kpk,ipro), &
679	& zintt(2,2,kpk,ipro), &
680	& zints(2,2,kpk,ipro), &
681	& zgdept(2,2,kpk,ipro), &
682	& zgdepw(2,2,kpk,ipro) &
683	& )
684
685	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
686	iobs = jobs - prodatqc%nprofup
687	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
688	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
689	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
690	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
691	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
692	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
693	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
694	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
695	END DO
696
697	! Initiialise depth arrays
698	zgdept = 0.0
699	zgdepw = 0.0
700
701	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
702	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
703	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
704	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
705	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
706	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdept(:,:,:), &
707	& zgdept )
708	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdepw(:,:,:), &
709	& zgdepw )
710
711	! At the end of the day also get interpolated means
712	IF ( idayend == 0 ) THEN
713
714	ALLOCATE( &
715	& zinmt(2,2,kpk,ipro), &
716	& zinms(2,2,kpk,ipro) &
717	& )
718
719	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
720	& prodatqc%vdmean(:,:,:,1), zinmt )
721	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
722	& prodatqc%vdmean(:,:,:,2), zinms )
723
724	ENDIF
725
726	! Return if no observations to process
727	! Has to be done after comm commands to ensure processors
728	! stay in sync
729	IF ( ipro == 0 ) RETURN
730
731	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
732
733	iobs = jobs - prodatqc%nprofup
734
735	IF ( kt /= prodatqc%mstp(jobs) ) THEN
736
737	IF(lwp) THEN
738	WRITE(numout,*)
739	WRITE(numout,*) ' E R R O R : Observation', &
740	& ' time step is not consistent with the', &
741	& ' model time step'
742	WRITE(numout,*) ' ========='
743	WRITE(numout,*)
744	WRITE(numout,*) ' Record = ', jobs, &
745	& ' kt = ', kt, &
746	& ' mstp = ', prodatqc%mstp(jobs), &
747	& ' ntyp = ', prodatqc%ntyp(jobs)
748	ENDIF
749	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
750	ENDIF
751
752	zlam = prodatqc%rlam(jobs)
753	zphi = prodatqc%rphi(jobs)
754
755	! Horizontal weights
756	! Only calculated once, for both T and S.
757	! Masked values are calculated later.
758
759	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
760	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
761
762	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
763	& zglam(:,:,iobs), zgphi(:,:,iobs), &
764	& zmask(:,:,1,iobs), zweig, zmsk_1 )
765
766	ENDIF
767
768	! IF zmsk_1 = 0; then ob is on land
769	IF (zmsk_1(1) < 0.1) THEN
770	WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'
771
772	ELSE
773
774	! Temperature
775
776	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
777
778	zobsk(:) = obfillflt
779
780	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
781
782	IF ( idayend == 0 ) THEN
783
784	! Daily averaged moored buoy (MRB) data
785
786	! vertically interpolate all 4 corners
787	ista = prodatqc%npvsta(jobs,1)
788	iend = prodatqc%npvend(jobs,1)
789	inum_obs = iend - ista + 1
790	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
791
792	DO iin=1,2
793	DO ijn=1,2
794
795
796
797	IF ( k1dint == 1 ) THEN
798	CALL obs_int_z1d_spl( kpk, &
799	& zinmt(iin,ijn,:,iobs), &
800	& zobs2k, zgdept(iin,ijn,:,iobs), &
801	& zmask(iin,ijn,:,iobs))
802	ENDIF
803
804	CALL obs_level_search(kpk, &
805	& zgdept(iin,ijn,:,iobs), &
806	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
807	& iv_indic)
808	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
809	& prodatqc%var(1)%vdep(ista:iend), &
810	& zinmt(iin,ijn,:,iobs), &
811	& zobs2k, interp_corner(iin,ijn,:), &
812	& zgdept(iin,ijn,:,iobs), &
813	& zmask(iin,ijn,:,iobs))
814
815	ENDDO
816	ENDDO
817
818
819	ELSE
820
821	CALL ctl_stop( ' A nonzero' // &
822	& ' number of profile T BUOY data should' // &
823	& ' only occur at the end of a given day' )
824
825	ENDIF
826
827	ELSE
828
829	! Point data
830
831	! vertically interpolate all 4 corners
832	ista = prodatqc%npvsta(jobs,1)
833	iend = prodatqc%npvend(jobs,1)
834	inum_obs = iend - ista + 1
835	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
836	DO iin=1,2
837	DO ijn=1,2
838
839
840	IF ( k1dint == 1 ) THEN
841	CALL obs_int_z1d_spl( kpk, &
842	& zintt(iin,ijn,:,iobs),&
843	& zobs2k, zgdept(iin,ijn,:,iobs), &
844	& zmask(iin,ijn,:,iobs))
845
846	ENDIF
847
848	CALL obs_level_search(kpk, &
849	& zgdept(iin,ijn,:,iobs),&
850	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
851	& iv_indic)
852	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
853	& prodatqc%var(1)%vdep(ista:iend), &
854	& zintt(iin,ijn,:,iobs), &
855	& zobs2k,interp_corner(iin,ijn,:), &
856	& zgdept(iin,ijn,:,iobs), &
857	& zmask(iin,ijn,:,iobs) )
858
859	ENDDO
860	ENDDO
861
862	ENDIF
863
864	!-------------------------------------------------------------
865	! Compute the horizontal interpolation for every profile level
866	!-------------------------------------------------------------
867
868	DO ikn=1,inum_obs
869	iend=ista+ikn-1
870
871	l_zweig(:,:,1) = 0._wp
872
873	! This code forces the horizontal weights to be
874	! zero IF the observation is below the bottom of the
875	! corners of the interpolation nodes, Or if it is in
876	! the mask. This is important for observations are near
877	! steep bathymetry
878	DO iin=1,2
879	DO ijn=1,2
880
881	depth_loop1: DO ik=kpk,2,-1
882	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
883
884	l_zweig(iin,ijn,1) = &
885	& zweig(iin,ijn,1) * &
886	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
887	& - prodatqc%var(1)%vdep(iend)),0._wp)
888
889	EXIT depth_loop1
890	ENDIF
891	ENDDO depth_loop1
892
893	ENDDO
894	ENDDO
895
896	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
897	& prodatqc%var(1)%vmod(iend:iend) )
898
899	! Set QC flag for any observations found below the bottom
900	! needed as the check here is more strict than that in obs_prep
901	IF (sum(l_zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
902
903	ENDDO
904
905
906	DEALLOCATE(interp_corner,iv_indic)
907
908	ENDIF
909
910
911	! Salinity
912
913	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
914
915	zobsk(:) = obfillflt
916
917	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
918
919	IF ( idayend == 0 ) THEN
920
921	! Daily averaged moored buoy (MRB) data
922
923	! vertically interpolate all 4 corners
924	ista = prodatqc%npvsta(jobs,2)
925	iend = prodatqc%npvend(jobs,2)
926	inum_obs = iend - ista + 1
927	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
928
929	DO iin=1,2
930	DO ijn=1,2
931
932
933
934	IF ( k1dint == 1 ) THEN
935	CALL obs_int_z1d_spl( kpk, &
936	& zinms(iin,ijn,:,iobs), &
937	& zobs2k, zgdept(iin,ijn,:,iobs), &
938	& zmask(iin,ijn,:,iobs))
939	ENDIF
940
941	CALL obs_level_search(kpk, &
942	& zgdept(iin,ijn,:,iobs), &
943	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
944	& iv_indic)
945	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
946	& prodatqc%var(2)%vdep(ista:iend), &
947	& zinms(iin,ijn,:,iobs), &
948	& zobs2k, interp_corner(iin,ijn,:), &
949	& zgdept(iin,ijn,:,iobs), &
950	& zmask(iin,ijn,:,iobs))
951
952	ENDDO
953	ENDDO
954
955
956	ELSE
957
958	CALL ctl_stop( ' A nonzero' // &
959	& ' number of profile T BUOY data should' // &
960	& ' only occur at the end of a given day' )
961
962	ENDIF
963
964	ELSE
965
966	! Point data
967
968	! vertically interpolate all 4 corners
969	ista = prodatqc%npvsta(jobs,2)
970	iend = prodatqc%npvend(jobs,2)
971	inum_obs = iend - ista + 1
972	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
973
974	DO iin=1,2
975	DO ijn=1,2
976
977
978	IF ( k1dint == 1 ) THEN
979	CALL obs_int_z1d_spl( kpk, &
980	& zints(iin,ijn,:,iobs),&
981	& zobs2k, zgdept(iin,ijn,:,iobs), &
982	& zmask(iin,ijn,:,iobs))
983
984	ENDIF
985
986	CALL obs_level_search(kpk, &
987	& zgdept(iin,ijn,:,iobs),&
988	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
989	& iv_indic)
990	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
991	& prodatqc%var(2)%vdep(ista:iend), &
992	& zints(iin,ijn,:,iobs), &
993	& zobs2k,interp_corner(iin,ijn,:), &
994	& zgdept(iin,ijn,:,iobs), &
995	& zmask(iin,ijn,:,iobs) )
996
997	ENDDO
998	ENDDO
999
1000	ENDIF
1001
1002	!-------------------------------------------------------------
1003	! Compute the horizontal interpolation for every profile level
1004	!-------------------------------------------------------------
1005
1006	DO ikn=1,inum_obs
1007	iend=ista+ikn-1
1008
1009	l_zweig(:,:,1) = 0._wp
1010
1011	! This code forces the horizontal weights to be
1012	! zero IF the observation is below the bottom of the
1013	! corners of the interpolation nodes, Or if it is in
1014	! the mask. This is important for observations are near
1015	! steep bathymetry
1016	DO iin=1,2
1017	DO ijn=1,2
1018
1019	depth_loop2: DO ik=kpk,2,-1
1020	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
1021
1022	l_zweig(iin,ijn,1) = &
1023	& zweig(iin,ijn,1) * &
1024	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
1025	& - prodatqc%var(2)%vdep(iend)),0._wp)
1026
1027	EXIT depth_loop2
1028	ENDIF
1029	ENDDO depth_loop2
1030
1031	ENDDO
1032	ENDDO
1033
1034	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
1035	& prodatqc%var(2)%vmod(iend:iend) )
1036
1037	! Set QC flag for any observations found below the bottom
1038	! needed as the check here is more strict than that in obs_prep
1039	IF (sum(l_zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
1040
1041	ENDDO
1042
1043
1044	DEALLOCATE(interp_corner,iv_indic)
1045
1046	ENDIF
1047
1048	ENDIF
1049
1050	END DO
1051
1052	! Deallocate the data for interpolation
1053	DEALLOCATE( &
1054	& igrdi, &
1055	& igrdj, &
1056	& zglam, &
1057	& zgphi, &
1058	& zmask, &
1059	& zintt, &
1060	& zints, &
1061	& zgdept,&
1062	& zgdepw &
1063	& )
1064	! At the end of the day also get interpolated means
1065	IF ( idayend == 0 ) THEN
1066	DEALLOCATE( &
1067	& zinmt, &
1068	& zinms &
1069	& )
1070	ENDIF
1071
1072	prodatqc%nprofup = prodatqc%nprofup + ipro
1073
1074	END SUBROUTINE obs_pro_sco_opt
1075
1076	SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
1077	& psshn, psshmask, k2dint )
1078	!!-----------------------------------------------------------------------
1079	!!
1080	!! * ROUTINE obs_sla_opt *
1081	!!
1082	!! ** Purpose : Compute the model counterpart of sea level anomaly
1083	!! data by interpolating from the model grid to the
1084	!! observation point.
1085	!!
1086	!! ** Method : Linearly interpolate to each observation point using
1087	!! the model values at the corners of the surrounding grid box.
1088	!!
1089	!! The now model SSH is first computed at the obs (lon, lat) point.
1090	!!
1091	!! Several horizontal interpolation schemes are available:
1092	!! - distance-weighted (great circle) (k2dint = 0)
1093	!! - distance-weighted (small angle) (k2dint = 1)
1094	!! - bilinear (geographical grid) (k2dint = 2)
1095	!! - bilinear (quadrilateral grid) (k2dint = 3)
1096	!! - polynomial (quadrilateral grid) (k2dint = 4)
1097	!!
1098	!! The sea level anomaly at the observation points is then computed
1099	!! by removing a mean dynamic topography (defined at the obs. point).
1100	!!
1101	!! ** Action :
1102	!!
1103	!! History :
1104	!! ! 07-03 (A. Weaver)
1105	!!-----------------------------------------------------------------------
1106
1107	!! * Modules used
1108	USE obs_surf_def ! Definition of storage space for surface observations
1109
1110	IMPLICIT NONE
1111
1112	!! * Arguments
1113	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening
1114	INTEGER, INTENT(IN) :: kt ! Time step
1115	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1116	INTEGER, INTENT(IN) :: kpj
1117	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1118	! (kit000-1 = restart time)
1119	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1120	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1121	& psshn, & ! Model SSH field
1122	& psshmask ! Land-sea mask
1123
1124	!! * Local declarations
1125	INTEGER :: ji
1126	INTEGER :: jj
1127	INTEGER :: jobs
1128	INTEGER :: inrc
1129	INTEGER :: isla
1130	INTEGER :: iobs
1131	REAL(KIND=wp) :: zlam
1132	REAL(KIND=wp) :: zphi
1133	REAL(KIND=wp) :: zext(1), zobsmask(1)
1134	REAL(kind=wp), DIMENSION(2,2,1) :: &
1135	& zweig
1136	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1137	& zmask, &
1138	& zsshl, &
1139	& zglam, &
1140	& zgphi
1141	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1142	& igrdi, &
1143	& igrdj
1144
1145	!------------------------------------------------------------------------
1146	! Local initialization
1147	!------------------------------------------------------------------------
1148	! ... Record and data counters
1149	inrc = kt - kit000 + 2
1150	isla = sladatqc%nsstp(inrc)
1151
1152	! Get the data for interpolation
1153
1154	ALLOCATE( &
1155	& igrdi(2,2,isla), &
1156	& igrdj(2,2,isla), &
1157	& zglam(2,2,isla), &
1158	& zgphi(2,2,isla), &
1159	& zmask(2,2,isla), &
1160	& zsshl(2,2,isla) &
1161	& )
1162
1163	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1164	iobs = jobs - sladatqc%nsurfup
1165	igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
1166	igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
1167	igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
1168	igrdj(1,2,iobs) = sladatqc%mj(jobs)
1169	igrdi(2,1,iobs) = sladatqc%mi(jobs)
1170	igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
1171	igrdi(2,2,iobs) = sladatqc%mi(jobs)
1172	igrdj(2,2,iobs) = sladatqc%mj(jobs)
1173	END DO
1174
1175	CALL obs_int_comm_2d( 2, 2, isla, &
1176	& igrdi, igrdj, glamt, zglam )
1177	CALL obs_int_comm_2d( 2, 2, isla, &
1178	& igrdi, igrdj, gphit, zgphi )
1179	CALL obs_int_comm_2d( 2, 2, isla, &
1180	& igrdi, igrdj, psshmask, zmask )
1181	CALL obs_int_comm_2d( 2, 2, isla, &
1182	& igrdi, igrdj, psshn, zsshl )
1183
1184	! Loop over observations
1185
1186	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1187
1188	iobs = jobs - sladatqc%nsurfup
1189
1190	IF ( kt /= sladatqc%mstp(jobs) ) THEN
1191
1192	IF(lwp) THEN
1193	WRITE(numout,*)
1194	WRITE(numout,*) ' E R R O R : Observation', &
1195	& ' time step is not consistent with the', &
1196	& ' model time step'
1197	WRITE(numout,*) ' ========='
1198	WRITE(numout,*)
1199	WRITE(numout,*) ' Record = ', jobs, &
1200	& ' kt = ', kt, &
1201	& ' mstp = ', sladatqc%mstp(jobs), &
1202	& ' ntyp = ', sladatqc%ntyp(jobs)
1203	ENDIF
1204	CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
1205
1206	ENDIF
1207
1208	zlam = sladatqc%rlam(jobs)
1209	zphi = sladatqc%rphi(jobs)
1210
1211	! Get weights to interpolate the model SSH to the observation point
1212	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1213	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1214	& zmask(:,:,iobs), zweig, zobsmask )
1215
1216
1217	! Interpolate the model SSH to the observation point
1218	CALL obs_int_h2d( 1, 1, &
1219	& zweig, zsshl(:,:,iobs), zext )
1220
1221	sladatqc%rext(jobs,1) = zext(1)
1222	! ... Remove the MDT at the observation point
1223	sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
1224
1225	END DO
1226
1227	! Deallocate the data for interpolation
1228	DEALLOCATE( &
1229	& igrdi, &
1230	& igrdj, &
1231	& zglam, &
1232	& zgphi, &
1233	& zmask, &
1234	& zsshl &
1235	& )
1236
1237	sladatqc%nsurfup = sladatqc%nsurfup + isla
1238
1239	END SUBROUTINE obs_sla_opt
1240
1241	SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
1242	& psstn, psstmask, k2dint, ld_nightav )
1243	!!-----------------------------------------------------------------------
1244	!!
1245	!! * ROUTINE obs_sst_opt *
1246	!!
1247	!! ** Purpose : Compute the model counterpart of surface temperature
1248	!! data by interpolating from the model grid to the
1249	!! observation point.
1250	!!
1251	!! ** Method : Linearly interpolate to each observation point using
1252	!! the model values at the corners of the surrounding grid box.
1253	!!
1254	!! The now model SST is first computed at the obs (lon, lat) point.
1255	!!
1256	!! Several horizontal interpolation schemes are available:
1257	!! - distance-weighted (great circle) (k2dint = 0)
1258	!! - distance-weighted (small angle) (k2dint = 1)
1259	!! - bilinear (geographical grid) (k2dint = 2)
1260	!! - bilinear (quadrilateral grid) (k2dint = 3)
1261	!! - polynomial (quadrilateral grid) (k2dint = 4)
1262	!!
1263	!!
1264	!! ** Action :
1265	!!
1266	!! History :
1267	!! ! 07-07 (S. Ricci ) : Original
1268	!!
1269	!!-----------------------------------------------------------------------
1270
1271	!! * Modules used
1272	USE obs_surf_def ! Definition of storage space for surface observations
1273	USE sbcdcy
1274
1275	IMPLICIT NONE
1276
1277	!! * Arguments
1278	TYPE(obs_surf), INTENT(INOUT) :: &
1279	& sstdatqc ! Subset of surface data not failing screening
1280	INTEGER, INTENT(IN) :: kt ! Time step
1281	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1282	INTEGER, INTENT(IN) :: kpj
1283	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1284	! (kit000-1 = restart time)
1285	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1286	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1287	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1288	& psstn, & ! Model SST field
1289	& psstmask ! Land-sea mask
1290
1291	!! * Local declarations
1292	INTEGER :: ji
1293	INTEGER :: jj
1294	INTEGER :: jobs
1295	INTEGER :: inrc
1296	INTEGER :: isst
1297	INTEGER :: iobs
1298	INTEGER :: idayend
1299	REAL(KIND=wp) :: zlam
1300	REAL(KIND=wp) :: zphi
1301	REAL(KIND=wp) :: zext(1), zobsmask(1)
1302	REAL(KIND=wp) :: zdaystp
1303	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1304	& icount_sstnight, &
1305	& imask_night
1306	REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1307	& zintmp, &
1308	& zouttmp, &
1309	& zmeanday ! to compute model sst in region of 24h daylight (pole)
1310	REAL(kind=wp), DIMENSION(2,2,1) :: &
1311	& zweig
1312	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1313	& zmask, &
1314	& zsstl, &
1315	& zsstm, &
1316	& zglam, &
1317	& zgphi
1318	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1319	& igrdi, &
1320	& igrdj
1321	LOGICAL, INTENT(IN) :: ld_nightav
1322
1323	!-----------------------------------------------------------------------
1324	! Local initialization
1325	!-----------------------------------------------------------------------
1326	! ... Record and data counters
1327	inrc = kt - kit000 + 2
1328	isst = sstdatqc%nsstp(inrc)
1329
1330	IF ( ld_nightav ) THEN
1331
1332	! Initialize array for night mean
1333
1334	IF ( kt .EQ. 0 ) THEN
1335	ALLOCATE ( icount_sstnight(kpi,kpj) )
1336	ALLOCATE ( imask_night(kpi,kpj) )
1337	ALLOCATE ( zintmp(kpi,kpj) )
1338	ALLOCATE ( zouttmp(kpi,kpj) )
1339	ALLOCATE ( zmeanday(kpi,kpj) )
1340	nday_qsr = -1 ! initialisation flag for nbc_dcy
1341	ENDIF
1342
1343	! Initialize daily mean for first timestep
1344	idayend = MOD( kt - kit000 + 1, kdaystp )
1345
1346	! Added kt == 0 test to catch restart case
1347	IF ( idayend == 1 .OR. kt == 0) THEN
1348	IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
1349	DO jj = 1, jpj
1350	DO ji = 1, jpi
1351	sstdatqc%vdmean(ji,jj) = 0.0
1352	zmeanday(ji,jj) = 0.0
1353	icount_sstnight(ji,jj) = 0
1354	END DO
1355	END DO
1356	ENDIF
1357
1358	zintmp(:,:) = 0.0
1359	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1360	imask_night(:,:) = INT( zouttmp(:,:) )
1361
1362	DO jj = 1, jpj
1363	DO ji = 1, jpi
1364	! Increment the temperature field for computing night mean and counter
1365	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1366	& + psstn(ji,jj)*imask_night(ji,jj)
1367	zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj)
1368	icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
1369	END DO
1370	END DO
1371
1372	! Compute the daily mean at the end of day
1373
1374	zdaystp = 1.0 / REAL( kdaystp )
1375
1376	IF ( idayend == 0 ) THEN
1377	DO jj = 1, jpj
1378	DO ji = 1, jpi
1379	! Test if "no night" point
1380	IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
1381	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1382	& / icount_sstnight(ji,jj)
1383	ELSE
1384	sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1385	ENDIF
1386	END DO
1387	END DO
1388	ENDIF
1389
1390	ENDIF
1391
1392	! Get the data for interpolation
1393
1394	ALLOCATE( &
1395	& igrdi(2,2,isst), &
1396	& igrdj(2,2,isst), &
1397	& zglam(2,2,isst), &
1398	& zgphi(2,2,isst), &
1399	& zmask(2,2,isst), &
1400	& zsstl(2,2,isst) &
1401	& )
1402
1403	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1404	iobs = jobs - sstdatqc%nsurfup
1405	igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
1406	igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
1407	igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
1408	igrdj(1,2,iobs) = sstdatqc%mj(jobs)
1409	igrdi(2,1,iobs) = sstdatqc%mi(jobs)
1410	igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
1411	igrdi(2,2,iobs) = sstdatqc%mi(jobs)
1412	igrdj(2,2,iobs) = sstdatqc%mj(jobs)
1413	END DO
1414
1415	CALL obs_int_comm_2d( 2, 2, isst, &
1416	& igrdi, igrdj, glamt, zglam )
1417	CALL obs_int_comm_2d( 2, 2, isst, &
1418	& igrdi, igrdj, gphit, zgphi )
1419	CALL obs_int_comm_2d( 2, 2, isst, &
1420	& igrdi, igrdj, psstmask, zmask )
1421	CALL obs_int_comm_2d( 2, 2, isst, &
1422	& igrdi, igrdj, psstn, zsstl )
1423
1424	! At the end of the day get interpolated means
1425	IF ( idayend == 0 .AND. ld_nightav ) THEN
1426
1427	ALLOCATE( &
1428	& zsstm(2,2,isst) &
1429	& )
1430
1431	CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
1432	& sstdatqc%vdmean(:,:), zsstm )
1433
1434	ENDIF
1435
1436	! Loop over observations
1437
1438	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1439
1440	iobs = jobs - sstdatqc%nsurfup
1441
1442	IF ( kt /= sstdatqc%mstp(jobs) ) THEN
1443
1444	IF(lwp) THEN
1445	WRITE(numout,*)
1446	WRITE(numout,*) ' E R R O R : Observation', &
1447	& ' time step is not consistent with the', &
1448	& ' model time step'
1449	WRITE(numout,*) ' ========='
1450	WRITE(numout,*)
1451	WRITE(numout,*) ' Record = ', jobs, &
1452	& ' kt = ', kt, &
1453	& ' mstp = ', sstdatqc%mstp(jobs), &
1454	& ' ntyp = ', sstdatqc%ntyp(jobs)
1455	ENDIF
1456	CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
1457
1458	ENDIF
1459
1460	zlam = sstdatqc%rlam(jobs)
1461	zphi = sstdatqc%rphi(jobs)
1462
1463	! Get weights to interpolate the model SST to the observation point
1464	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1465	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1466	& zmask(:,:,iobs), zweig, zobsmask )
1467
1468	! Interpolate the model SST to the observation point
1469
1470	IF ( ld_nightav ) THEN
1471
1472	IF ( idayend == 0 ) THEN
1473	! Daily averaged/diurnal cycle of SST data
1474	CALL obs_int_h2d( 1, 1, &
1475	& zweig, zsstm(:,:,iobs), zext )
1476	ELSE
1477	CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // &
1478	& ' number of night SST data should' // &
1479	& ' only occur at the end of a given day' )
1480	ENDIF
1481
1482	ELSE
1483
1484	CALL obs_int_h2d( 1, 1, &
1485	& zweig, zsstl(:,:,iobs), zext )
1486
1487	ENDIF
1488	sstdatqc%rmod(jobs,1) = zext(1)
1489
1490	END DO
1491
1492	! Deallocate the data for interpolation
1493	DEALLOCATE( &
1494	& igrdi, &
1495	& igrdj, &
1496	& zglam, &
1497	& zgphi, &
1498	& zmask, &
1499	& zsstl &
1500	& )
1501
1502	! At the end of the day also get interpolated means
1503	IF ( idayend == 0 .AND. ld_nightav ) THEN
1504	DEALLOCATE( &
1505	& zsstm &
1506	& )
1507	ENDIF
1508
1509	sstdatqc%nsurfup = sstdatqc%nsurfup + isst
1510
1511	END SUBROUTINE obs_sst_opt
1512
1513	SUBROUTINE obs_sss_opt
1514	!!-----------------------------------------------------------------------
1515	!!
1516	!! * ROUTINE obs_sss_opt *
1517	!!
1518	!! ** Purpose : Compute the model counterpart of sea surface salinity
1519	!! data by interpolating from the model grid to the
1520	!! observation point.
1521	!!
1522	!! ** Method :
1523	!!
1524	!! ** Action :
1525	!!
1526	!! History :
1527	!! ! ??-??
1528	!!-----------------------------------------------------------------------
1529
1530	IMPLICIT NONE
1531
1532	END SUBROUTINE obs_sss_opt
1533
1534	SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
1535	& pseaicen, pseaicemask, k2dint )
1536
1537	!!-----------------------------------------------------------------------
1538	!!
1539	!! * ROUTINE obs_seaice_opt *
1540	!!
1541	!! ** Purpose : Compute the model counterpart of surface temperature
1542	!! data by interpolating from the model grid to the
1543	!! observation point.
1544	!!
1545	!! ** Method : Linearly interpolate to each observation point using
1546	!! the model values at the corners of the surrounding grid box.
1547	!!
1548	!! The now model sea ice is first computed at the obs (lon, lat) point.
1549	!!
1550	!! Several horizontal interpolation schemes are available:
1551	!! - distance-weighted (great circle) (k2dint = 0)
1552	!! - distance-weighted (small angle) (k2dint = 1)
1553	!! - bilinear (geographical grid) (k2dint = 2)
1554	!! - bilinear (quadrilateral grid) (k2dint = 3)
1555	!! - polynomial (quadrilateral grid) (k2dint = 4)
1556	!!
1557	!!
1558	!! ** Action :
1559	!!
1560	!! History :
1561	!! ! 07-07 (S. Ricci ) : Original
1562	!!
1563	!!-----------------------------------------------------------------------
1564
1565	!! * Modules used
1566	USE obs_surf_def ! Definition of storage space for surface observations
1567
1568	IMPLICIT NONE
1569
1570	!! * Arguments
1571	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening
1572	INTEGER, INTENT(IN) :: kt ! Time step
1573	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1574	INTEGER, INTENT(IN) :: kpj
1575	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1576	! (kit000-1 = restart time)
1577	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1578	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1579	& pseaicen, & ! Model sea ice field
1580	& pseaicemask ! Land-sea mask
1581
1582	!! * Local declarations
1583	INTEGER :: ji
1584	INTEGER :: jj
1585	INTEGER :: jobs
1586	INTEGER :: inrc
1587	INTEGER :: iseaice
1588	INTEGER :: iobs
1589
1590	REAL(KIND=wp) :: zlam
1591	REAL(KIND=wp) :: zphi
1592	REAL(KIND=wp) :: zext(1), zobsmask(1)
1593	REAL(kind=wp), DIMENSION(2,2,1) :: &
1594	& zweig
1595	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1596	& zmask, &
1597	& zseaicel, &
1598	& zglam, &
1599	& zgphi
1600	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1601	& igrdi, &
1602	& igrdj
1603
1604	!------------------------------------------------------------------------
1605	! Local initialization
1606	!------------------------------------------------------------------------
1607	! ... Record and data counters
1608	inrc = kt - kit000 + 2
1609	iseaice = seaicedatqc%nsstp(inrc)
1610
1611	! Get the data for interpolation
1612
1613	ALLOCATE( &
1614	& igrdi(2,2,iseaice), &
1615	& igrdj(2,2,iseaice), &
1616	& zglam(2,2,iseaice), &
1617	& zgphi(2,2,iseaice), &
1618	& zmask(2,2,iseaice), &
1619	& zseaicel(2,2,iseaice) &
1620	& )
1621
1622	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1623	iobs = jobs - seaicedatqc%nsurfup
1624	igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1625	igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1626	igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1627	igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1628	igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1629	igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1630	igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1631	igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1632	END DO
1633
1634	CALL obs_int_comm_2d( 2, 2, iseaice, &
1635	& igrdi, igrdj, glamt, zglam )
1636	CALL obs_int_comm_2d( 2, 2, iseaice, &
1637	& igrdi, igrdj, gphit, zgphi )
1638	CALL obs_int_comm_2d( 2, 2, iseaice, &
1639	& igrdi, igrdj, pseaicemask, zmask )
1640	CALL obs_int_comm_2d( 2, 2, iseaice, &
1641	& igrdi, igrdj, pseaicen, zseaicel )
1642
1643	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1644
1645	iobs = jobs - seaicedatqc%nsurfup
1646
1647	IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1648
1649	IF(lwp) THEN
1650	WRITE(numout,*)
1651	WRITE(numout,*) ' E R R O R : Observation', &
1652	& ' time step is not consistent with the', &
1653	& ' model time step'
1654	WRITE(numout,*) ' ========='
1655	WRITE(numout,*)
1656	WRITE(numout,*) ' Record = ', jobs, &
1657	& ' kt = ', kt, &
1658	& ' mstp = ', seaicedatqc%mstp(jobs), &
1659	& ' ntyp = ', seaicedatqc%ntyp(jobs)
1660	ENDIF
1661	CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1662
1663	ENDIF
1664
1665	zlam = seaicedatqc%rlam(jobs)
1666	zphi = seaicedatqc%rphi(jobs)
1667
1668	! Get weights to interpolate the model sea ice to the observation point
1669	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1670	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1671	& zmask(:,:,iobs), zweig, zobsmask )
1672
1673	! ... Interpolate the model sea ice to the observation point
1674	CALL obs_int_h2d( 1, 1, &
1675	& zweig, zseaicel(:,:,iobs), zext )
1676
1677	seaicedatqc%rmod(jobs,1) = zext(1)
1678
1679	END DO
1680
1681	! Deallocate the data for interpolation
1682	DEALLOCATE( &
1683	& igrdi, &
1684	& igrdj, &
1685	& zglam, &
1686	& zgphi, &
1687	& zmask, &
1688	& zseaicel &
1689	& )
1690
1691	seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1692
1693	END SUBROUTINE obs_seaice_opt
1694
1695	SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1696	& pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1697	& ld_dailyav )
1698	!!-----------------------------------------------------------------------
1699	!!
1700	!! * ROUTINE obs_vel_opt *
1701	!!
1702	!! ** Purpose : Compute the model counterpart of velocity profile
1703	!! data by interpolating from the model grid to the
1704	!! observation point.
1705	!!
1706	!! ** Method : Linearly interpolate zonal and meridional components of velocity
1707	!! to each observation point using the model values at the corners of
1708	!! the surrounding grid box. The model velocity components are on a
1709	!! staggered C- grid.
1710	!!
1711	!! For velocity data from the TAO array, the model equivalent is
1712	!! a daily mean velocity field. So, we first compute
1713	!! the mean, then interpolate only at the end of the day.
1714	!!
1715	!! ** Action :
1716	!!
1717	!! History :
1718	!! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles
1719	!! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1720	!!-----------------------------------------------------------------------
1721
1722	!! * Modules used
1723	USE obs_profiles_def ! Definition of storage space for profile obs.
1724
1725	IMPLICIT NONE
1726
1727	!! * Arguments
1728	TYPE(obs_prof), INTENT(INOUT) :: &
1729	& prodatqc ! Subset of profile data not failing screening
1730	INTEGER, INTENT(IN) :: kt ! Time step
1731	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1732	INTEGER, INTENT(IN) :: kpj
1733	INTEGER, INTENT(IN) :: kpk
1734	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1735	! (kit000-1 = restart time)
1736	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
1737	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1738	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1739	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1740	& pun, & ! Model zonal component of velocity
1741	& pvn, & ! Model meridional component of velocity
1742	& pumask, & ! Land-sea mask
1743	& pvmask ! Land-sea mask
1744	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1745	& pgdept ! Model array of depth levels
1746	LOGICAL, INTENT(IN) :: ld_dailyav
1747
1748	!! * Local declarations
1749	INTEGER :: ji
1750	INTEGER :: jj
1751	INTEGER :: jk
1752	INTEGER :: jobs
1753	INTEGER :: inrc
1754	INTEGER :: ipro
1755	INTEGER :: idayend
1756	INTEGER :: ista
1757	INTEGER :: iend
1758	INTEGER :: iobs
1759	INTEGER, DIMENSION(imaxavtypes) :: &
1760	& idailyavtypes
1761	REAL(KIND=wp) :: zlam
1762	REAL(KIND=wp) :: zphi
1763	REAL(KIND=wp) :: zdaystp
1764	REAL(KIND=wp), DIMENSION(kpk) :: &
1765	& zobsmasku, &
1766	& zobsmaskv, &
1767	& zobsmask, &
1768	& zobsk, &
1769	& zobs2k
1770	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1771	& zweigu,zweigv
1772	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1773	& zumask, zvmask, &
1774	& zintu, &
1775	& zintv, &
1776	& zinmu, &
1777	& zinmv
1778	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1779	& zglamu, zglamv, &
1780	& zgphiu, zgphiv
1781	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1782	& igrdiu, &
1783	& igrdju, &
1784	& igrdiv, &
1785	& igrdjv
1786
1787	!------------------------------------------------------------------------
1788	! Local initialization
1789	!------------------------------------------------------------------------
1790	! ... Record and data counters
1791	inrc = kt - kit000 + 2
1792	ipro = prodatqc%npstp(inrc)
1793
1794	! Initialize daily mean for first timestep
1795	idayend = MOD( kt - kit000 + 1, kdaystp )
1796
1797	! Added kt == 0 test to catch restart case
1798	IF ( idayend == 1 .OR. kt == 0) THEN
1799	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1800	prodatqc%vdmean(:,:,:,1) = 0.0
1801	prodatqc%vdmean(:,:,:,2) = 0.0
1802	ENDIF
1803
1804	! Increment the zonal velocity field for computing daily mean
1805	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1806	! Increment the meridional velocity field for computing daily mean
1807	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1808
1809	! Compute the daily mean at the end of day
1810	zdaystp = 1.0 / REAL( kdaystp )
1811	IF ( idayend == 0 ) THEN
1812	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1813	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1814	ENDIF
1815
1816	! Get the data for interpolation
1817	ALLOCATE( &
1818	& igrdiu(2,2,ipro), &
1819	& igrdju(2,2,ipro), &
1820	& igrdiv(2,2,ipro), &
1821	& igrdjv(2,2,ipro), &
1822	& zglamu(2,2,ipro), zglamv(2,2,ipro), &
1823	& zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1824	& zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1825	& zintu(2,2,kpk,ipro), &
1826	& zintv(2,2,kpk,ipro) &
1827	& )
1828
1829	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1830	iobs = jobs - prodatqc%nprofup
1831	igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1832	igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1833	igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1834	igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1835	igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1836	igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1837	igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1838	igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1839	igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1840	igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1841	igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1842	igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1843	igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1844	igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1845	igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1846	igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1847	END DO
1848
1849	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1850	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1851	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1852	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1853
1854	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1855	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1856	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1857	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1858
1859	! At the end of the day also get interpolated means
1860	IF ( idayend == 0 ) THEN
1861
1862	ALLOCATE( &
1863	& zinmu(2,2,kpk,ipro), &
1864	& zinmv(2,2,kpk,ipro) &
1865	& )
1866
1867	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1868	& prodatqc%vdmean(:,:,:,1), zinmu )
1869	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1870	& prodatqc%vdmean(:,:,:,2), zinmv )
1871
1872	ENDIF
1873
1874	! loop over observations
1875
1876	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1877
1878	iobs = jobs - prodatqc%nprofup
1879
1880	IF ( kt /= prodatqc%mstp(jobs) ) THEN
1881
1882	IF(lwp) THEN
1883	WRITE(numout,*)
1884	WRITE(numout,*) ' E R R O R : Observation', &
1885	& ' time step is not consistent with the', &
1886	& ' model time step'
1887	WRITE(numout,*) ' ========='
1888	WRITE(numout,*)
1889	WRITE(numout,*) ' Record = ', jobs, &
1890	& ' kt = ', kt, &
1891	& ' mstp = ', prodatqc%mstp(jobs), &
1892	& ' ntyp = ', prodatqc%ntyp(jobs)
1893	ENDIF
1894	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1895	ENDIF
1896
1897	zlam = prodatqc%rlam(jobs)
1898	zphi = prodatqc%rphi(jobs)
1899
1900	! Initialize observation masks
1901
1902	zobsmasku(:) = 0.0
1903	zobsmaskv(:) = 0.0
1904
1905	! Horizontal weights and vertical mask
1906
1907	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1908
1909	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1910	& zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1911	& zumask(:,:,:,iobs), zweigu, zobsmasku )
1912
1913	ENDIF
1914
1915
1916	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1917
1918	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1919	& zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1920	& zvmask(:,:,:,iobs), zweigv, zobsmaskv )
1921
1922	ENDIF
1923
1924	! Ensure that the vertical mask on u and v are consistent.
1925
1926	zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1927
1928	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1929
1930	zobsk(:) = obfillflt
1931
1932	IF ( ld_dailyav ) THEN
1933
1934	IF ( idayend == 0 ) THEN
1935
1936	! Daily averaged data
1937
1938	CALL obs_int_h2d( kpk, kpk, &
1939	& zweigu, zinmu(:,:,:,iobs), zobsk )
1940
1941
1942	ELSE
1943
1944	CALL ctl_stop( ' A nonzero' // &
1945	& ' number of U profile data should' // &
1946	& ' only occur at the end of a given day' )
1947
1948	ENDIF
1949
1950	ELSE
1951
1952	! Point data
1953
1954	CALL obs_int_h2d( kpk, kpk, &
1955	& zweigu, zintu(:,:,:,iobs), zobsk )
1956
1957	ENDIF
1958
1959	!-------------------------------------------------------------
1960	! Compute vertical second-derivative of the interpolating
1961	! polynomial at obs points
1962	!-------------------------------------------------------------
1963
1964	IF ( k1dint == 1 ) THEN
1965	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1966	& pgdept, zobsmask )
1967	ENDIF
1968
1969	!-----------------------------------------------------------------
1970	! Vertical interpolation to the observation point
1971	!-----------------------------------------------------------------
1972	ista = prodatqc%npvsta(jobs,1)
1973	iend = prodatqc%npvend(jobs,1)
1974	CALL obs_int_z1d( kpk, &
1975	& prodatqc%var(1)%mvk(ista:iend), &
1976	& k1dint, iend - ista + 1, &
1977	& prodatqc%var(1)%vdep(ista:iend), &
1978	& zobsk, zobs2k, &
1979	& prodatqc%var(1)%vmod(ista:iend), &
1980	& pgdept, zobsmask )
1981
1982	ENDIF
1983
1984	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1985
1986	zobsk(:) = obfillflt
1987
1988	IF ( ld_dailyav ) THEN
1989
1990	IF ( idayend == 0 ) THEN
1991
1992	! Daily averaged data
1993
1994	CALL obs_int_h2d( kpk, kpk, &
1995	& zweigv, zinmv(:,:,:,iobs), zobsk )
1996
1997	ELSE
1998
1999	CALL ctl_stop( ' A nonzero' // &
2000	& ' number of V profile data should' // &
2001	& ' only occur at the end of a given day' )
2002
2003	ENDIF
2004
2005	ELSE
2006
2007	! Point data
2008
2009	CALL obs_int_h2d( kpk, kpk, &
2010	& zweigv, zintv(:,:,:,iobs), zobsk )
2011
2012	ENDIF
2013
2014
2015	!-------------------------------------------------------------
2016	! Compute vertical second-derivative of the interpolating
2017	! polynomial at obs points
2018	!-------------------------------------------------------------
2019
2020	IF ( k1dint == 1 ) THEN
2021	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
2022	& pgdept, zobsmask )
2023	ENDIF
2024
2025	!----------------------------------------------------------------
2026	! Vertical interpolation to the observation point
2027	!----------------------------------------------------------------
2028	ista = prodatqc%npvsta(jobs,2)
2029	iend = prodatqc%npvend(jobs,2)
2030	CALL obs_int_z1d( kpk, &
2031	& prodatqc%var(2)%mvk(ista:iend),&
2032	& k1dint, iend - ista + 1, &
2033	& prodatqc%var(2)%vdep(ista:iend),&
2034	& zobsk, zobs2k, &
2035	& prodatqc%var(2)%vmod(ista:iend),&
2036	& pgdept, zobsmask )
2037
2038	ENDIF
2039
2040	END DO
2041
2042	! Deallocate the data for interpolation
2043	DEALLOCATE( &
2044	& igrdiu, &
2045	& igrdju, &
2046	& igrdiv, &
2047	& igrdjv, &
2048	& zglamu, zglamv, &
2049	& zgphiu, zgphiv, &
2050	& zumask, zvmask, &
2051	& zintu, &
2052	& zintv &
2053	& )
2054	! At the end of the day also get interpolated means
2055	IF ( idayend == 0 ) THEN
2056	DEALLOCATE( &
2057	& zinmu, &
2058	& zinmv &
2059	& )
2060	ENDIF
2061
2062	prodatqc%nprofup = prodatqc%nprofup + ipro
2063
2064	END SUBROUTINE obs_vel_opt
2065
2066	SUBROUTINE obs_logchl_opt( logchldatqc, kt, kpi, kpj, kit000, &
2067	& plogchln, plogchlmask, k2dint )
2068
2069	!!-----------------------------------------------------------------------
2070	!!
2071	!! * ROUTINE obs_logchl_opt *
2072	!!
2073	!! ** Purpose : Compute the model counterpart of logchl
2074	!! data by interpolating from the model grid to the
2075	!! observation point.
2076	!!
2077	!! ** Method : Linearly interpolate to each observation point using
2078	!! the model values at the corners of the surrounding grid box.
2079	!!
2080	!! The now model logchl is first computed at the obs (lon, lat) point.
2081	!!
2082	!! Several horizontal interpolation schemes are available:
2083	!! - distance-weighted (great circle) (k2dint = 0)
2084	!! - distance-weighted (small angle) (k2dint = 1)
2085	!! - bilinear (geographical grid) (k2dint = 2)
2086	!! - bilinear (quadrilateral grid) (k2dint = 3)
2087	!! - polynomial (quadrilateral grid) (k2dint = 4)
2088	!!
2089	!!
2090	!! ** Action :
2091	!!
2092	!! History :
2093	!!
2094	!!-----------------------------------------------------------------------
2095
2096	!! * Modules used
2097	USE obs_surf_def ! Definition of storage space for surface observations
2098
2099	IMPLICIT NONE
2100
2101	!! * Arguments
2102	TYPE(obs_surf), INTENT(INOUT) :: logchldatqc ! Subset of surface data not failing screening
2103	INTEGER, INTENT(IN) :: kt ! Time step
2104	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
2105	INTEGER, INTENT(IN) :: kpj
2106	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
2107	! (kit000-1 = restart time)
2108	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
2109	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
2110	& plogchln, & ! Model logchl field
2111	& plogchlmask ! Land-sea mask
2112
2113	!! * Local declarations
2114	INTEGER :: ji
2115	INTEGER :: jj
2116	INTEGER :: jobs
2117	INTEGER :: inrc
2118	INTEGER :: ilogchl
2119	INTEGER :: iobs
2120
2121	REAL(KIND=wp) :: zlam
2122	REAL(KIND=wp) :: zphi
2123	REAL(KIND=wp) :: zext(1), zobsmask(1)
2124	REAL(kind=wp), DIMENSION(2,2,1) :: &
2125	& zweig
2126	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2127	& zmask, &
2128	& zlogchll, &
2129	& zglam, &
2130	& zgphi
2131	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
2132	& igrdi, &
2133	& igrdj
2134
2135	!------------------------------------------------------------------------
2136	! Local initialization
2137	!------------------------------------------------------------------------
2138	! ... Record and data counters
2139	inrc = kt - kit000 + 2
2140	ilogchl = logchldatqc%nsstp(inrc)
2141
2142	! Get the data for interpolation
2143
2144	ALLOCATE( &
2145	& igrdi(2,2,ilogchl), &
2146	& igrdj(2,2,ilogchl), &
2147	& zglam(2,2,ilogchl), &
2148	& zgphi(2,2,ilogchl), &
2149	& zmask(2,2,ilogchl), &
2150	& zlogchll(2,2,ilogchl) &
2151	& )
2152
2153	DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl
2154	iobs = jobs - logchldatqc%nsurfup
2155	igrdi(1,1,iobs) = logchldatqc%mi(jobs)-1
2156	igrdj(1,1,iobs) = logchldatqc%mj(jobs)-1
2157	igrdi(1,2,iobs) = logchldatqc%mi(jobs)-1
2158	igrdj(1,2,iobs) = logchldatqc%mj(jobs)
2159	igrdi(2,1,iobs) = logchldatqc%mi(jobs)
2160	igrdj(2,1,iobs) = logchldatqc%mj(jobs)-1
2161	igrdi(2,2,iobs) = logchldatqc%mi(jobs)
2162	igrdj(2,2,iobs) = logchldatqc%mj(jobs)
2163	END DO
2164
2165	CALL obs_int_comm_2d( 2, 2, ilogchl, &
2166	& igrdi, igrdj, glamt, zglam )
2167	CALL obs_int_comm_2d( 2, 2, ilogchl, &
2168	& igrdi, igrdj, gphit, zgphi )
2169	CALL obs_int_comm_2d( 2, 2, ilogchl, &
2170	& igrdi, igrdj, plogchlmask, zmask )
2171	CALL obs_int_comm_2d( 2, 2, ilogchl, &
2172	& igrdi, igrdj, plogchln, zlogchll )
2173
2174	DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl
2175
2176	iobs = jobs - logchldatqc%nsurfup
2177
2178	IF ( kt /= logchldatqc%mstp(jobs) ) THEN
2179
2180	IF(lwp) THEN
2181	WRITE(numout,*)
2182	WRITE(numout,*) ' E R R O R : Observation', &
2183	& ' time step is not consistent with the', &
2184	& ' model time step'
2185	WRITE(numout,*) ' ========='
2186	WRITE(numout,*)
2187	WRITE(numout,*) ' Record = ', jobs, &
2188	& ' kt = ', kt, &
2189	& ' mstp = ', logchldatqc%mstp(jobs), &
2190	& ' ntyp = ', logchldatqc%ntyp(jobs)
2191	ENDIF
2192	CALL ctl_stop( 'obs_logchl_opt', 'Inconsistent time' )
2193
2194	ENDIF
2195
2196	zlam = logchldatqc%rlam(jobs)
2197	zphi = logchldatqc%rphi(jobs)
2198
2199	! Get weights to interpolate the model logchl to the observation point
2200	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
2201	& zglam(:,:,iobs), zgphi(:,:,iobs), &
2202	& zmask(:,:,iobs), zweig, zobsmask )
2203
2204	! ... Interpolate the model logchl to the observation point
2205	CALL obs_int_h2d( 1, 1, &
2206	& zweig, zlogchll(:,:,iobs), zext )
2207
2208	logchldatqc%rmod(jobs,1) = zext(1)
2209
2210	END DO
2211
2212	! Deallocate the data for interpolation
2213	DEALLOCATE( &
2214	& igrdi, &
2215	& igrdj, &
2216	& zglam, &
2217	& zgphi, &
2218	& zmask, &
2219	& zlogchll &
2220	& )
2221
2222	logchldatqc%nsurfup = logchldatqc%nsurfup + ilogchl
2223
2224	END SUBROUTINE obs_logchl_opt
2225
2226	SUBROUTINE obs_spm_opt( spmdatqc, kt, kpi, kpj, kit000, &
2227	& pspmn, pspmmask, k2dint )
2228
2229	!!-----------------------------------------------------------------------
2230	!!
2231	!! * ROUTINE obs_spm_opt *
2232	!!
2233	!! ** Purpose : Compute the model counterpart of spm
2234	!! data by interpolating from the model grid to the
2235	!! observation point.
2236	!!
2237	!! ** Method : Linearly interpolate to each observation point using
2238	!! the model values at the corners of the surrounding grid box.
2239	!!
2240	!! The now model spm is first computed at the obs (lon, lat) point.
2241	!!
2242	!! Several horizontal interpolation schemes are available:
2243	!! - distance-weighted (great circle) (k2dint = 0)
2244	!! - distance-weighted (small angle) (k2dint = 1)
2245	!! - bilinear (geographical grid) (k2dint = 2)
2246	!! - bilinear (quadrilateral grid) (k2dint = 3)
2247	!! - polynomial (quadrilateral grid) (k2dint = 4)
2248	!!
2249	!!
2250	!! ** Action :
2251	!!
2252	!! History :
2253	!!
2254	!!-----------------------------------------------------------------------
2255
2256	!! * Modules used
2257	USE obs_surf_def ! Definition of storage space for surface observations
2258
2259	IMPLICIT NONE
2260
2261	!! * Arguments
2262	TYPE(obs_surf), INTENT(INOUT) :: spmdatqc ! Subset of surface data not failing screening
2263	INTEGER, INTENT(IN) :: kt ! Time step
2264	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
2265	INTEGER, INTENT(IN) :: kpj
2266	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
2267	! (kit000-1 = restart time)
2268	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
2269	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
2270	& pspmn, & ! Model spm field
2271	& pspmmask ! Land-sea mask
2272
2273	!! * Local declarations
2274	INTEGER :: ji
2275	INTEGER :: jj
2276	INTEGER :: jobs
2277	INTEGER :: inrc
2278	INTEGER :: ispm
2279	INTEGER :: iobs
2280
2281	REAL(KIND=wp) :: zlam
2282	REAL(KIND=wp) :: zphi
2283	REAL(KIND=wp) :: zext(1), zobsmask(1)
2284	REAL(kind=wp), DIMENSION(2,2,1) :: &
2285	& zweig
2286	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2287	& zmask, &
2288	& zspml, &
2289	& zglam, &
2290	& zgphi
2291	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
2292	& igrdi, &
2293	& igrdj
2294
2295	!------------------------------------------------------------------------
2296	! Local initialization
2297	!------------------------------------------------------------------------
2298	! ... Record and data counters
2299	inrc = kt - kit000 + 2
2300	ispm = spmdatqc%nsstp(inrc)
2301
2302	! Get the data for interpolation
2303
2304	ALLOCATE( &
2305	& igrdi(2,2,ispm), &
2306	& igrdj(2,2,ispm), &
2307	& zglam(2,2,ispm), &
2308	& zgphi(2,2,ispm), &
2309	& zmask(2,2,ispm), &
2310	& zspml(2,2,ispm) &
2311	& )
2312
2313	DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm
2314	iobs = jobs - spmdatqc%nsurfup
2315	igrdi(1,1,iobs) = spmdatqc%mi(jobs)-1
2316	igrdj(1,1,iobs) = spmdatqc%mj(jobs)-1
2317	igrdi(1,2,iobs) = spmdatqc%mi(jobs)-1
2318	igrdj(1,2,iobs) = spmdatqc%mj(jobs)
2319	igrdi(2,1,iobs) = spmdatqc%mi(jobs)
2320	igrdj(2,1,iobs) = spmdatqc%mj(jobs)-1
2321	igrdi(2,2,iobs) = spmdatqc%mi(jobs)
2322	igrdj(2,2,iobs) = spmdatqc%mj(jobs)
2323	END DO
2324
2325	CALL obs_int_comm_2d( 2, 2, ispm, &
2326	& igrdi, igrdj, glamt, zglam )
2327	CALL obs_int_comm_2d( 2, 2, ispm, &
2328	& igrdi, igrdj, gphit, zgphi )
2329	CALL obs_int_comm_2d( 2, 2, ispm, &
2330	& igrdi, igrdj, pspmmask, zmask )
2331	CALL obs_int_comm_2d( 2, 2, ispm, &
2332	& igrdi, igrdj, pspmn, zspml )
2333
2334	DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm
2335
2336	iobs = jobs - spmdatqc%nsurfup
2337
2338	IF ( kt /= spmdatqc%mstp(jobs) ) THEN
2339
2340	IF(lwp) THEN
2341	WRITE(numout,*)
2342	WRITE(numout,*) ' E R R O R : Observation', &
2343	& ' time step is not consistent with the', &
2344	& ' model time step'
2345	WRITE(numout,*) ' ========='
2346	WRITE(numout,*)
2347	WRITE(numout,*) ' Record = ', jobs, &
2348	& ' kt = ', kt, &
2349	& ' mstp = ', spmdatqc%mstp(jobs), &
2350	& ' ntyp = ', spmdatqc%ntyp(jobs)
2351	ENDIF
2352	CALL ctl_stop( 'obs_spm_opt', 'Inconsistent time' )
2353
2354	ENDIF
2355
2356	zlam = spmdatqc%rlam(jobs)
2357	zphi = spmdatqc%rphi(jobs)
2358
2359	! Get weights to interpolate the model spm to the observation point
2360	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
2361	& zglam(:,:,iobs), zgphi(:,:,iobs), &
2362	& zmask(:,:,iobs), zweig, zobsmask )
2363
2364	! ... Interpolate the model spm to the observation point
2365	CALL obs_int_h2d( 1, 1, &
2366	& zweig, zspml(:,:,iobs), zext )
2367
2368	spmdatqc%rmod(jobs,1) = zext(1)
2369
2370	END DO
2371
2372	! Deallocate the data for interpolation
2373	DEALLOCATE( &
2374	& igrdi, &
2375	& igrdj, &
2376	& zglam, &
2377	& zgphi, &
2378	& zmask, &
2379	& zspml &
2380	& )
2381
2382	spmdatqc%nsurfup = spmdatqc%nsurfup + ispm
2383
2384	END SUBROUTINE obs_spm_opt
2385
2386	SUBROUTINE obs_fco2_opt( fco2datqc, kt, kpi, kpj, kit000, &
2387	& pfco2n, pfco2mask, k2dint )
2388
2389	!!-----------------------------------------------------------------------
2390	!!
2391	!! * ROUTINE obs_fco2_opt *
2392	!!
2393	!! ** Purpose : Compute the model counterpart of fco2
2394	!! data by interpolating from the model grid to the
2395	!! observation point.
2396	!!
2397	!! ** Method : Linearly interpolate to each observation point using
2398	!! the model values at the corners of the surrounding grid box.
2399	!!
2400	!! The now model fco2 is first computed at the obs (lon, lat) point.
2401	!!
2402	!! Several horizontal interpolation schemes are available:
2403	!! - distance-weighted (great circle) (k2dint = 0)
2404	!! - distance-weighted (small angle) (k2dint = 1)
2405	!! - bilinear (geographical grid) (k2dint = 2)
2406	!! - bilinear (quadrilateral grid) (k2dint = 3)
2407	!! - polynomial (quadrilateral grid) (k2dint = 4)
2408	!!
2409	!!
2410	!! ** Action :
2411	!!
2412	!! History :
2413	!!
2414	!!-----------------------------------------------------------------------
2415
2416	!! * Modules used
2417	USE obs_surf_def ! Definition of storage space for surface observations
2418
2419	IMPLICIT NONE
2420
2421	!! * Arguments
2422	TYPE(obs_surf), INTENT(INOUT) :: fco2datqc ! Subset of surface data not failing screening
2423	INTEGER, INTENT(IN) :: kt ! Time step
2424	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
2425	INTEGER, INTENT(IN) :: kpj
2426	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
2427	! (kit000-1 = restart time)
2428	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
2429	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
2430	& pfco2n, & ! Model fco2 field
2431	& pfco2mask ! Land-sea mask
2432
2433	!! * Local declarations
2434	INTEGER :: ji
2435	INTEGER :: jj
2436	INTEGER :: jobs
2437	INTEGER :: inrc
2438	INTEGER :: ifco2
2439	INTEGER :: iobs
2440
2441	REAL(KIND=wp) :: zlam
2442	REAL(KIND=wp) :: zphi
2443	REAL(KIND=wp) :: zext(1), zobsmask(1)
2444	REAL(kind=wp), DIMENSION(2,2,1) :: &
2445	& zweig
2446	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2447	& zmask, &
2448	& zfco2l, &
2449	& zglam, &
2450	& zgphi
2451	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
2452	& igrdi, &
2453	& igrdj
2454
2455	!------------------------------------------------------------------------
2456	! Local initialization
2457	!------------------------------------------------------------------------
2458	! ... Record and data counters
2459	inrc = kt - kit000 + 2
2460	ifco2 = fco2datqc%nsstp(inrc)
2461
2462	! Get the data for interpolation
2463
2464	ALLOCATE( &
2465	& igrdi(2,2,ifco2), &
2466	& igrdj(2,2,ifco2), &
2467	& zglam(2,2,ifco2), &
2468	& zgphi(2,2,ifco2), &
2469	& zmask(2,2,ifco2), &
2470	& zfco2l(2,2,ifco2) &
2471	& )
2472
2473	DO jobs = fco2datqc%nsurfup + 1, fco2datqc%nsurfup + ifco2
2474	iobs = jobs - fco2datqc%nsurfup
2475	igrdi(1,1,iobs) = fco2datqc%mi(jobs)-1
2476	igrdj(1,1,iobs) = fco2datqc%mj(jobs)-1
2477	igrdi(1,2,iobs) = fco2datqc%mi(jobs)-1
2478	igrdj(1,2,iobs) = fco2datqc%mj(jobs)
2479	igrdi(2,1,iobs) = fco2datqc%mi(jobs)
2480	igrdj(2,1,iobs) = fco2datqc%mj(jobs)-1
2481	igrdi(2,2,iobs) = fco2datqc%mi(jobs)
2482	igrdj(2,2,iobs) = fco2datqc%mj(jobs)
2483	END DO
2484
2485	CALL obs_int_comm_2d( 2, 2, ifco2, &
2486	& igrdi, igrdj, glamt, zglam )
2487	CALL obs_int_comm_2d( 2, 2, ifco2, &
2488	& igrdi, igrdj, gphit, zgphi )
2489	CALL obs_int_comm_2d( 2, 2, ifco2, &
2490	& igrdi, igrdj, pfco2mask, zmask )
2491	CALL obs_int_comm_2d( 2, 2, ifco2, &
2492	& igrdi, igrdj, pfco2n, zfco2l )
2493
2494	DO jobs = fco2datqc%nsurfup + 1, fco2datqc%nsurfup + ifco2
2495
2496	iobs = jobs - fco2datqc%nsurfup
2497
2498	IF ( kt /= fco2datqc%mstp(jobs) ) THEN
2499
2500	IF(lwp) THEN
2501	WRITE(numout,*)
2502	WRITE(numout,*) ' E R R O R : Observation', &
2503	& ' time step is not consistent with the', &
2504	& ' model time step'
2505	WRITE(numout,*) ' ========='
2506	WRITE(numout,*)
2507	WRITE(numout,*) ' Record = ', jobs, &
2508	& ' kt = ', kt, &
2509	& ' mstp = ', fco2datqc%mstp(jobs), &
2510	& ' ntyp = ', fco2datqc%ntyp(jobs)
2511	ENDIF
2512	CALL ctl_stop( 'obs_fco2_opt', 'Inconsistent time' )
2513
2514	ENDIF
2515
2516	zlam = fco2datqc%rlam(jobs)
2517	zphi = fco2datqc%rphi(jobs)
2518
2519	! Get weights to interpolate the model fco2 to the observation point
2520	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
2521	& zglam(:,:,iobs), zgphi(:,:,iobs), &
2522	& zmask(:,:,iobs), zweig, zobsmask )
2523
2524	! ... Interpolate the model fco2 to the observation point
2525	CALL obs_int_h2d( 1, 1, &
2526	& zweig, zfco2l(:,:,iobs), zext )
2527
2528	fco2datqc%rmod(jobs,1) = zext(1)
2529
2530	END DO
2531
2532	! Deallocate the data for interpolation
2533	DEALLOCATE( &
2534	& igrdi, &
2535	& igrdj, &
2536	& zglam, &
2537	& zgphi, &
2538	& zmask, &
2539	& zfco2l &
2540	& )
2541
2542	fco2datqc%nsurfup = fco2datqc%nsurfup + ifco2
2543
2544	END SUBROUTINE obs_fco2_opt
2545
2546	SUBROUTINE obs_pco2_opt( pco2datqc, kt, kpi, kpj, kit000, &
2547	& ppco2n, ppco2mask, k2dint )
2548
2549	!!-----------------------------------------------------------------------
2550	!!
2551	!! * ROUTINE obs_pco2_opt *
2552	!!
2553	!! ** Purpose : Compute the model counterpart of pco2
2554	!! data by interpolating from the model grid to the
2555	!! observation point.
2556	!!
2557	!! ** Method : Linearly interpolate to each observation point using
2558	!! the model values at the corners of the surrounding grid box.
2559	!!
2560	!! The now model pco2 is first computed at the obs (lon, lat) point.
2561	!!
2562	!! Several horizontal interpolation schemes are available:
2563	!! - distance-weighted (great circle) (k2dint = 0)
2564	!! - distance-weighted (small angle) (k2dint = 1)
2565	!! - bilinear (geographical grid) (k2dint = 2)
2566	!! - bilinear (quadrilateral grid) (k2dint = 3)
2567	!! - polynomial (quadrilateral grid) (k2dint = 4)
2568	!!
2569	!!
2570	!! ** Action :
2571	!!
2572	!! History :
2573	!!
2574	!!-----------------------------------------------------------------------
2575
2576	!! * Modules used
2577	USE obs_surf_def ! Definition of storage space for surface observations
2578
2579	IMPLICIT NONE
2580
2581	!! * Arguments
2582	TYPE(obs_surf), INTENT(INOUT) :: pco2datqc ! Subset of surface data not failing screening
2583	INTEGER, INTENT(IN) :: kt ! Time step
2584	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
2585	INTEGER, INTENT(IN) :: kpj
2586	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
2587	! (kit000-1 = restart time)
2588	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
2589	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
2590	& ppco2n, & ! Model pco2 field
2591	& ppco2mask ! Land-sea mask
2592
2593	!! * Local declarations
2594	INTEGER :: ji
2595	INTEGER :: jj
2596	INTEGER :: jobs
2597	INTEGER :: inrc
2598	INTEGER :: ipco2
2599	INTEGER :: iobs
2600
2601	REAL(KIND=wp) :: zlam
2602	REAL(KIND=wp) :: zphi
2603	REAL(KIND=wp) :: zext(1), zobsmask(1)
2604	REAL(kind=wp), DIMENSION(2,2,1) :: &
2605	& zweig
2606	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2607	& zmask, &
2608	& zpco2l, &
2609	& zglam, &
2610	& zgphi
2611	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
2612	& igrdi, &
2613	& igrdj
2614
2615	!------------------------------------------------------------------------
2616	! Local initialization
2617	!------------------------------------------------------------------------
2618	! ... Record and data counters
2619	inrc = kt - kit000 + 2
2620	ipco2 = pco2datqc%nsstp(inrc)
2621
2622	! Get the data for interpolation
2623
2624	ALLOCATE( &
2625	& igrdi(2,2,ipco2), &
2626	& igrdj(2,2,ipco2), &
2627	& zglam(2,2,ipco2), &
2628	& zgphi(2,2,ipco2), &
2629	& zmask(2,2,ipco2), &
2630	& zpco2l(2,2,ipco2) &
2631	& )
2632
2633	DO jobs = pco2datqc%nsurfup + 1, pco2datqc%nsurfup + ipco2
2634	iobs = jobs - pco2datqc%nsurfup
2635	igrdi(1,1,iobs) = pco2datqc%mi(jobs)-1
2636	igrdj(1,1,iobs) = pco2datqc%mj(jobs)-1
2637	igrdi(1,2,iobs) = pco2datqc%mi(jobs)-1
2638	igrdj(1,2,iobs) = pco2datqc%mj(jobs)
2639	igrdi(2,1,iobs) = pco2datqc%mi(jobs)
2640	igrdj(2,1,iobs) = pco2datqc%mj(jobs)-1
2641	igrdi(2,2,iobs) = pco2datqc%mi(jobs)
2642	igrdj(2,2,iobs) = pco2datqc%mj(jobs)
2643	END DO
2644
2645	CALL obs_int_comm_2d( 2, 2, ipco2, &
2646	& igrdi, igrdj, glamt, zglam )
2647	CALL obs_int_comm_2d( 2, 2, ipco2, &
2648	& igrdi, igrdj, gphit, zgphi )
2649	CALL obs_int_comm_2d( 2, 2, ipco2, &
2650	& igrdi, igrdj, ppco2mask, zmask )
2651	CALL obs_int_comm_2d( 2, 2, ipco2, &
2652	& igrdi, igrdj, ppco2n, zpco2l )
2653
2654	DO jobs = pco2datqc%nsurfup + 1, pco2datqc%nsurfup + ipco2
2655
2656	iobs = jobs - pco2datqc%nsurfup
2657
2658	IF ( kt /= pco2datqc%mstp(jobs) ) THEN
2659
2660	IF(lwp) THEN
2661	WRITE(numout,*)
2662	WRITE(numout,*) ' E R R O R : Observation', &
2663	& ' time step is not consistent with the', &
2664	& ' model time step'
2665	WRITE(numout,*) ' ========='
2666	WRITE(numout,*)
2667	WRITE(numout,*) ' Record = ', jobs, &
2668	& ' kt = ', kt, &
2669	& ' mstp = ', pco2datqc%mstp(jobs), &
2670	& ' ntyp = ', pco2datqc%ntyp(jobs)
2671	ENDIF
2672	CALL ctl_stop( 'obs_pco2_opt', 'Inconsistent time' )
2673
2674	ENDIF
2675
2676	zlam = pco2datqc%rlam(jobs)
2677	zphi = pco2datqc%rphi(jobs)
2678
2679	! Get weights to interpolate the model pco2 to the observation point
2680	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
2681	& zglam(:,:,iobs), zgphi(:,:,iobs), &
2682	& zmask(:,:,iobs), zweig, zobsmask )
2683
2684	! ... Interpolate the model pco2 to the observation point
2685	CALL obs_int_h2d( 1, 1, &
2686	& zweig, zpco2l(:,:,iobs), zext )
2687
2688	pco2datqc%rmod(jobs,1) = zext(1)
2689
2690	END DO
2691
2692	! Deallocate the data for interpolation
2693	DEALLOCATE( &
2694	& igrdi, &
2695	& igrdj, &
2696	& zglam, &
2697	& zgphi, &
2698	& zmask, &
2699	& zpco2l &
2700	& )
2701
2702	pco2datqc%nsurfup = pco2datqc%nsurfup + ipco2
2703
2704	END SUBROUTINE obs_pco2_opt
2705
2706	END MODULE obs_oper
2707

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

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

Context Navigation

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

Download in other formats: