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

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

obs_oper.F90 in branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

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

Last change on this file since 6494 was 6301, checked in by djlea, 8 years ago
Update to QC out observations which are caught by the strict bathymetry check in obs_oper. Reverting some changes to obs_prep which are no longer needed.
File size: 76.0 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,kpj,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	! Initiialise 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	! Set QC flag for any observations found below the bottom
888	! needed as the check here is more strict than that in obs_prep
889	IF (sum(l_zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
890
891	ENDDO
892
893
894	DEALLOCATE(interp_corner,iv_indic)
895
896	ENDIF
897
898
899	! Salinity
900
901	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
902
903	zobsk(:) = obfillflt
904
905	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
906
907	IF ( idayend == 0 ) THEN
908
909	! Daily averaged moored buoy (MRB) data
910
911	! vertically interpolate all 4 corners
912	ista = prodatqc%npvsta(jobs,2)
913	iend = prodatqc%npvend(jobs,2)
914	inum_obs = iend - ista + 1
915	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
916
917	DO iin=1,2
918	DO ijn=1,2
919
920
921
922	IF ( k1dint == 1 ) THEN
923	CALL obs_int_z1d_spl( kpk, &
924	& zinms(iin,ijn,:,iobs), &
925	& zobs2k, zgdept(iin,ijn,:,iobs), &
926	& zmask(iin,ijn,:,iobs))
927	ENDIF
928
929	CALL obs_level_search(kpk, &
930	& zgdept(iin,ijn,:,iobs), &
931	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
932	& iv_indic)
933	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
934	& prodatqc%var(2)%vdep(ista:iend), &
935	& zinms(iin,ijn,:,iobs), &
936	& zobs2k, interp_corner(iin,ijn,:), &
937	& zgdept(iin,ijn,:,iobs), &
938	& zmask(iin,ijn,:,iobs))
939
940	ENDDO
941	ENDDO
942
943
944	ELSE
945
946	CALL ctl_stop( ' A nonzero' // &
947	& ' number of profile T BUOY data should' // &
948	& ' only occur at the end of a given day' )
949
950	ENDIF
951
952	ELSE
953
954	! Point data
955
956	! vertically interpolate all 4 corners
957	ista = prodatqc%npvsta(jobs,2)
958	iend = prodatqc%npvend(jobs,2)
959	inum_obs = iend - ista + 1
960	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
961
962	DO iin=1,2
963	DO ijn=1,2
964
965
966	IF ( k1dint == 1 ) THEN
967	CALL obs_int_z1d_spl( kpk, &
968	& zints(iin,ijn,:,iobs),&
969	& zobs2k, zgdept(iin,ijn,:,iobs), &
970	& zmask(iin,ijn,:,iobs))
971
972	ENDIF
973
974	CALL obs_level_search(kpk, &
975	& zgdept(iin,ijn,:,iobs),&
976	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
977	& iv_indic)
978	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
979	& prodatqc%var(2)%vdep(ista:iend), &
980	& zints(iin,ijn,:,iobs), &
981	& zobs2k,interp_corner(iin,ijn,:), &
982	& zgdept(iin,ijn,:,iobs), &
983	& zmask(iin,ijn,:,iobs) )
984
985	ENDDO
986	ENDDO
987
988	ENDIF
989
990	!-------------------------------------------------------------
991	! Compute the horizontal interpolation for every profile level
992	!-------------------------------------------------------------
993
994	DO ikn=1,inum_obs
995	iend=ista+ikn-1
996
997	l_zweig(:,:,1) = 0._wp
998
999	! This code forces the horizontal weights to be
1000	! zero IF the observation is below the bottom of the
1001	! corners of the interpolation nodes, Or if it is in
1002	! the mask. This is important for observations are near
1003	! steep bathymetry
1004	DO iin=1,2
1005	DO ijn=1,2
1006
1007	depth_loop2: DO ik=kpk,2,-1
1008	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
1009
1010	l_zweig(iin,ijn,1) = &
1011	& zweig(iin,ijn,1) * &
1012	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
1013	& - prodatqc%var(2)%vdep(iend)),0._wp)
1014
1015	EXIT depth_loop2
1016	ENDIF
1017	ENDDO depth_loop2
1018
1019	ENDDO
1020	ENDDO
1021
1022	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
1023	& prodatqc%var(2)%vmod(iend:iend) )
1024
1025	! Set QC flag for any observations found below the bottom
1026	! needed as the check here is more strict than that in obs_prep
1027	IF (sum(l_zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
1028
1029	ENDDO
1030
1031
1032	DEALLOCATE(interp_corner,iv_indic)
1033
1034	ENDIF
1035
1036	ENDIF
1037
1038	END DO
1039
1040	! Deallocate the data for interpolation
1041	DEALLOCATE( &
1042	& igrdi, &
1043	& igrdj, &
1044	& zglam, &
1045	& zgphi, &
1046	& zmask, &
1047	& zintt, &
1048	& zints, &
1049	& zgdept,&
1050	& zgdepw &
1051	& )
1052	! At the end of the day also get interpolated means
1053	IF ( idayend == 0 ) THEN
1054	DEALLOCATE( &
1055	& zinmt, &
1056	& zinms &
1057	& )
1058	ENDIF
1059
1060	prodatqc%nprofup = prodatqc%nprofup + ipro
1061
1062	END SUBROUTINE obs_pro_sco_opt
1063
1064	SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
1065	& psshn, psshmask, k2dint )
1066	!!-----------------------------------------------------------------------
1067	!!
1068	!! * ROUTINE obs_sla_opt *
1069	!!
1070	!! ** Purpose : Compute the model counterpart of sea level anomaly
1071	!! data by interpolating from the model grid to the
1072	!! observation point.
1073	!!
1074	!! ** Method : Linearly interpolate to each observation point using
1075	!! the model values at the corners of the surrounding grid box.
1076	!!
1077	!! The now model SSH is first computed at the obs (lon, lat) point.
1078	!!
1079	!! Several horizontal interpolation schemes are available:
1080	!! - distance-weighted (great circle) (k2dint = 0)
1081	!! - distance-weighted (small angle) (k2dint = 1)
1082	!! - bilinear (geographical grid) (k2dint = 2)
1083	!! - bilinear (quadrilateral grid) (k2dint = 3)
1084	!! - polynomial (quadrilateral grid) (k2dint = 4)
1085	!!
1086	!! The sea level anomaly at the observation points is then computed
1087	!! by removing a mean dynamic topography (defined at the obs. point).
1088	!!
1089	!! ** Action :
1090	!!
1091	!! History :
1092	!! ! 07-03 (A. Weaver)
1093	!!-----------------------------------------------------------------------
1094
1095	!! * Modules used
1096	USE obs_surf_def ! Definition of storage space for surface observations
1097
1098	IMPLICIT NONE
1099
1100	!! * Arguments
1101	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening
1102	INTEGER, INTENT(IN) :: kt ! Time step
1103	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1104	INTEGER, INTENT(IN) :: kpj
1105	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1106	! (kit000-1 = restart time)
1107	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1108	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1109	& psshn, & ! Model SSH field
1110	& psshmask ! Land-sea mask
1111
1112	!! * Local declarations
1113	INTEGER :: ji
1114	INTEGER :: jj
1115	INTEGER :: jobs
1116	INTEGER :: inrc
1117	INTEGER :: isla
1118	INTEGER :: iobs
1119	REAL(KIND=wp) :: zlam
1120	REAL(KIND=wp) :: zphi
1121	REAL(KIND=wp) :: zext(1), zobsmask(1)
1122	REAL(kind=wp), DIMENSION(2,2,1) :: &
1123	& zweig
1124	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1125	& zmask, &
1126	& zsshl, &
1127	& zglam, &
1128	& zgphi
1129	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1130	& igrdi, &
1131	& igrdj
1132
1133	!------------------------------------------------------------------------
1134	! Local initialization
1135	!------------------------------------------------------------------------
1136	! ... Record and data counters
1137	inrc = kt - kit000 + 2
1138	isla = sladatqc%nsstp(inrc)
1139
1140	! Get the data for interpolation
1141
1142	ALLOCATE( &
1143	& igrdi(2,2,isla), &
1144	& igrdj(2,2,isla), &
1145	& zglam(2,2,isla), &
1146	& zgphi(2,2,isla), &
1147	& zmask(2,2,isla), &
1148	& zsshl(2,2,isla) &
1149	& )
1150
1151	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1152	iobs = jobs - sladatqc%nsurfup
1153	igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
1154	igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
1155	igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
1156	igrdj(1,2,iobs) = sladatqc%mj(jobs)
1157	igrdi(2,1,iobs) = sladatqc%mi(jobs)
1158	igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
1159	igrdi(2,2,iobs) = sladatqc%mi(jobs)
1160	igrdj(2,2,iobs) = sladatqc%mj(jobs)
1161	END DO
1162
1163	CALL obs_int_comm_2d( 2, 2, isla, &
1164	& igrdi, igrdj, glamt, zglam )
1165	CALL obs_int_comm_2d( 2, 2, isla, &
1166	& igrdi, igrdj, gphit, zgphi )
1167	CALL obs_int_comm_2d( 2, 2, isla, &
1168	& igrdi, igrdj, psshmask, zmask )
1169	CALL obs_int_comm_2d( 2, 2, isla, &
1170	& igrdi, igrdj, psshn, zsshl )
1171
1172	! Loop over observations
1173
1174	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1175
1176	iobs = jobs - sladatqc%nsurfup
1177
1178	IF ( kt /= sladatqc%mstp(jobs) ) THEN
1179
1180	IF(lwp) THEN
1181	WRITE(numout,*)
1182	WRITE(numout,*) ' E R R O R : Observation', &
1183	& ' time step is not consistent with the', &
1184	& ' model time step'
1185	WRITE(numout,*) ' ========='
1186	WRITE(numout,*)
1187	WRITE(numout,*) ' Record = ', jobs, &
1188	& ' kt = ', kt, &
1189	& ' mstp = ', sladatqc%mstp(jobs), &
1190	& ' ntyp = ', sladatqc%ntyp(jobs)
1191	ENDIF
1192	CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
1193
1194	ENDIF
1195
1196	zlam = sladatqc%rlam(jobs)
1197	zphi = sladatqc%rphi(jobs)
1198
1199	! Get weights to interpolate the model SSH to the observation point
1200	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1201	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1202	& zmask(:,:,iobs), zweig, zobsmask )
1203
1204
1205	! Interpolate the model SSH to the observation point
1206	CALL obs_int_h2d( 1, 1, &
1207	& zweig, zsshl(:,:,iobs), zext )
1208
1209	sladatqc%rext(jobs,1) = zext(1)
1210	! ... Remove the MDT at the observation point
1211	sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
1212
1213	END DO
1214
1215	! Deallocate the data for interpolation
1216	DEALLOCATE( &
1217	& igrdi, &
1218	& igrdj, &
1219	& zglam, &
1220	& zgphi, &
1221	& zmask, &
1222	& zsshl &
1223	& )
1224
1225	sladatqc%nsurfup = sladatqc%nsurfup + isla
1226
1227	END SUBROUTINE obs_sla_opt
1228
1229	SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
1230	& psstn, psstmask, k2dint, ld_nightav )
1231	!!-----------------------------------------------------------------------
1232	!!
1233	!! * ROUTINE obs_sst_opt *
1234	!!
1235	!! ** Purpose : Compute the model counterpart of surface temperature
1236	!! data by interpolating from the model grid to the
1237	!! observation point.
1238	!!
1239	!! ** Method : Linearly interpolate to each observation point using
1240	!! the model values at the corners of the surrounding grid box.
1241	!!
1242	!! The now model SST is first computed at the obs (lon, lat) point.
1243	!!
1244	!! Several horizontal interpolation schemes are available:
1245	!! - distance-weighted (great circle) (k2dint = 0)
1246	!! - distance-weighted (small angle) (k2dint = 1)
1247	!! - bilinear (geographical grid) (k2dint = 2)
1248	!! - bilinear (quadrilateral grid) (k2dint = 3)
1249	!! - polynomial (quadrilateral grid) (k2dint = 4)
1250	!!
1251	!!
1252	!! ** Action :
1253	!!
1254	!! History :
1255	!! ! 07-07 (S. Ricci ) : Original
1256	!!
1257	!!-----------------------------------------------------------------------
1258
1259	!! * Modules used
1260	USE obs_surf_def ! Definition of storage space for surface observations
1261	USE sbcdcy
1262
1263	IMPLICIT NONE
1264
1265	!! * Arguments
1266	TYPE(obs_surf), INTENT(INOUT) :: &
1267	& sstdatqc ! Subset of surface data not failing screening
1268	INTEGER, INTENT(IN) :: kt ! Time step
1269	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1270	INTEGER, INTENT(IN) :: kpj
1271	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1272	! (kit000-1 = restart time)
1273	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1274	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1275	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1276	& psstn, & ! Model SST field
1277	& psstmask ! Land-sea mask
1278
1279	!! * Local declarations
1280	INTEGER :: ji
1281	INTEGER :: jj
1282	INTEGER :: jobs
1283	INTEGER :: inrc
1284	INTEGER :: isst
1285	INTEGER :: iobs
1286	INTEGER :: idayend
1287	REAL(KIND=wp) :: zlam
1288	REAL(KIND=wp) :: zphi
1289	REAL(KIND=wp) :: zext(1), zobsmask(1)
1290	REAL(KIND=wp) :: zdaystp
1291	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1292	& icount_sstnight, &
1293	& imask_night
1294	REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1295	& zintmp, &
1296	& zouttmp, &
1297	& zmeanday ! to compute model sst in region of 24h daylight (pole)
1298	REAL(kind=wp), DIMENSION(2,2,1) :: &
1299	& zweig
1300	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1301	& zmask, &
1302	& zsstl, &
1303	& zsstm, &
1304	& zglam, &
1305	& zgphi
1306	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1307	& igrdi, &
1308	& igrdj
1309	LOGICAL, INTENT(IN) :: ld_nightav
1310
1311	!-----------------------------------------------------------------------
1312	! Local initialization
1313	!-----------------------------------------------------------------------
1314	! ... Record and data counters
1315	inrc = kt - kit000 + 2
1316	isst = sstdatqc%nsstp(inrc)
1317
1318	IF ( ld_nightav ) THEN
1319
1320	! Initialize array for night mean
1321
1322	IF ( kt .EQ. 0 ) THEN
1323	ALLOCATE ( icount_sstnight(kpi,kpj) )
1324	ALLOCATE ( imask_night(kpi,kpj) )
1325	ALLOCATE ( zintmp(kpi,kpj) )
1326	ALLOCATE ( zouttmp(kpi,kpj) )
1327	ALLOCATE ( zmeanday(kpi,kpj) )
1328	nday_qsr = -1 ! initialisation flag for nbc_dcy
1329	ENDIF
1330
1331	! Initialize daily mean for first timestep
1332	idayend = MOD( kt - kit000 + 1, kdaystp )
1333
1334	! Added kt == 0 test to catch restart case
1335	IF ( idayend == 1 .OR. kt == 0) THEN
1336	IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
1337	DO jj = 1, jpj
1338	DO ji = 1, jpi
1339	sstdatqc%vdmean(ji,jj) = 0.0
1340	zmeanday(ji,jj) = 0.0
1341	icount_sstnight(ji,jj) = 0
1342	END DO
1343	END DO
1344	ENDIF
1345
1346	zintmp(:,:) = 0.0
1347	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1348	imask_night(:,:) = INT( zouttmp(:,:) )
1349
1350	DO jj = 1, jpj
1351	DO ji = 1, jpi
1352	! Increment the temperature field for computing night mean and counter
1353	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1354	& + psstn(ji,jj)*imask_night(ji,jj)
1355	zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj)
1356	icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
1357	END DO
1358	END DO
1359
1360	! Compute the daily mean at the end of day
1361
1362	zdaystp = 1.0 / REAL( kdaystp )
1363
1364	IF ( idayend == 0 ) THEN
1365	DO jj = 1, jpj
1366	DO ji = 1, jpi
1367	! Test if "no night" point
1368	IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
1369	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1370	& / icount_sstnight(ji,jj)
1371	ELSE
1372	sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1373	ENDIF
1374	END DO
1375	END DO
1376	ENDIF
1377
1378	ENDIF
1379
1380	! Get the data for interpolation
1381
1382	ALLOCATE( &
1383	& igrdi(2,2,isst), &
1384	& igrdj(2,2,isst), &
1385	& zglam(2,2,isst), &
1386	& zgphi(2,2,isst), &
1387	& zmask(2,2,isst), &
1388	& zsstl(2,2,isst) &
1389	& )
1390
1391	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1392	iobs = jobs - sstdatqc%nsurfup
1393	igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
1394	igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
1395	igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
1396	igrdj(1,2,iobs) = sstdatqc%mj(jobs)
1397	igrdi(2,1,iobs) = sstdatqc%mi(jobs)
1398	igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
1399	igrdi(2,2,iobs) = sstdatqc%mi(jobs)
1400	igrdj(2,2,iobs) = sstdatqc%mj(jobs)
1401	END DO
1402
1403	CALL obs_int_comm_2d( 2, 2, isst, &
1404	& igrdi, igrdj, glamt, zglam )
1405	CALL obs_int_comm_2d( 2, 2, isst, &
1406	& igrdi, igrdj, gphit, zgphi )
1407	CALL obs_int_comm_2d( 2, 2, isst, &
1408	& igrdi, igrdj, psstmask, zmask )
1409	CALL obs_int_comm_2d( 2, 2, isst, &
1410	& igrdi, igrdj, psstn, zsstl )
1411
1412	! At the end of the day get interpolated means
1413	IF ( idayend == 0 .AND. ld_nightav ) THEN
1414
1415	ALLOCATE( &
1416	& zsstm(2,2,isst) &
1417	& )
1418
1419	CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
1420	& sstdatqc%vdmean(:,:), zsstm )
1421
1422	ENDIF
1423
1424	! Loop over observations
1425
1426	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1427
1428	iobs = jobs - sstdatqc%nsurfup
1429
1430	IF ( kt /= sstdatqc%mstp(jobs) ) THEN
1431
1432	IF(lwp) THEN
1433	WRITE(numout,*)
1434	WRITE(numout,*) ' E R R O R : Observation', &
1435	& ' time step is not consistent with the', &
1436	& ' model time step'
1437	WRITE(numout,*) ' ========='
1438	WRITE(numout,*)
1439	WRITE(numout,*) ' Record = ', jobs, &
1440	& ' kt = ', kt, &
1441	& ' mstp = ', sstdatqc%mstp(jobs), &
1442	& ' ntyp = ', sstdatqc%ntyp(jobs)
1443	ENDIF
1444	CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
1445
1446	ENDIF
1447
1448	zlam = sstdatqc%rlam(jobs)
1449	zphi = sstdatqc%rphi(jobs)
1450
1451	! Get weights to interpolate the model SST to the observation point
1452	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1453	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1454	& zmask(:,:,iobs), zweig, zobsmask )
1455
1456	! Interpolate the model SST to the observation point
1457
1458	IF ( ld_nightav ) THEN
1459
1460	IF ( idayend == 0 ) THEN
1461	! Daily averaged/diurnal cycle of SST data
1462	CALL obs_int_h2d( 1, 1, &
1463	& zweig, zsstm(:,:,iobs), zext )
1464	ELSE
1465	CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // &
1466	& ' number of night SST data should' // &
1467	& ' only occur at the end of a given day' )
1468	ENDIF
1469
1470	ELSE
1471
1472	CALL obs_int_h2d( 1, 1, &
1473	& zweig, zsstl(:,:,iobs), zext )
1474
1475	ENDIF
1476	sstdatqc%rmod(jobs,1) = zext(1)
1477
1478	END DO
1479
1480	! Deallocate the data for interpolation
1481	DEALLOCATE( &
1482	& igrdi, &
1483	& igrdj, &
1484	& zglam, &
1485	& zgphi, &
1486	& zmask, &
1487	& zsstl &
1488	& )
1489
1490	! At the end of the day also get interpolated means
1491	IF ( idayend == 0 .AND. ld_nightav ) THEN
1492	DEALLOCATE( &
1493	& zsstm &
1494	& )
1495	ENDIF
1496
1497	sstdatqc%nsurfup = sstdatqc%nsurfup + isst
1498
1499	END SUBROUTINE obs_sst_opt
1500
1501	SUBROUTINE obs_sss_opt
1502	!!-----------------------------------------------------------------------
1503	!!
1504	!! * ROUTINE obs_sss_opt *
1505	!!
1506	!! ** Purpose : Compute the model counterpart of sea surface salinity
1507	!! data by interpolating from the model grid to the
1508	!! observation point.
1509	!!
1510	!! ** Method :
1511	!!
1512	!! ** Action :
1513	!!
1514	!! History :
1515	!! ! ??-??
1516	!!-----------------------------------------------------------------------
1517
1518	IMPLICIT NONE
1519
1520	END SUBROUTINE obs_sss_opt
1521
1522	SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
1523	& pseaicen, pseaicemask, k2dint )
1524
1525	!!-----------------------------------------------------------------------
1526	!!
1527	!! * ROUTINE obs_seaice_opt *
1528	!!
1529	!! ** Purpose : Compute the model counterpart of surface temperature
1530	!! data by interpolating from the model grid to the
1531	!! observation point.
1532	!!
1533	!! ** Method : Linearly interpolate to each observation point using
1534	!! the model values at the corners of the surrounding grid box.
1535	!!
1536	!! The now model sea ice is first computed at the obs (lon, lat) point.
1537	!!
1538	!! Several horizontal interpolation schemes are available:
1539	!! - distance-weighted (great circle) (k2dint = 0)
1540	!! - distance-weighted (small angle) (k2dint = 1)
1541	!! - bilinear (geographical grid) (k2dint = 2)
1542	!! - bilinear (quadrilateral grid) (k2dint = 3)
1543	!! - polynomial (quadrilateral grid) (k2dint = 4)
1544	!!
1545	!!
1546	!! ** Action :
1547	!!
1548	!! History :
1549	!! ! 07-07 (S. Ricci ) : Original
1550	!!
1551	!!-----------------------------------------------------------------------
1552
1553	!! * Modules used
1554	USE obs_surf_def ! Definition of storage space for surface observations
1555
1556	IMPLICIT NONE
1557
1558	!! * Arguments
1559	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening
1560	INTEGER, INTENT(IN) :: kt ! Time step
1561	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1562	INTEGER, INTENT(IN) :: kpj
1563	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1564	! (kit000-1 = restart time)
1565	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1566	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1567	& pseaicen, & ! Model sea ice field
1568	& pseaicemask ! Land-sea mask
1569
1570	!! * Local declarations
1571	INTEGER :: ji
1572	INTEGER :: jj
1573	INTEGER :: jobs
1574	INTEGER :: inrc
1575	INTEGER :: iseaice
1576	INTEGER :: iobs
1577
1578	REAL(KIND=wp) :: zlam
1579	REAL(KIND=wp) :: zphi
1580	REAL(KIND=wp) :: zext(1), zobsmask(1)
1581	REAL(kind=wp), DIMENSION(2,2,1) :: &
1582	& zweig
1583	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1584	& zmask, &
1585	& zseaicel, &
1586	& zglam, &
1587	& zgphi
1588	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1589	& igrdi, &
1590	& igrdj
1591
1592	!------------------------------------------------------------------------
1593	! Local initialization
1594	!------------------------------------------------------------------------
1595	! ... Record and data counters
1596	inrc = kt - kit000 + 2
1597	iseaice = seaicedatqc%nsstp(inrc)
1598
1599	! Get the data for interpolation
1600
1601	ALLOCATE( &
1602	& igrdi(2,2,iseaice), &
1603	& igrdj(2,2,iseaice), &
1604	& zglam(2,2,iseaice), &
1605	& zgphi(2,2,iseaice), &
1606	& zmask(2,2,iseaice), &
1607	& zseaicel(2,2,iseaice) &
1608	& )
1609
1610	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1611	iobs = jobs - seaicedatqc%nsurfup
1612	igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1613	igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1614	igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1615	igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1616	igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1617	igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1618	igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1619	igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1620	END DO
1621
1622	CALL obs_int_comm_2d( 2, 2, iseaice, &
1623	& igrdi, igrdj, glamt, zglam )
1624	CALL obs_int_comm_2d( 2, 2, iseaice, &
1625	& igrdi, igrdj, gphit, zgphi )
1626	CALL obs_int_comm_2d( 2, 2, iseaice, &
1627	& igrdi, igrdj, pseaicemask, zmask )
1628	CALL obs_int_comm_2d( 2, 2, iseaice, &
1629	& igrdi, igrdj, pseaicen, zseaicel )
1630
1631	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1632
1633	iobs = jobs - seaicedatqc%nsurfup
1634
1635	IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1636
1637	IF(lwp) THEN
1638	WRITE(numout,*)
1639	WRITE(numout,*) ' E R R O R : Observation', &
1640	& ' time step is not consistent with the', &
1641	& ' model time step'
1642	WRITE(numout,*) ' ========='
1643	WRITE(numout,*)
1644	WRITE(numout,*) ' Record = ', jobs, &
1645	& ' kt = ', kt, &
1646	& ' mstp = ', seaicedatqc%mstp(jobs), &
1647	& ' ntyp = ', seaicedatqc%ntyp(jobs)
1648	ENDIF
1649	CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1650
1651	ENDIF
1652
1653	zlam = seaicedatqc%rlam(jobs)
1654	zphi = seaicedatqc%rphi(jobs)
1655
1656	! Get weights to interpolate the model sea ice to the observation point
1657	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1658	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1659	& zmask(:,:,iobs), zweig, zobsmask )
1660
1661	! ... Interpolate the model sea ice to the observation point
1662	CALL obs_int_h2d( 1, 1, &
1663	& zweig, zseaicel(:,:,iobs), zext )
1664
1665	seaicedatqc%rmod(jobs,1) = zext(1)
1666
1667	END DO
1668
1669	! Deallocate the data for interpolation
1670	DEALLOCATE( &
1671	& igrdi, &
1672	& igrdj, &
1673	& zglam, &
1674	& zgphi, &
1675	& zmask, &
1676	& zseaicel &
1677	& )
1678
1679	seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1680
1681	END SUBROUTINE obs_seaice_opt
1682
1683	SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1684	& pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1685	& ld_dailyav )
1686	!!-----------------------------------------------------------------------
1687	!!
1688	!! * ROUTINE obs_vel_opt *
1689	!!
1690	!! ** Purpose : Compute the model counterpart of velocity profile
1691	!! data by interpolating from the model grid to the
1692	!! observation point.
1693	!!
1694	!! ** Method : Linearly interpolate zonal and meridional components of velocity
1695	!! to each observation point using the model values at the corners of
1696	!! the surrounding grid box. The model velocity components are on a
1697	!! staggered C- grid.
1698	!!
1699	!! For velocity data from the TAO array, the model equivalent is
1700	!! a daily mean velocity field. So, we first compute
1701	!! the mean, then interpolate only at the end of the day.
1702	!!
1703	!! ** Action :
1704	!!
1705	!! History :
1706	!! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles
1707	!! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1708	!!-----------------------------------------------------------------------
1709
1710	!! * Modules used
1711	USE obs_profiles_def ! Definition of storage space for profile obs.
1712
1713	IMPLICIT NONE
1714
1715	!! * Arguments
1716	TYPE(obs_prof), INTENT(INOUT) :: &
1717	& prodatqc ! Subset of profile data not failing screening
1718	INTEGER, INTENT(IN) :: kt ! Time step
1719	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1720	INTEGER, INTENT(IN) :: kpj
1721	INTEGER, INTENT(IN) :: kpk
1722	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1723	! (kit000-1 = restart time)
1724	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
1725	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1726	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1727	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1728	& pun, & ! Model zonal component of velocity
1729	& pvn, & ! Model meridional component of velocity
1730	& pumask, & ! Land-sea mask
1731	& pvmask ! Land-sea mask
1732	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1733	& pgdept ! Model array of depth levels
1734	LOGICAL, INTENT(IN) :: ld_dailyav
1735
1736	!! * Local declarations
1737	INTEGER :: ji
1738	INTEGER :: jj
1739	INTEGER :: jk
1740	INTEGER :: jobs
1741	INTEGER :: inrc
1742	INTEGER :: ipro
1743	INTEGER :: idayend
1744	INTEGER :: ista
1745	INTEGER :: iend
1746	INTEGER :: iobs
1747	INTEGER, DIMENSION(imaxavtypes) :: &
1748	& idailyavtypes
1749	REAL(KIND=wp) :: zlam
1750	REAL(KIND=wp) :: zphi
1751	REAL(KIND=wp) :: zdaystp
1752	REAL(KIND=wp), DIMENSION(kpk) :: &
1753	& zobsmasku, &
1754	& zobsmaskv, &
1755	& zobsmask, &
1756	& zobsk, &
1757	& zobs2k
1758	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1759	& zweigu,zweigv
1760	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1761	& zumask, zvmask, &
1762	& zintu, &
1763	& zintv, &
1764	& zinmu, &
1765	& zinmv
1766	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1767	& zglamu, zglamv, &
1768	& zgphiu, zgphiv
1769	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1770	& igrdiu, &
1771	& igrdju, &
1772	& igrdiv, &
1773	& igrdjv
1774
1775	!------------------------------------------------------------------------
1776	! Local initialization
1777	!------------------------------------------------------------------------
1778	! ... Record and data counters
1779	inrc = kt - kit000 + 2
1780	ipro = prodatqc%npstp(inrc)
1781
1782	! Initialize daily mean for first timestep
1783	idayend = MOD( kt - kit000 + 1, kdaystp )
1784
1785	! Added kt == 0 test to catch restart case
1786	IF ( idayend == 1 .OR. kt == 0) THEN
1787	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1788	prodatqc%vdmean(:,:,:,1) = 0.0
1789	prodatqc%vdmean(:,:,:,2) = 0.0
1790	ENDIF
1791
1792	! Increment the zonal velocity field for computing daily mean
1793	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1794	! Increment the meridional velocity field for computing daily mean
1795	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1796
1797	! Compute the daily mean at the end of day
1798	zdaystp = 1.0 / REAL( kdaystp )
1799	IF ( idayend == 0 ) THEN
1800	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1801	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1802	ENDIF
1803
1804	! Get the data for interpolation
1805	ALLOCATE( &
1806	& igrdiu(2,2,ipro), &
1807	& igrdju(2,2,ipro), &
1808	& igrdiv(2,2,ipro), &
1809	& igrdjv(2,2,ipro), &
1810	& zglamu(2,2,ipro), zglamv(2,2,ipro), &
1811	& zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1812	& zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1813	& zintu(2,2,kpk,ipro), &
1814	& zintv(2,2,kpk,ipro) &
1815	& )
1816
1817	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1818	iobs = jobs - prodatqc%nprofup
1819	igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1820	igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1821	igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1822	igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1823	igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1824	igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1825	igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1826	igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1827	igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1828	igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1829	igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1830	igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1831	igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1832	igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1833	igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1834	igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1835	END DO
1836
1837	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1838	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1839	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1840	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1841
1842	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1843	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1844	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1845	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1846
1847	! At the end of the day also get interpolated means
1848	IF ( idayend == 0 ) THEN
1849
1850	ALLOCATE( &
1851	& zinmu(2,2,kpk,ipro), &
1852	& zinmv(2,2,kpk,ipro) &
1853	& )
1854
1855	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1856	& prodatqc%vdmean(:,:,:,1), zinmu )
1857	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1858	& prodatqc%vdmean(:,:,:,2), zinmv )
1859
1860	ENDIF
1861
1862	! loop over observations
1863
1864	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1865
1866	iobs = jobs - prodatqc%nprofup
1867
1868	IF ( kt /= prodatqc%mstp(jobs) ) THEN
1869
1870	IF(lwp) THEN
1871	WRITE(numout,*)
1872	WRITE(numout,*) ' E R R O R : Observation', &
1873	& ' time step is not consistent with the', &
1874	& ' model time step'
1875	WRITE(numout,*) ' ========='
1876	WRITE(numout,*)
1877	WRITE(numout,*) ' Record = ', jobs, &
1878	& ' kt = ', kt, &
1879	& ' mstp = ', prodatqc%mstp(jobs), &
1880	& ' ntyp = ', prodatqc%ntyp(jobs)
1881	ENDIF
1882	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1883	ENDIF
1884
1885	zlam = prodatqc%rlam(jobs)
1886	zphi = prodatqc%rphi(jobs)
1887
1888	! Initialize observation masks
1889
1890	zobsmasku(:) = 0.0
1891	zobsmaskv(:) = 0.0
1892
1893	! Horizontal weights and vertical mask
1894
1895	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1896
1897	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1898	& zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1899	& zumask(:,:,:,iobs), zweigu, zobsmasku )
1900
1901	ENDIF
1902
1903
1904	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1905
1906	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1907	& zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1908	& zvmask(:,:,:,iobs), zweigv, zobsmaskv )
1909
1910	ENDIF
1911
1912	! Ensure that the vertical mask on u and v are consistent.
1913
1914	zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1915
1916	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1917
1918	zobsk(:) = obfillflt
1919
1920	IF ( ld_dailyav ) THEN
1921
1922	IF ( idayend == 0 ) THEN
1923
1924	! Daily averaged data
1925
1926	CALL obs_int_h2d( kpk, kpk, &
1927	& zweigu, zinmu(:,:,:,iobs), zobsk )
1928
1929
1930	ELSE
1931
1932	CALL ctl_stop( ' A nonzero' // &
1933	& ' number of U profile data should' // &
1934	& ' only occur at the end of a given day' )
1935
1936	ENDIF
1937
1938	ELSE
1939
1940	! Point data
1941
1942	CALL obs_int_h2d( kpk, kpk, &
1943	& zweigu, zintu(:,:,:,iobs), zobsk )
1944
1945	ENDIF
1946
1947	!-------------------------------------------------------------
1948	! Compute vertical second-derivative of the interpolating
1949	! polynomial at obs points
1950	!-------------------------------------------------------------
1951
1952	IF ( k1dint == 1 ) THEN
1953	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1954	& pgdept, zobsmask )
1955	ENDIF
1956
1957	!-----------------------------------------------------------------
1958	! Vertical interpolation to the observation point
1959	!-----------------------------------------------------------------
1960	ista = prodatqc%npvsta(jobs,1)
1961	iend = prodatqc%npvend(jobs,1)
1962	CALL obs_int_z1d( kpk, &
1963	& prodatqc%var(1)%mvk(ista:iend), &
1964	& k1dint, iend - ista + 1, &
1965	& prodatqc%var(1)%vdep(ista:iend), &
1966	& zobsk, zobs2k, &
1967	& prodatqc%var(1)%vmod(ista:iend), &
1968	& pgdept, zobsmask )
1969
1970	ENDIF
1971
1972	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1973
1974	zobsk(:) = obfillflt
1975
1976	IF ( ld_dailyav ) THEN
1977
1978	IF ( idayend == 0 ) THEN
1979
1980	! Daily averaged data
1981
1982	CALL obs_int_h2d( kpk, kpk, &
1983	& zweigv, zinmv(:,:,:,iobs), zobsk )
1984
1985	ELSE
1986
1987	CALL ctl_stop( ' A nonzero' // &
1988	& ' number of V profile data should' // &
1989	& ' only occur at the end of a given day' )
1990
1991	ENDIF
1992
1993	ELSE
1994
1995	! Point data
1996
1997	CALL obs_int_h2d( kpk, kpk, &
1998	& zweigv, zintv(:,:,:,iobs), zobsk )
1999
2000	ENDIF
2001
2002
2003	!-------------------------------------------------------------
2004	! Compute vertical second-derivative of the interpolating
2005	! polynomial at obs points
2006	!-------------------------------------------------------------
2007
2008	IF ( k1dint == 1 ) THEN
2009	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
2010	& pgdept, zobsmask )
2011	ENDIF
2012
2013	!----------------------------------------------------------------
2014	! Vertical interpolation to the observation point
2015	!----------------------------------------------------------------
2016	ista = prodatqc%npvsta(jobs,2)
2017	iend = prodatqc%npvend(jobs,2)
2018	CALL obs_int_z1d( kpk, &
2019	& prodatqc%var(2)%mvk(ista:iend),&
2020	& k1dint, iend - ista + 1, &
2021	& prodatqc%var(2)%vdep(ista:iend),&
2022	& zobsk, zobs2k, &
2023	& prodatqc%var(2)%vmod(ista:iend),&
2024	& pgdept, zobsmask )
2025
2026	ENDIF
2027
2028	END DO
2029
2030	! Deallocate the data for interpolation
2031	DEALLOCATE( &
2032	& igrdiu, &
2033	& igrdju, &
2034	& igrdiv, &
2035	& igrdjv, &
2036	& zglamu, zglamv, &
2037	& zgphiu, zgphiv, &
2038	& zumask, zvmask, &
2039	& zintu, &
2040	& zintv &
2041	& )
2042	! At the end of the day also get interpolated means
2043	IF ( idayend == 0 ) THEN
2044	DEALLOCATE( &
2045	& zinmu, &
2046	& zinmv &
2047	& )
2048	ENDIF
2049
2050	prodatqc%nprofup = prodatqc%nprofup + ipro
2051
2052	END SUBROUTINE obs_vel_opt
2053
2054	END MODULE obs_oper
2055

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

Download in other formats: