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/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 5979

Last change on this file since 5979 was 5979, checked in by timgraham, 8 years ago
Bug fixes following last merge
Property svn:keywords set to `Id`
File size: 75.5 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	!!----------------------------------------------------------------------
26
27	!! * Modules used
28	USE par_kind, ONLY : & ! Precision variables
29	& wp
30	USE in_out_manager ! I/O manager
31	USE obs_inter_sup ! Interpolation support
32	USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the observation pt
33	& obs_int_h2d, &
34	& obs_int_h2d_init
35	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt
36	& obs_int_z1d, &
37	& obs_int_z1d_spl
38	USE obs_const, ONLY : &
39	& obfillflt ! Fillvalue
40	USE dom_oce, ONLY : &
41	& glamt, glamu, glamv, &
42	& gphit, gphiu, gphiv, &
43	#if defined key_vvl
44	& gdept_n
45	#else
46	& gdept_0
47	#endif
48	USE lib_mpp, ONLY : &
49	& ctl_warn, ctl_stop
50	USE obs_grid, ONLY : &
51	& obs_level_search
52
53	IMPLICIT NONE
54
55	!! * Routine accessibility
56	PRIVATE
57
58	PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations
59	& obs_pro_sco_opt, & ! Compute the model counterpart of profile observations
60	! in generalised vertical coordinates
61	& obs_sla_opt, & ! Compute the model counterpart of SLA observations
62	& obs_sst_opt, & ! Compute the model counterpart of SST observations
63	& obs_sss_opt, & ! Compute the model counterpart of SSS observations
64	& obs_seaice_opt, &
65	& obs_vel_opt ! Compute the model counterpart of velocity profile data
66
67	INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
68
69	!!----------------------------------------------------------------------
70	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
71	!! $Id$
72	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
73	!!----------------------------------------------------------------------
74
75	!! * Substitutions
76	# include "domzgr_substitute.h90"
77	CONTAINS
78
79	SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
80	& ptn, psn, pgdept, ptmask, k1dint, k2dint, &
81	& kdailyavtypes )
82	!!-----------------------------------------------------------------------
83	!!
84	!! * ROUTINE obs_pro_opt *
85	!!
86	!! ** Purpose : Compute the model counterpart of profiles
87	!! data by interpolating from the model grid to the
88	!! observation point.
89	!!
90	!! ** Method : Linearly interpolate to each observation point using
91	!! the model values at the corners of the surrounding grid box.
92	!!
93	!! First, a vertical profile of horizontally interpolated model
94	!! now temperatures is computed at the obs (lon, lat) point.
95	!! Several horizontal interpolation schemes are available:
96	!! - distance-weighted (great circle) (k2dint = 0)
97	!! - distance-weighted (small angle) (k2dint = 1)
98	!! - bilinear (geographical grid) (k2dint = 2)
99	!! - bilinear (quadrilateral grid) (k2dint = 3)
100	!! - polynomial (quadrilateral grid) (k2dint = 4)
101	!!
102	!! Next, the vertical temperature profile is interpolated to the
103	!! data depth points. Two vertical interpolation schemes are
104	!! available:
105	!! - linear (k1dint = 0)
106	!! - Cubic spline (k1dint = 1)
107	!!
108	!! For the cubic spline the 2nd derivative of the interpolating
109	!! polynomial is computed before entering the vertical interpolation
110	!! routine.
111	!!
112	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
113	!! a daily mean model temperature field. So, we first compute
114	!! the mean, then interpolate only at the end of the day.
115	!!
116	!! Note: the in situ temperature observations must be converted
117	!! to potential temperature (the model variable) prior to
118	!! assimilation.
119	!!??????????????????????????????????????????????????????????????
120	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
121	!!??????????????????????????????????????????????????????????????
122	!!
123	!! ** Action :
124	!!
125	!! History :
126	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
127	!! ! 06-03 (G. Smith) NEMOVAR migration
128	!! ! 06-10 (A. Weaver) Cleanup
129	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
130	!! ! 07-03 (K. Mogensen) General handling of profiles
131	!!-----------------------------------------------------------------------
132
133	!! * Modules used
134	USE obs_profiles_def ! Definition of storage space for profile obs.
135
136	IMPLICIT NONE
137
138	!! * Arguments
139	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
140	INTEGER, INTENT(IN) :: kt ! Time step
141	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
142	INTEGER, INTENT(IN) :: kpj
143	INTEGER, INTENT(IN) :: kpk
144	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
145	! (kit000-1 = restart time)
146	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
147	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
148	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
149	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
150	& ptn, & ! Model temperature field
151	& psn, & ! Model salinity field
152	& ptmask ! Land-sea mask
153	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
154	& pgdept ! Model array of depth levels
155	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
156	& kdailyavtypes! Types for daily averages
157	!! * Local declarations
158	INTEGER :: ji
159	INTEGER :: jj
160	INTEGER :: jk
161	INTEGER :: jobs
162	INTEGER :: inrc
163	INTEGER :: ipro
164	INTEGER :: idayend
165	INTEGER :: ista
166	INTEGER :: iend
167	INTEGER :: iobs
168	INTEGER, DIMENSION(imaxavtypes) :: &
169	& idailyavtypes
170	REAL(KIND=wp) :: zlam
171	REAL(KIND=wp) :: zphi
172	REAL(KIND=wp) :: zdaystp
173	REAL(KIND=wp), DIMENSION(kpk) :: &
174	& zobsmask, &
175	& zobsk, &
176	& zobs2k
177	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
178	& zweig
179	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
180	& zmask, &
181	& zintt, &
182	& zints, &
183	& zinmt, &
184	& zinms
185	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
186	& zglam, &
187	& zgphi
188	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
189	& igrdi, &
190	& igrdj
191
192	!------------------------------------------------------------------------
193	! Local initialization
194	!------------------------------------------------------------------------
195	! ... Record and data counters
196	inrc = kt - kit000 + 2
197	ipro = prodatqc%npstp(inrc)
198
199	! Daily average types
200	IF ( PRESENT(kdailyavtypes) ) THEN
201	idailyavtypes(:) = kdailyavtypes(:)
202	ELSE
203	idailyavtypes(:) = -1
204	ENDIF
205
206	! Initialize daily mean for first timestep
207	idayend = MOD( kt - kit000 + 1, kdaystp )
208
209	! Added kt == 0 test to catch restart case
210	IF ( idayend == 1 .OR. kt == 0) THEN
211	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
212	DO jk = 1, jpk
213	DO jj = 1, jpj
214	DO ji = 1, jpi
215	prodatqc%vdmean(ji,jj,jk,1) = 0.0
216	prodatqc%vdmean(ji,jj,jk,2) = 0.0
217	END DO
218	END DO
219	END DO
220	ENDIF
221
222	DO jk = 1, jpk
223	DO jj = 1, jpj
224	DO ji = 1, jpi
225	! Increment the temperature field for computing daily mean
226	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
227	& + ptn(ji,jj,jk)
228	! Increment the salinity field for computing daily mean
229	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
230	& + psn(ji,jj,jk)
231	END DO
232	END DO
233	END DO
234
235	! Compute the daily mean at the end of day
236	zdaystp = 1.0 / REAL( kdaystp )
237	IF ( idayend == 0 ) THEN
238	DO jk = 1, jpk
239	DO jj = 1, jpj
240	DO ji = 1, jpi
241	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
242	& * zdaystp
243	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
244	& * zdaystp
245	END DO
246	END DO
247	END DO
248	ENDIF
249
250	! Get the data for interpolation
251	ALLOCATE( &
252	& igrdi(2,2,ipro), &
253	& igrdj(2,2,ipro), &
254	& zglam(2,2,ipro), &
255	& zgphi(2,2,ipro), &
256	& zmask(2,2,kpk,ipro), &
257	& zintt(2,2,kpk,ipro), &
258	& zints(2,2,kpk,ipro) &
259	& )
260
261	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
262	iobs = jobs - prodatqc%nprofup
263	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
264	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
265	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
266	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
267	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
268	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
269	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
270	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
271	END DO
272
273	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
274	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
275	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
276	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
277	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
278
279	! At the end of the day also get interpolated means
280	IF ( idayend == 0 ) THEN
281
282	ALLOCATE( &
283	& zinmt(2,2,kpk,ipro), &
284	& zinms(2,2,kpk,ipro) &
285	& )
286
287	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
288	& prodatqc%vdmean(:,:,:,1), zinmt )
289	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
290	& prodatqc%vdmean(:,:,:,2), zinms )
291
292	ENDIF
293
294	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
295
296	iobs = jobs - prodatqc%nprofup
297
298	IF ( kt /= prodatqc%mstp(jobs) ) THEN
299
300	IF(lwp) THEN
301	WRITE(numout,*)
302	WRITE(numout,*) ' E R R O R : Observation', &
303	& ' time step is not consistent with the', &
304	& ' model time step'
305	WRITE(numout,*) ' ========='
306	WRITE(numout,*)
307	WRITE(numout,*) ' Record = ', jobs, &
308	& ' kt = ', kt, &
309	& ' mstp = ', prodatqc%mstp(jobs), &
310	& ' ntyp = ', prodatqc%ntyp(jobs)
311	ENDIF
312	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
313	ENDIF
314
315	zlam = prodatqc%rlam(jobs)
316	zphi = prodatqc%rphi(jobs)
317
318	! Horizontal weights and vertical mask
319
320	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
321	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
322
323	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
324	& zglam(:,:,iobs), zgphi(:,:,iobs), &
325	& zmask(:,:,:,iobs), zweig, zobsmask )
326
327	ENDIF
328
329	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
330
331	zobsk(:) = obfillflt
332
333	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
334
335	IF ( idayend == 0 ) THEN
336
337	! Daily averaged moored buoy (MRB) data
338
339	CALL obs_int_h2d( kpk, kpk, &
340	& zweig, zinmt(:,:,:,iobs), zobsk )
341
342
343	ELSE
344
345	CALL ctl_stop( ' A nonzero' // &
346	& ' number of profile T BUOY data should' // &
347	& ' only occur at the end of a given day' )
348
349	ENDIF
350
351	ELSE
352
353	! Point data
354
355	CALL obs_int_h2d( kpk, kpk, &
356	& zweig, zintt(:,:,:,iobs), zobsk )
357
358	ENDIF
359
360	!-------------------------------------------------------------
361	! Compute vertical second-derivative of the interpolating
362	! polynomial at obs points
363	!-------------------------------------------------------------
364
365	IF ( k1dint == 1 ) THEN
366	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
367	& pgdept, zobsmask )
368	ENDIF
369
370	!-----------------------------------------------------------------
371	! Vertical interpolation to the observation point
372	!-----------------------------------------------------------------
373	ista = prodatqc%npvsta(jobs,1)
374	iend = prodatqc%npvend(jobs,1)
375	CALL obs_int_z1d( kpk, &
376	& prodatqc%var(1)%mvk(ista:iend), &
377	& k1dint, iend - ista + 1, &
378	& prodatqc%var(1)%vdep(ista:iend), &
379	& zobsk, zobs2k, &
380	& prodatqc%var(1)%vmod(ista:iend), &
381	& pgdept, zobsmask )
382
383	ENDIF
384
385	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
386
387	zobsk(:) = obfillflt
388
389	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
390
391	IF ( idayend == 0 ) THEN
392
393	! Daily averaged moored buoy (MRB) data
394
395	CALL obs_int_h2d( kpk, kpk, &
396	& zweig, zinms(:,:,:,iobs), zobsk )
397
398	ELSE
399
400	CALL ctl_stop( ' A nonzero' // &
401	& ' number of profile S BUOY data should' // &
402	& ' only occur at the end of a given day' )
403
404	ENDIF
405
406	ELSE
407
408	! Point data
409
410	CALL obs_int_h2d( kpk, kpk, &
411	& zweig, zints(:,:,:,iobs), zobsk )
412
413	ENDIF
414
415
416	!-------------------------------------------------------------
417	! Compute vertical second-derivative of the interpolating
418	! polynomial at obs points
419	!-------------------------------------------------------------
420
421	IF ( k1dint == 1 ) THEN
422	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
423	& pgdept, zobsmask )
424	ENDIF
425
426	!----------------------------------------------------------------
427	! Vertical interpolation to the observation point
428	!----------------------------------------------------------------
429	ista = prodatqc%npvsta(jobs,2)
430	iend = prodatqc%npvend(jobs,2)
431	CALL obs_int_z1d( kpk, &
432	& prodatqc%var(2)%mvk(ista:iend),&
433	& k1dint, iend - ista + 1, &
434	& prodatqc%var(2)%vdep(ista:iend),&
435	& zobsk, zobs2k, &
436	& prodatqc%var(2)%vmod(ista:iend),&
437	& pgdept, zobsmask )
438
439	ENDIF
440
441	END DO
442
443	! Deallocate the data for interpolation
444	DEALLOCATE( &
445	& igrdi, &
446	& igrdj, &
447	& zglam, &
448	& zgphi, &
449	& zmask, &
450	& zintt, &
451	& zints &
452	& )
453	! At the end of the day also get interpolated means
454	IF ( idayend == 0 ) THEN
455	DEALLOCATE( &
456	& zinmt, &
457	& zinms &
458	& )
459	ENDIF
460
461	prodatqc%nprofup = prodatqc%nprofup + ipro
462
463	END SUBROUTINE obs_pro_opt
464
465	SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
466	& ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, &
467	& kdailyavtypes )
468	!!-----------------------------------------------------------------------
469	!!
470	!! * ROUTINE obs_pro_opt *
471	!!
472	!! ** Purpose : Compute the model counterpart of profiles
473	!! data by interpolating from the model grid to the
474	!! observation point. Generalised vertical coordinate version
475	!!
476	!! ** Method : Linearly interpolate to each observation point using
477	!! the model values at the corners of the surrounding grid box.
478	!!
479	!! First, model values on the model grid are interpolated vertically to the
480	!! Depths of the profile observations. Two vertical interpolation schemes are
481	!! available:
482	!! - linear (k1dint = 0)
483	!! - Cubic spline (k1dint = 1)
484	!!
485	!!
486	!! Secondly the interpolated values are interpolated horizontally to the
487	!! obs (lon, lat) point.
488	!! Several horizontal interpolation schemes are available:
489	!! - distance-weighted (great circle) (k2dint = 0)
490	!! - distance-weighted (small angle) (k2dint = 1)
491	!! - bilinear (geographical grid) (k2dint = 2)
492	!! - bilinear (quadrilateral grid) (k2dint = 3)
493	!! - polynomial (quadrilateral grid) (k2dint = 4)
494	!!
495	!! For the cubic spline the 2nd derivative of the interpolating
496	!! polynomial is computed before entering the vertical interpolation
497	!! routine.
498	!!
499	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
500	!! a daily mean model temperature field. So, we first compute
501	!! the mean, then interpolate only at the end of the day.
502	!!
503	!! This is the procedure to be used with generalised vertical model
504	!! coordinates (ie s-coordinates. It is ~4x slower than the equivalent
505	!! horizontal then vertical interpolation algorithm, but can deal with situations
506	!! where the model levels are not flat.
507	!! ONLY PERFORMED if ln_sco=.TRUE.
508	!!
509	!! Note: the in situ temperature observations must be converted
510	!! to potential temperature (the model variable) prior to
511	!! assimilation.
512	!!??????????????????????????????????????????????????????????????
513	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
514	!!??????????????????????????????????????????????????????????????
515	!!
516	!! ** Action :
517	!!
518	!! History :
519	!! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised
520	!! vertical coordinates
521	!!-----------------------------------------------------------------------
522
523	!! * Modules used
524	USE obs_profiles_def ! Definition of storage space for profile obs.
525	USE dom_oce, ONLY : &
526	#if defined key_vvl
527	& gdepw_n
528	#else
529	& gdepw_0
530	#endif
531
532	IMPLICIT NONE
533
534	!! * Arguments
535	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
536	INTEGER, INTENT(IN) :: kt ! Time step
537	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
538	INTEGER, INTENT(IN) :: kpj
539	INTEGER, INTENT(IN) :: kpk
540	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
541	! (kit000-1 = restart time)
542	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
543	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
544	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
545	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
546	& ptn, & ! Model temperature field
547	& psn, & ! Model salinity field
548	& ptmask ! Land-sea mask
549	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,jpj,kpk) :: &
550	& pgdept ! Model array of depth T levels
551	& pgdepw ! Model array of depth W levels
552	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
553	& kdailyavtypes ! Types for daily averages
554
555	!! * Local declarations
556	INTEGER :: ji
557	INTEGER :: jj
558	INTEGER :: jk
559	INTEGER :: iico, ijco
560	INTEGER :: jobs
561	INTEGER :: inrc
562	INTEGER :: ipro
563	INTEGER :: idayend
564	INTEGER :: ista
565	INTEGER :: iend
566	INTEGER :: iobs
567	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
568	INTEGER, DIMENSION(imaxavtypes) :: &
569	& idailyavtypes
570	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
571	& igrdi, &
572	& igrdj
573	INTEGER :: &
574	& inum_obs
575	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
576	REAL(KIND=wp) :: zlam
577	REAL(KIND=wp) :: zphi
578	REAL(KIND=wp) :: zdaystp
579	REAL(KIND=wp), DIMENSION(kpk) :: &
580	& zobsmask, &
581	& zobsk, &
582	& zobs2k
583	REAL(KIND=wp), DIMENSION(2,2,1) :: &
584	& zweig, &
585	& l_zweig
586	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
587	& zmask, &
588	& zintt, &
589	& zints, &
590	& zinmt, &
591	& zgdept,&
592	& zgdepw,&
593	& zinms
594	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
595	& zglam, &
596	& zgphi
597	REAL(KIND=wp), DIMENSION(1) :: zmsk_1
598	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
599
600	!------------------------------------------------------------------------
601	! Local initialization
602	!------------------------------------------------------------------------
603	! ... Record and data counters
604	inrc = kt - kit000 + 2
605	ipro = prodatqc%npstp(inrc)
606
607	! Daily average types
608	IF ( PRESENT(kdailyavtypes) ) THEN
609	idailyavtypes(:) = kdailyavtypes(:)
610	ELSE
611	idailyavtypes(:) = -1
612	ENDIF
613
614	! Initialize daily mean for first time-step
615	idayend = MOD( kt - kit000 + 1, kdaystp )
616
617	! Added kt == 0 test to catch restart case
618	IF ( idayend == 1 .OR. kt == 0) THEN
619
620	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
621	DO jk = 1, jpk
622	DO jj = 1, jpj
623	DO ji = 1, jpi
624	prodatqc%vdmean(ji,jj,jk,1) = 0.0
625	prodatqc%vdmean(ji,jj,jk,2) = 0.0
626	END DO
627	END DO
628	END DO
629
630	ENDIF
631
632	DO jk = 1, jpk
633	DO jj = 1, jpj
634	DO ji = 1, jpi
635	! Increment the temperature field for computing daily mean
636	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
637	& + ptn(ji,jj,jk)
638	! Increment the salinity field for computing daily mean
639	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
640	& + psn(ji,jj,jk)
641	END DO
642	END DO
643	END DO
644
645	! Compute the daily mean at the end of day
646	zdaystp = 1.0 / REAL( kdaystp )
647	IF ( idayend == 0 ) THEN
648	DO jk = 1, jpk
649	DO jj = 1, jpj
650	DO ji = 1, jpi
651	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
652	& * zdaystp
653	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
654	& * zdaystp
655	END DO
656	END DO
657	END DO
658	ENDIF
659
660	! Get the data for interpolation
661	ALLOCATE( &
662	& igrdi(2,2,ipro), &
663	& igrdj(2,2,ipro), &
664	& zglam(2,2,ipro), &
665	& zgphi(2,2,ipro), &
666	& zmask(2,2,kpk,ipro), &
667	& zintt(2,2,kpk,ipro), &
668	& zints(2,2,kpk,ipro), &
669	& zgdept(2,2,kpk,ipro), &
670	& zgdepw(2,2,kpk,ipro) &
671	& )
672
673	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
674	iobs = jobs - prodatqc%nprofup
675	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
676	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
677	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
678	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
679	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
680	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
681	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
682	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
683	END DO
684
685	! Initialise depth arrays
686	zgdept = 0.0
687	zgdepw = 0.0
688
689	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
690	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
691	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
692	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
693	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
694	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdept(:,:,:), &
695	& zgdept )
696	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdepw(:,:,:), &
697	& zgdepw )
698
699	! At the end of the day also get interpolated means
700	IF ( idayend == 0 ) THEN
701
702	ALLOCATE( &
703	& zinmt(2,2,kpk,ipro), &
704	& zinms(2,2,kpk,ipro) &
705	& )
706
707	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
708	& prodatqc%vdmean(:,:,:,1), zinmt )
709	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
710	& prodatqc%vdmean(:,:,:,2), zinms )
711
712	ENDIF
713
714	! Return if no observations to process
715	! Has to be done after comm commands to ensure processors
716	! stay in sync
717	IF ( ipro == 0 ) RETURN
718
719	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
720
721	iobs = jobs - prodatqc%nprofup
722
723	IF ( kt /= prodatqc%mstp(jobs) ) THEN
724
725	IF(lwp) THEN
726	WRITE(numout,*)
727	WRITE(numout,*) ' E R R O R : Observation', &
728	& ' time step is not consistent with the', &
729	& ' model time step'
730	WRITE(numout,*) ' ========='
731	WRITE(numout,*)
732	WRITE(numout,*) ' Record = ', jobs, &
733	& ' kt = ', kt, &
734	& ' mstp = ', prodatqc%mstp(jobs), &
735	& ' ntyp = ', prodatqc%ntyp(jobs)
736	ENDIF
737	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
738	ENDIF
739
740	zlam = prodatqc%rlam(jobs)
741	zphi = prodatqc%rphi(jobs)
742
743	! Horizontal weights
744	! Only calculated once, for both T and S.
745	! Masked values are calculated later.
746
747	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
748	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
749
750	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
751	& zglam(:,:,iobs), zgphi(:,:,iobs), &
752	& zmask(:,:,1,iobs), zweig, zmsk_1 )
753
754	ENDIF
755
756	! IF zmsk_1 = 0; then ob is on land
757	IF (zmsk_1(1) < 0.1) THEN
758	WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'
759
760	ELSE
761
762	! Temperature
763
764	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
765
766	zobsk(:) = obfillflt
767
768	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
769
770	IF ( idayend == 0 ) THEN
771
772	! Daily averaged moored buoy (MRB) data
773
774	! vertically interpolate all 4 corners
775	ista = prodatqc%npvsta(jobs,1)
776	iend = prodatqc%npvend(jobs,1)
777	inum_obs = iend - ista + 1
778	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
779
780	DO iin=1,2
781	DO ijn=1,2
782
783
784
785	IF ( k1dint == 1 ) THEN
786	CALL obs_int_z1d_spl( kpk, &
787	& zinmt(iin,ijn,:,iobs), &
788	& zobs2k, zgdept(iin,ijn,:,iobs), &
789	& zmask(iin,ijn,:,iobs))
790	ENDIF
791
792	CALL obs_level_search(kpk, &
793	& zgdept(iin,ijn,:,iobs), &
794	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
795	& iv_indic)
796	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
797	& prodatqc%var(1)%vdep(ista:iend), &
798	& zinmt(iin,ijn,:,iobs), &
799	& zobs2k, interp_corner(iin,ijn,:), &
800	& zgdept(iin,ijn,:,iobs), &
801	& zmask(iin,ijn,:,iobs))
802
803	ENDDO
804	ENDDO
805
806
807	ELSE
808
809	CALL ctl_stop( ' A nonzero' // &
810	& ' number of profile T BUOY data should' // &
811	& ' only occur at the end of a given day' )
812
813	ENDIF
814
815	ELSE
816
817	! Point data
818
819	! vertically interpolate all 4 corners
820	ista = prodatqc%npvsta(jobs,1)
821	iend = prodatqc%npvend(jobs,1)
822	inum_obs = iend - ista + 1
823	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
824	DO iin=1,2
825	DO ijn=1,2
826
827
828	IF ( k1dint == 1 ) THEN
829	CALL obs_int_z1d_spl( kpk, &
830	& zintt(iin,ijn,:,iobs),&
831	& zobs2k, zgdept(iin,ijn,:,iobs), &
832	& zmask(iin,ijn,:,iobs))
833
834	ENDIF
835
836	CALL obs_level_search(kpk, &
837	& zgdept(iin,ijn,:,iobs),&
838	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
839	& iv_indic)
840	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
841	& prodatqc%var(1)%vdep(ista:iend), &
842	& zintt(iin,ijn,:,iobs), &
843	& zobs2k,interp_corner(iin,ijn,:), &
844	& zgdept(iin,ijn,:,iobs), &
845	& zmask(iin,ijn,:,iobs) )
846
847	ENDDO
848	ENDDO
849
850	ENDIF
851
852	!-------------------------------------------------------------
853	! Compute the horizontal interpolation for every profile level
854	!-------------------------------------------------------------
855
856	DO ikn=1,inum_obs
857	iend=ista+ikn-1
858
859	l_zweig(:,:,1) = 0._wp
860
861	! This code forces the horizontal weights to be
862	! zero IF the observation is below the bottom of the
863	! corners of the interpolation nodes, Or if it is in
864	! the mask. This is important for observations are near
865	! steep bathymetry
866	DO iin=1,2
867	DO ijn=1,2
868
869	depth_loop1: DO ik=kpk,2,-1
870	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
871
872	l_zweig(iin,ijn,1) = &
873	& zweig(iin,ijn,1) * &
874	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
875	& - prodatqc%var(1)%vdep(iend)),0._wp)
876
877	EXIT depth_loop1
878	ENDIF
879	ENDDO depth_loop1
880
881	ENDDO
882	ENDDO
883
884	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
885	& prodatqc%var(1)%vmod(iend:iend) )
886
887	ENDDO
888
889
890	DEALLOCATE(interp_corner,iv_indic)
891
892	ENDIF
893
894
895	! Salinity
896
897	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
898
899	zobsk(:) = obfillflt
900
901	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
902
903	IF ( idayend == 0 ) THEN
904
905	! Daily averaged moored buoy (MRB) data
906
907	! vertically interpolate all 4 corners
908	ista = prodatqc%npvsta(iobs,2)
909	iend = prodatqc%npvend(iobs,2)
910	inum_obs = iend - ista + 1
911	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
912
913	DO iin=1,2
914	DO ijn=1,2
915
916
917
918	IF ( k1dint == 1 ) THEN
919	CALL obs_int_z1d_spl( kpk, &
920	& zinms(iin,ijn,:,iobs), &
921	& zobs2k, zgdept(iin,ijn,:,iobs), &
922	& zmask(iin,ijn,:,iobs))
923	ENDIF
924
925	CALL obs_level_search(kpk, &
926	& zgdept(iin,ijn,:,iobs), &
927	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
928	& iv_indic)
929	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
930	& prodatqc%var(2)%vdep(ista:iend), &
931	& zinms(iin,ijn,:,iobs), &
932	& zobs2k, interp_corner(iin,ijn,:), &
933	& zgdept(iin,ijn,:,iobs), &
934	& zmask(iin,ijn,:,iobs))
935
936	ENDDO
937	ENDDO
938
939
940	ELSE
941
942	CALL ctl_stop( ' A nonzero' // &
943	& ' number of profile T BUOY data should' // &
944	& ' only occur at the end of a given day' )
945
946	ENDIF
947
948	ELSE
949
950	! Point data
951
952	! vertically interpolate all 4 corners
953	ista = prodatqc%npvsta(jobs,2)
954	iend = prodatqc%npvend(jobs,2)
955	inum_obs = iend - ista + 1
956	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
957
958	DO iin=1,2
959	DO ijn=1,2
960
961
962	IF ( k1dint == 1 ) THEN
963	CALL obs_int_z1d_spl( kpk, &
964	& zints(iin,ijn,:,iobs),&
965	& zobs2k, zgdept(iin,ijn,:,iobs), &
966	& zmask(iin,ijn,:,iobs))
967
968	ENDIF
969
970	CALL obs_level_search(kpk, &
971	& zgdept(iin,ijn,:,iobs),&
972	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
973	& iv_indic)
974	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
975	& prodatqc%var(2)%vdep(ista:iend), &
976	& zints(iin,ijn,:,iobs), &
977	& zobs2k,interp_corner(iin,ijn,:), &
978	& zgdept(iin,ijn,:,iobs), &
979	& zmask(iin,ijn,:,iobs) )
980
981	ENDDO
982	ENDDO
983
984	ENDIF
985
986	!-------------------------------------------------------------
987	! Compute the horizontal interpolation for every profile level
988	!-------------------------------------------------------------
989
990	DO ikn=1,inum_obs
991	iend=ista+ikn-1
992
993	l_zweig(:,:,1) = 0._wp
994
995	! This code forces the horizontal weights to be
996	! zero IF the observation is below the bottom of the
997	! corners of the interpolation nodes, Or if it is in
998	! the mask. This is important for observations are near
999	! steep bathymetry
1000	DO iin=1,2
1001	DO ijn=1,2
1002
1003	depth_loop2: DO ik=kpk,2,-1
1004	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
1005
1006	l_zweig(iin,ijn,1) = &
1007	& zweig(iin,ijn,1) * &
1008	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
1009	& - prodatqc%var(2)%vdep(iend)),0._wp)
1010
1011	EXIT depth_loop2
1012	ENDIF
1013	ENDDO depth_loop2
1014
1015	ENDDO
1016	ENDDO
1017
1018	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
1019	& prodatqc%var(2)%vmod(iend:iend) )
1020
1021	ENDDO
1022
1023
1024	DEALLOCATE(interp_corner,iv_indic)
1025
1026	ENDIF
1027
1028	ENDIF
1029
1030	END DO
1031
1032	! Deallocate the data for interpolation
1033	DEALLOCATE( &
1034	& igrdi, &
1035	& igrdj, &
1036	& zglam, &
1037	& zgphi, &
1038	& zmask, &
1039	& zintt, &
1040	& zints, &
1041	& zgdept,&
1042	& zgdepw &
1043	& )
1044	! At the end of the day also get interpolated means
1045	IF ( idayend == 0 ) THEN
1046	DEALLOCATE( &
1047	& zinmt, &
1048	& zinms &
1049	& )
1050	ENDIF
1051
1052	prodatqc%nprofup = prodatqc%nprofup + ipro
1053
1054	END SUBROUTINE obs_pro_sco_opt
1055
1056	SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
1057	& psshn, psshmask, k2dint )
1058	!!-----------------------------------------------------------------------
1059	!!
1060	!! * ROUTINE obs_sla_opt *
1061	!!
1062	!! ** Purpose : Compute the model counterpart of sea level anomaly
1063	!! data by interpolating from the model grid to the
1064	!! observation point.
1065	!!
1066	!! ** Method : Linearly interpolate to each observation point using
1067	!! the model values at the corners of the surrounding grid box.
1068	!!
1069	!! The now model SSH is first computed at the obs (lon, lat) point.
1070	!!
1071	!! Several horizontal interpolation schemes are available:
1072	!! - distance-weighted (great circle) (k2dint = 0)
1073	!! - distance-weighted (small angle) (k2dint = 1)
1074	!! - bilinear (geographical grid) (k2dint = 2)
1075	!! - bilinear (quadrilateral grid) (k2dint = 3)
1076	!! - polynomial (quadrilateral grid) (k2dint = 4)
1077	!!
1078	!! The sea level anomaly at the observation points is then computed
1079	!! by removing a mean dynamic topography (defined at the obs. point).
1080	!!
1081	!! ** Action :
1082	!!
1083	!! History :
1084	!! ! 07-03 (A. Weaver)
1085	!!-----------------------------------------------------------------------
1086
1087	!! * Modules used
1088	USE obs_surf_def ! Definition of storage space for surface observations
1089
1090	IMPLICIT NONE
1091
1092	!! * Arguments
1093	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening
1094	INTEGER, INTENT(IN) :: kt ! Time step
1095	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1096	INTEGER, INTENT(IN) :: kpj
1097	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1098	! (kit000-1 = restart time)
1099	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1100	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1101	& psshn, & ! Model SSH field
1102	& psshmask ! Land-sea mask
1103
1104	!! * Local declarations
1105	INTEGER :: ji
1106	INTEGER :: jj
1107	INTEGER :: jobs
1108	INTEGER :: inrc
1109	INTEGER :: isla
1110	INTEGER :: iobs
1111	REAL(KIND=wp) :: zlam
1112	REAL(KIND=wp) :: zphi
1113	REAL(KIND=wp) :: zext(1), zobsmask(1)
1114	REAL(kind=wp), DIMENSION(2,2,1) :: &
1115	& zweig
1116	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1117	& zmask, &
1118	& zsshl, &
1119	& zglam, &
1120	& zgphi
1121	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1122	& igrdi, &
1123	& igrdj
1124
1125	!------------------------------------------------------------------------
1126	! Local initialization
1127	!------------------------------------------------------------------------
1128	! ... Record and data counters
1129	inrc = kt - kit000 + 2
1130	isla = sladatqc%nsstp(inrc)
1131
1132	! Get the data for interpolation
1133
1134	ALLOCATE( &
1135	& igrdi(2,2,isla), &
1136	& igrdj(2,2,isla), &
1137	& zglam(2,2,isla), &
1138	& zgphi(2,2,isla), &
1139	& zmask(2,2,isla), &
1140	& zsshl(2,2,isla) &
1141	& )
1142
1143	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1144	iobs = jobs - sladatqc%nsurfup
1145	igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
1146	igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
1147	igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
1148	igrdj(1,2,iobs) = sladatqc%mj(jobs)
1149	igrdi(2,1,iobs) = sladatqc%mi(jobs)
1150	igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
1151	igrdi(2,2,iobs) = sladatqc%mi(jobs)
1152	igrdj(2,2,iobs) = sladatqc%mj(jobs)
1153	END DO
1154
1155	CALL obs_int_comm_2d( 2, 2, isla, &
1156	& igrdi, igrdj, glamt, zglam )
1157	CALL obs_int_comm_2d( 2, 2, isla, &
1158	& igrdi, igrdj, gphit, zgphi )
1159	CALL obs_int_comm_2d( 2, 2, isla, &
1160	& igrdi, igrdj, psshmask, zmask )
1161	CALL obs_int_comm_2d( 2, 2, isla, &
1162	& igrdi, igrdj, psshn, zsshl )
1163
1164	! Loop over observations
1165
1166	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1167
1168	iobs = jobs - sladatqc%nsurfup
1169
1170	IF ( kt /= sladatqc%mstp(jobs) ) THEN
1171
1172	IF(lwp) THEN
1173	WRITE(numout,*)
1174	WRITE(numout,*) ' E R R O R : Observation', &
1175	& ' time step is not consistent with the', &
1176	& ' model time step'
1177	WRITE(numout,*) ' ========='
1178	WRITE(numout,*)
1179	WRITE(numout,*) ' Record = ', jobs, &
1180	& ' kt = ', kt, &
1181	& ' mstp = ', sladatqc%mstp(jobs), &
1182	& ' ntyp = ', sladatqc%ntyp(jobs)
1183	ENDIF
1184	CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
1185
1186	ENDIF
1187
1188	zlam = sladatqc%rlam(jobs)
1189	zphi = sladatqc%rphi(jobs)
1190
1191	! Get weights to interpolate the model SSH to the observation point
1192	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1193	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1194	& zmask(:,:,iobs), zweig, zobsmask )
1195
1196
1197	! Interpolate the model SSH to the observation point
1198	CALL obs_int_h2d( 1, 1, &
1199	& zweig, zsshl(:,:,iobs), zext )
1200
1201	sladatqc%rext(jobs,1) = zext(1)
1202	! ... Remove the MDT at the observation point
1203	sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
1204
1205	END DO
1206
1207	! Deallocate the data for interpolation
1208	DEALLOCATE( &
1209	& igrdi, &
1210	& igrdj, &
1211	& zglam, &
1212	& zgphi, &
1213	& zmask, &
1214	& zsshl &
1215	& )
1216
1217	sladatqc%nsurfup = sladatqc%nsurfup + isla
1218
1219	END SUBROUTINE obs_sla_opt
1220
1221	SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
1222	& psstn, psstmask, k2dint, ld_nightav )
1223	!!-----------------------------------------------------------------------
1224	!!
1225	!! * ROUTINE obs_sst_opt *
1226	!!
1227	!! ** Purpose : Compute the model counterpart of surface temperature
1228	!! data by interpolating from the model grid to the
1229	!! observation point.
1230	!!
1231	!! ** Method : Linearly interpolate to each observation point using
1232	!! the model values at the corners of the surrounding grid box.
1233	!!
1234	!! The now model SST is first computed at the obs (lon, lat) point.
1235	!!
1236	!! Several horizontal interpolation schemes are available:
1237	!! - distance-weighted (great circle) (k2dint = 0)
1238	!! - distance-weighted (small angle) (k2dint = 1)
1239	!! - bilinear (geographical grid) (k2dint = 2)
1240	!! - bilinear (quadrilateral grid) (k2dint = 3)
1241	!! - polynomial (quadrilateral grid) (k2dint = 4)
1242	!!
1243	!!
1244	!! ** Action :
1245	!!
1246	!! History :
1247	!! ! 07-07 (S. Ricci ) : Original
1248	!!
1249	!!-----------------------------------------------------------------------
1250
1251	!! * Modules used
1252	USE obs_surf_def ! Definition of storage space for surface observations
1253	USE sbcdcy
1254
1255	IMPLICIT NONE
1256
1257	!! * Arguments
1258	TYPE(obs_surf), INTENT(INOUT) :: &
1259	& sstdatqc ! Subset of surface data not failing screening
1260	INTEGER, INTENT(IN) :: kt ! Time step
1261	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1262	INTEGER, INTENT(IN) :: kpj
1263	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1264	! (kit000-1 = restart time)
1265	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1266	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1267	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1268	& psstn, & ! Model SST field
1269	& psstmask ! Land-sea mask
1270
1271	!! * Local declarations
1272	INTEGER :: ji
1273	INTEGER :: jj
1274	INTEGER :: jobs
1275	INTEGER :: inrc
1276	INTEGER :: isst
1277	INTEGER :: iobs
1278	INTEGER :: idayend
1279	REAL(KIND=wp) :: zlam
1280	REAL(KIND=wp) :: zphi
1281	REAL(KIND=wp) :: zext(1), zobsmask(1)
1282	REAL(KIND=wp) :: zdaystp
1283	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1284	& icount_sstnight, &
1285	& imask_night
1286	REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1287	& zintmp, &
1288	& zouttmp, &
1289	& zmeanday ! to compute model sst in region of 24h daylight (pole)
1290	REAL(kind=wp), DIMENSION(2,2,1) :: &
1291	& zweig
1292	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1293	& zmask, &
1294	& zsstl, &
1295	& zsstm, &
1296	& zglam, &
1297	& zgphi
1298	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1299	& igrdi, &
1300	& igrdj
1301	LOGICAL, INTENT(IN) :: ld_nightav
1302
1303	!-----------------------------------------------------------------------
1304	! Local initialization
1305	!-----------------------------------------------------------------------
1306	! ... Record and data counters
1307	inrc = kt - kit000 + 2
1308	isst = sstdatqc%nsstp(inrc)
1309
1310	IF ( ld_nightav ) THEN
1311
1312	! Initialize array for night mean
1313
1314	IF ( kt .EQ. 0 ) THEN
1315	ALLOCATE ( icount_sstnight(kpi,kpj) )
1316	ALLOCATE ( imask_night(kpi,kpj) )
1317	ALLOCATE ( zintmp(kpi,kpj) )
1318	ALLOCATE ( zouttmp(kpi,kpj) )
1319	ALLOCATE ( zmeanday(kpi,kpj) )
1320	nday_qsr = -1 ! initialisation flag for nbc_dcy
1321	ENDIF
1322
1323	! Initialize daily mean for first timestep
1324	idayend = MOD( kt - kit000 + 1, kdaystp )
1325
1326	! Added kt == 0 test to catch restart case
1327	IF ( idayend == 1 .OR. kt == 0) THEN
1328	IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
1329	DO jj = 1, jpj
1330	DO ji = 1, jpi
1331	sstdatqc%vdmean(ji,jj) = 0.0
1332	zmeanday(ji,jj) = 0.0
1333	icount_sstnight(ji,jj) = 0
1334	END DO
1335	END DO
1336	ENDIF
1337
1338	zintmp(:,:) = 0.0
1339	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1340	imask_night(:,:) = INT( zouttmp(:,:) )
1341
1342	DO jj = 1, jpj
1343	DO ji = 1, jpi
1344	! Increment the temperature field for computing night mean and counter
1345	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1346	& + psstn(ji,jj)*imask_night(ji,jj)
1347	zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj)
1348	icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
1349	END DO
1350	END DO
1351
1352	! Compute the daily mean at the end of day
1353
1354	zdaystp = 1.0 / REAL( kdaystp )
1355
1356	IF ( idayend == 0 ) THEN
1357	DO jj = 1, jpj
1358	DO ji = 1, jpi
1359	! Test if "no night" point
1360	IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
1361	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1362	& / icount_sstnight(ji,jj)
1363	ELSE
1364	sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1365	ENDIF
1366	END DO
1367	END DO
1368	ENDIF
1369
1370	ENDIF
1371
1372	! Get the data for interpolation
1373
1374	ALLOCATE( &
1375	& igrdi(2,2,isst), &
1376	& igrdj(2,2,isst), &
1377	& zglam(2,2,isst), &
1378	& zgphi(2,2,isst), &
1379	& zmask(2,2,isst), &
1380	& zsstl(2,2,isst) &
1381	& )
1382
1383	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1384	iobs = jobs - sstdatqc%nsurfup
1385	igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
1386	igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
1387	igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
1388	igrdj(1,2,iobs) = sstdatqc%mj(jobs)
1389	igrdi(2,1,iobs) = sstdatqc%mi(jobs)
1390	igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
1391	igrdi(2,2,iobs) = sstdatqc%mi(jobs)
1392	igrdj(2,2,iobs) = sstdatqc%mj(jobs)
1393	END DO
1394
1395	CALL obs_int_comm_2d( 2, 2, isst, &
1396	& igrdi, igrdj, glamt, zglam )
1397	CALL obs_int_comm_2d( 2, 2, isst, &
1398	& igrdi, igrdj, gphit, zgphi )
1399	CALL obs_int_comm_2d( 2, 2, isst, &
1400	& igrdi, igrdj, psstmask, zmask )
1401	CALL obs_int_comm_2d( 2, 2, isst, &
1402	& igrdi, igrdj, psstn, zsstl )
1403
1404	! At the end of the day get interpolated means
1405	IF ( idayend == 0 .AND. ld_nightav ) THEN
1406
1407	ALLOCATE( &
1408	& zsstm(2,2,isst) &
1409	& )
1410
1411	CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
1412	& sstdatqc%vdmean(:,:), zsstm )
1413
1414	ENDIF
1415
1416	! Loop over observations
1417
1418	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1419
1420	iobs = jobs - sstdatqc%nsurfup
1421
1422	IF ( kt /= sstdatqc%mstp(jobs) ) THEN
1423
1424	IF(lwp) THEN
1425	WRITE(numout,*)
1426	WRITE(numout,*) ' E R R O R : Observation', &
1427	& ' time step is not consistent with the', &
1428	& ' model time step'
1429	WRITE(numout,*) ' ========='
1430	WRITE(numout,*)
1431	WRITE(numout,*) ' Record = ', jobs, &
1432	& ' kt = ', kt, &
1433	& ' mstp = ', sstdatqc%mstp(jobs), &
1434	& ' ntyp = ', sstdatqc%ntyp(jobs)
1435	ENDIF
1436	CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
1437
1438	ENDIF
1439
1440	zlam = sstdatqc%rlam(jobs)
1441	zphi = sstdatqc%rphi(jobs)
1442
1443	! Get weights to interpolate the model SST to the observation point
1444	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1445	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1446	& zmask(:,:,iobs), zweig, zobsmask )
1447
1448	! Interpolate the model SST to the observation point
1449
1450	IF ( ld_nightav ) THEN
1451
1452	IF ( idayend == 0 ) THEN
1453	! Daily averaged/diurnal cycle of SST data
1454	CALL obs_int_h2d( 1, 1, &
1455	& zweig, zsstm(:,:,iobs), zext )
1456	ELSE
1457	CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // &
1458	& ' number of night SST data should' // &
1459	& ' only occur at the end of a given day' )
1460	ENDIF
1461
1462	ELSE
1463
1464	CALL obs_int_h2d( 1, 1, &
1465	& zweig, zsstl(:,:,iobs), zext )
1466
1467	ENDIF
1468	sstdatqc%rmod(jobs,1) = zext(1)
1469
1470	END DO
1471
1472	! Deallocate the data for interpolation
1473	DEALLOCATE( &
1474	& igrdi, &
1475	& igrdj, &
1476	& zglam, &
1477	& zgphi, &
1478	& zmask, &
1479	& zsstl &
1480	& )
1481
1482	! At the end of the day also get interpolated means
1483	IF ( idayend == 0 .AND. ld_nightav ) THEN
1484	DEALLOCATE( &
1485	& zsstm &
1486	& )
1487	ENDIF
1488
1489	sstdatqc%nsurfup = sstdatqc%nsurfup + isst
1490
1491	END SUBROUTINE obs_sst_opt
1492
1493	SUBROUTINE obs_sss_opt
1494	!!-----------------------------------------------------------------------
1495	!!
1496	!! * ROUTINE obs_sss_opt *
1497	!!
1498	!! ** Purpose : Compute the model counterpart of sea surface salinity
1499	!! data by interpolating from the model grid to the
1500	!! observation point.
1501	!!
1502	!! ** Method :
1503	!!
1504	!! ** Action :
1505	!!
1506	!! History :
1507	!! ! ??-??
1508	!!-----------------------------------------------------------------------
1509
1510	IMPLICIT NONE
1511
1512	END SUBROUTINE obs_sss_opt
1513
1514	SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
1515	& pseaicen, pseaicemask, k2dint )
1516
1517	!!-----------------------------------------------------------------------
1518	!!
1519	!! * ROUTINE obs_seaice_opt *
1520	!!
1521	!! ** Purpose : Compute the model counterpart of surface temperature
1522	!! data by interpolating from the model grid to the
1523	!! observation point.
1524	!!
1525	!! ** Method : Linearly interpolate to each observation point using
1526	!! the model values at the corners of the surrounding grid box.
1527	!!
1528	!! The now model sea ice is first computed at the obs (lon, lat) point.
1529	!!
1530	!! Several horizontal interpolation schemes are available:
1531	!! - distance-weighted (great circle) (k2dint = 0)
1532	!! - distance-weighted (small angle) (k2dint = 1)
1533	!! - bilinear (geographical grid) (k2dint = 2)
1534	!! - bilinear (quadrilateral grid) (k2dint = 3)
1535	!! - polynomial (quadrilateral grid) (k2dint = 4)
1536	!!
1537	!!
1538	!! ** Action :
1539	!!
1540	!! History :
1541	!! ! 07-07 (S. Ricci ) : Original
1542	!!
1543	!!-----------------------------------------------------------------------
1544
1545	!! * Modules used
1546	USE obs_surf_def ! Definition of storage space for surface observations
1547
1548	IMPLICIT NONE
1549
1550	!! * Arguments
1551	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening
1552	INTEGER, INTENT(IN) :: kt ! Time step
1553	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1554	INTEGER, INTENT(IN) :: kpj
1555	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1556	! (kit000-1 = restart time)
1557	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1558	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1559	& pseaicen, & ! Model sea ice field
1560	& pseaicemask ! Land-sea mask
1561
1562	!! * Local declarations
1563	INTEGER :: ji
1564	INTEGER :: jj
1565	INTEGER :: jobs
1566	INTEGER :: inrc
1567	INTEGER :: iseaice
1568	INTEGER :: iobs
1569
1570	REAL(KIND=wp) :: zlam
1571	REAL(KIND=wp) :: zphi
1572	REAL(KIND=wp) :: zext(1), zobsmask(1)
1573	REAL(kind=wp), DIMENSION(2,2,1) :: &
1574	& zweig
1575	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1576	& zmask, &
1577	& zseaicel, &
1578	& zglam, &
1579	& zgphi
1580	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1581	& igrdi, &
1582	& igrdj
1583
1584	!------------------------------------------------------------------------
1585	! Local initialization
1586	!------------------------------------------------------------------------
1587	! ... Record and data counters
1588	inrc = kt - kit000 + 2
1589	iseaice = seaicedatqc%nsstp(inrc)
1590
1591	! Get the data for interpolation
1592
1593	ALLOCATE( &
1594	& igrdi(2,2,iseaice), &
1595	& igrdj(2,2,iseaice), &
1596	& zglam(2,2,iseaice), &
1597	& zgphi(2,2,iseaice), &
1598	& zmask(2,2,iseaice), &
1599	& zseaicel(2,2,iseaice) &
1600	& )
1601
1602	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1603	iobs = jobs - seaicedatqc%nsurfup
1604	igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1605	igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1606	igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1607	igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1608	igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1609	igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1610	igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1611	igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1612	END DO
1613
1614	CALL obs_int_comm_2d( 2, 2, iseaice, &
1615	& igrdi, igrdj, glamt, zglam )
1616	CALL obs_int_comm_2d( 2, 2, iseaice, &
1617	& igrdi, igrdj, gphit, zgphi )
1618	CALL obs_int_comm_2d( 2, 2, iseaice, &
1619	& igrdi, igrdj, pseaicemask, zmask )
1620	CALL obs_int_comm_2d( 2, 2, iseaice, &
1621	& igrdi, igrdj, pseaicen, zseaicel )
1622
1623	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1624
1625	iobs = jobs - seaicedatqc%nsurfup
1626
1627	IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1628
1629	IF(lwp) THEN
1630	WRITE(numout,*)
1631	WRITE(numout,*) ' E R R O R : Observation', &
1632	& ' time step is not consistent with the', &
1633	& ' model time step'
1634	WRITE(numout,*) ' ========='
1635	WRITE(numout,*)
1636	WRITE(numout,*) ' Record = ', jobs, &
1637	& ' kt = ', kt, &
1638	& ' mstp = ', seaicedatqc%mstp(jobs), &
1639	& ' ntyp = ', seaicedatqc%ntyp(jobs)
1640	ENDIF
1641	CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1642
1643	ENDIF
1644
1645	zlam = seaicedatqc%rlam(jobs)
1646	zphi = seaicedatqc%rphi(jobs)
1647
1648	! Get weights to interpolate the model sea ice to the observation point
1649	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1650	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1651	& zmask(:,:,iobs), zweig, zobsmask )
1652
1653	! ... Interpolate the model sea ice to the observation point
1654	CALL obs_int_h2d( 1, 1, &
1655	& zweig, zseaicel(:,:,iobs), zext )
1656
1657	seaicedatqc%rmod(jobs,1) = zext(1)
1658
1659	END DO
1660
1661	! Deallocate the data for interpolation
1662	DEALLOCATE( &
1663	& igrdi, &
1664	& igrdj, &
1665	& zglam, &
1666	& zgphi, &
1667	& zmask, &
1668	& zseaicel &
1669	& )
1670
1671	seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1672
1673	END SUBROUTINE obs_seaice_opt
1674
1675	SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1676	& pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1677	& ld_dailyav )
1678	!!-----------------------------------------------------------------------
1679	!!
1680	!! * ROUTINE obs_vel_opt *
1681	!!
1682	!! ** Purpose : Compute the model counterpart of velocity profile
1683	!! data by interpolating from the model grid to the
1684	!! observation point.
1685	!!
1686	!! ** Method : Linearly interpolate zonal and meridional components of velocity
1687	!! to each observation point using the model values at the corners of
1688	!! the surrounding grid box. The model velocity components are on a
1689	!! staggered C- grid.
1690	!!
1691	!! For velocity data from the TAO array, the model equivalent is
1692	!! a daily mean velocity field. So, we first compute
1693	!! the mean, then interpolate only at the end of the day.
1694	!!
1695	!! ** Action :
1696	!!
1697	!! History :
1698	!! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles
1699	!! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1700	!!-----------------------------------------------------------------------
1701
1702	!! * Modules used
1703	USE obs_profiles_def ! Definition of storage space for profile obs.
1704
1705	IMPLICIT NONE
1706
1707	!! * Arguments
1708	TYPE(obs_prof), INTENT(INOUT) :: &
1709	& prodatqc ! Subset of profile data not failing screening
1710	INTEGER, INTENT(IN) :: kt ! Time step
1711	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1712	INTEGER, INTENT(IN) :: kpj
1713	INTEGER, INTENT(IN) :: kpk
1714	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1715	! (kit000-1 = restart time)
1716	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
1717	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1718	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1719	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1720	& pun, & ! Model zonal component of velocity
1721	& pvn, & ! Model meridional component of velocity
1722	& pumask, & ! Land-sea mask
1723	& pvmask ! Land-sea mask
1724	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1725	& pgdept ! Model array of depth levels
1726	LOGICAL, INTENT(IN) :: ld_dailyav
1727
1728	!! * Local declarations
1729	INTEGER :: ji
1730	INTEGER :: jj
1731	INTEGER :: jk
1732	INTEGER :: jobs
1733	INTEGER :: inrc
1734	INTEGER :: ipro
1735	INTEGER :: idayend
1736	INTEGER :: ista
1737	INTEGER :: iend
1738	INTEGER :: iobs
1739	INTEGER, DIMENSION(imaxavtypes) :: &
1740	& idailyavtypes
1741	REAL(KIND=wp) :: zlam
1742	REAL(KIND=wp) :: zphi
1743	REAL(KIND=wp) :: zdaystp
1744	REAL(KIND=wp), DIMENSION(kpk) :: &
1745	& zobsmasku, &
1746	& zobsmaskv, &
1747	& zobsmask, &
1748	& zobsk, &
1749	& zobs2k
1750	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1751	& zweigu,zweigv
1752	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1753	& zumask, zvmask, &
1754	& zintu, &
1755	& zintv, &
1756	& zinmu, &
1757	& zinmv
1758	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1759	& zglamu, zglamv, &
1760	& zgphiu, zgphiv
1761	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1762	& igrdiu, &
1763	& igrdju, &
1764	& igrdiv, &
1765	& igrdjv
1766
1767	!------------------------------------------------------------------------
1768	! Local initialization
1769	!------------------------------------------------------------------------
1770	! ... Record and data counters
1771	inrc = kt - kit000 + 2
1772	ipro = prodatqc%npstp(inrc)
1773
1774	! Initialize daily mean for first timestep
1775	idayend = MOD( kt - kit000 + 1, kdaystp )
1776
1777	! Added kt == 0 test to catch restart case
1778	IF ( idayend == 1 .OR. kt == 0) THEN
1779	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1780	prodatqc%vdmean(:,:,:,1) = 0.0
1781	prodatqc%vdmean(:,:,:,2) = 0.0
1782	ENDIF
1783
1784	! Increment the zonal velocity field for computing daily mean
1785	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1786	! Increment the meridional velocity field for computing daily mean
1787	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1788
1789	! Compute the daily mean at the end of day
1790	zdaystp = 1.0 / REAL( kdaystp )
1791	IF ( idayend == 0 ) THEN
1792	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1793	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1794	ENDIF
1795
1796	! Get the data for interpolation
1797	ALLOCATE( &
1798	& igrdiu(2,2,ipro), &
1799	& igrdju(2,2,ipro), &
1800	& igrdiv(2,2,ipro), &
1801	& igrdjv(2,2,ipro), &
1802	& zglamu(2,2,ipro), zglamv(2,2,ipro), &
1803	& zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1804	& zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1805	& zintu(2,2,kpk,ipro), &
1806	& zintv(2,2,kpk,ipro) &
1807	& )
1808
1809	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1810	iobs = jobs - prodatqc%nprofup
1811	igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1812	igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1813	igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1814	igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1815	igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1816	igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1817	igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1818	igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1819	igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1820	igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1821	igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1822	igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1823	igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1824	igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1825	igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1826	igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1827	END DO
1828
1829	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1830	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1831	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1832	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1833
1834	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1835	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1836	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1837	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1838
1839	! At the end of the day also get interpolated means
1840	IF ( idayend == 0 ) THEN
1841
1842	ALLOCATE( &
1843	& zinmu(2,2,kpk,ipro), &
1844	& zinmv(2,2,kpk,ipro) &
1845	& )
1846
1847	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1848	& prodatqc%vdmean(:,:,:,1), zinmu )
1849	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1850	& prodatqc%vdmean(:,:,:,2), zinmv )
1851
1852	ENDIF
1853
1854	! loop over observations
1855
1856	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1857
1858	iobs = jobs - prodatqc%nprofup
1859
1860	IF ( kt /= prodatqc%mstp(jobs) ) THEN
1861
1862	IF(lwp) THEN
1863	WRITE(numout,*)
1864	WRITE(numout,*) ' E R R O R : Observation', &
1865	& ' time step is not consistent with the', &
1866	& ' model time step'
1867	WRITE(numout,*) ' ========='
1868	WRITE(numout,*)
1869	WRITE(numout,*) ' Record = ', jobs, &
1870	& ' kt = ', kt, &
1871	& ' mstp = ', prodatqc%mstp(jobs), &
1872	& ' ntyp = ', prodatqc%ntyp(jobs)
1873	ENDIF
1874	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1875	ENDIF
1876
1877	zlam = prodatqc%rlam(jobs)
1878	zphi = prodatqc%rphi(jobs)
1879
1880	! Initialize observation masks
1881
1882	zobsmasku(:) = 0.0
1883	zobsmaskv(:) = 0.0
1884
1885	! Horizontal weights and vertical mask
1886
1887	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1888
1889	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1890	& zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1891	& zumask(:,:,:,iobs), zweigu, zobsmasku )
1892
1893	ENDIF
1894
1895
1896	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1897
1898	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1899	& zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1900	& zvmask(:,:,:,iobs), zweigv, zobsmaskv )
1901
1902	ENDIF
1903
1904	! Ensure that the vertical mask on u and v are consistent.
1905
1906	zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1907
1908	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1909
1910	zobsk(:) = obfillflt
1911
1912	IF ( ld_dailyav ) THEN
1913
1914	IF ( idayend == 0 ) THEN
1915
1916	! Daily averaged data
1917
1918	CALL obs_int_h2d( kpk, kpk, &
1919	& zweigu, zinmu(:,:,:,iobs), zobsk )
1920
1921
1922	ELSE
1923
1924	CALL ctl_stop( ' A nonzero' // &
1925	& ' number of U profile data should' // &
1926	& ' only occur at the end of a given day' )
1927
1928	ENDIF
1929
1930	ELSE
1931
1932	! Point data
1933
1934	CALL obs_int_h2d( kpk, kpk, &
1935	& zweigu, zintu(:,:,:,iobs), zobsk )
1936
1937	ENDIF
1938
1939	!-------------------------------------------------------------
1940	! Compute vertical second-derivative of the interpolating
1941	! polynomial at obs points
1942	!-------------------------------------------------------------
1943
1944	IF ( k1dint == 1 ) THEN
1945	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1946	& pgdept, zobsmask )
1947	ENDIF
1948
1949	!-----------------------------------------------------------------
1950	! Vertical interpolation to the observation point
1951	!-----------------------------------------------------------------
1952	ista = prodatqc%npvsta(jobs,1)
1953	iend = prodatqc%npvend(jobs,1)
1954	CALL obs_int_z1d( kpk, &
1955	& prodatqc%var(1)%mvk(ista:iend), &
1956	& k1dint, iend - ista + 1, &
1957	& prodatqc%var(1)%vdep(ista:iend), &
1958	& zobsk, zobs2k, &
1959	& prodatqc%var(1)%vmod(ista:iend), &
1960	& pgdept, zobsmask )
1961
1962	ENDIF
1963
1964	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1965
1966	zobsk(:) = obfillflt
1967
1968	IF ( ld_dailyav ) THEN
1969
1970	IF ( idayend == 0 ) THEN
1971
1972	! Daily averaged data
1973
1974	CALL obs_int_h2d( kpk, kpk, &
1975	& zweigv, zinmv(:,:,:,iobs), zobsk )
1976
1977	ELSE
1978
1979	CALL ctl_stop( ' A nonzero' // &
1980	& ' number of V profile data should' // &
1981	& ' only occur at the end of a given day' )
1982
1983	ENDIF
1984
1985	ELSE
1986
1987	! Point data
1988
1989	CALL obs_int_h2d( kpk, kpk, &
1990	& zweigv, zintv(:,:,:,iobs), zobsk )
1991
1992	ENDIF
1993
1994
1995	!-------------------------------------------------------------
1996	! Compute vertical second-derivative of the interpolating
1997	! polynomial at obs points
1998	!-------------------------------------------------------------
1999
2000	IF ( k1dint == 1 ) THEN
2001	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
2002	& pgdept, zobsmask )
2003	ENDIF
2004
2005	!----------------------------------------------------------------
2006	! Vertical interpolation to the observation point
2007	!----------------------------------------------------------------
2008	ista = prodatqc%npvsta(jobs,2)
2009	iend = prodatqc%npvend(jobs,2)
2010	CALL obs_int_z1d( kpk, &
2011	& prodatqc%var(2)%mvk(ista:iend),&
2012	& k1dint, iend - ista + 1, &
2013	& prodatqc%var(2)%vdep(ista:iend),&
2014	& zobsk, zobs2k, &
2015	& prodatqc%var(2)%vmod(ista:iend),&
2016	& pgdept, zobsmask )
2017
2018	ENDIF
2019
2020	END DO
2021
2022	! Deallocate the data for interpolation
2023	DEALLOCATE( &
2024	& igrdiu, &
2025	& igrdju, &
2026	& igrdiv, &
2027	& igrdjv, &
2028	& zglamu, zglamv, &
2029	& zgphiu, zgphiv, &
2030	& zumask, zvmask, &
2031	& zintu, &
2032	& zintv &
2033	& )
2034	! At the end of the day also get interpolated means
2035	IF ( idayend == 0 ) THEN
2036	DEALLOCATE( &
2037	& zinmu, &
2038	& zinmv &
2039	& )
2040	ENDIF
2041
2042	prodatqc%nprofup = prodatqc%nprofup + ipro
2043
2044	END SUBROUTINE obs_vel_opt
2045
2046	END MODULE obs_oper
2047

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

Download in other formats: