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/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 4746

Last change on this file since 4746 was 4746, checked in by jwhile, 10 years ago
Updated to allow interpolation in s-coordinates
Property svn:keywords set to `Id`
File size: 75.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, 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 levels
551	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
552	& kdailyavtypes ! Types for daily averages
553	!! * Local declarations
554	INTEGER :: ji
555	INTEGER :: jj
556	INTEGER :: jk
557	INTEGER :: iico, ijco
558	INTEGER :: jobs
559	INTEGER :: inrc
560	INTEGER :: ipro
561	INTEGER :: idayend
562	INTEGER :: ista
563	INTEGER :: iend
564	INTEGER :: iobs
565	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
566	INTEGER, DIMENSION(imaxavtypes) :: &
567	& idailyavtypes
568	REAL(KIND=wp) :: zlam
569	REAL(KIND=wp) :: zphi
570	REAL(KIND=wp) :: zdaystp
571	REAL(KIND=wp), DIMENSION(kpk) :: &
572	& zobsmask, &
573	& zobsk, &
574	& zobs2k
575	REAL(KIND=wp), DIMENSION(2,2,1) :: &
576	& zweig, &
577	& l_zweig
578	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
579	& zmask, &
580	& zintt, &
581	& zints, &
582	& zinmt, &
583	& zgdept,&
584	& zgdepw,&
585	& zinms
586	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
587	& zglam, &
588	& zgphi
589	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
590	& igrdi, &
591	& igrdj
592	INTEGER :: &
593	& inum_obs
594	REAL(KIND=wp), DIMENSION(1) :: zmsk_1
595	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
596	INTEGER, ALLOCATABLE, DIMENSION(:) :: v_indic
597
598	!------------------------------------------------------------------------
599	! Local initialization
600	!------------------------------------------------------------------------
601	! ... Record and data counters
602	inrc = kt - kit000 + 2
603	ipro = prodatqc%npstp(inrc)
604
605	! Daily average types
606	IF ( PRESENT(kdailyavtypes) ) THEN
607	idailyavtypes(:) = kdailyavtypes(:)
608	ELSE
609	idailyavtypes(:) = -1
610	ENDIF
611
612	! Initialize daily mean for first time-step
613	idayend = MOD( kt - kit000 + 1, kdaystp )
614
615	! Added kt == 0 test to catch restart case
616	IF ( idayend == 1 .OR. kt == 0) THEN
617
618	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
619	DO jk = 1, jpk
620	DO jj = 1, jpj
621	DO ji = 1, jpi
622	prodatqc%vdmean(ji,jj,jk,1) = 0.0
623	prodatqc%vdmean(ji,jj,jk,2) = 0.0
624	END DO
625	END DO
626	END DO
627
628	ENDIF
629
630	DO jk = 1, jpk
631	DO jj = 1, jpj
632	DO ji = 1, jpi
633	! Increment the temperature field for computing daily mean
634	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
635	& + ptn(ji,jj,jk)
636	! Increment the salinity field for computing daily mean
637	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
638	& + psn(ji,jj,jk)
639	END DO
640	END DO
641	END DO
642
643	! Compute the daily mean at the end of day
644	zdaystp = 1.0 / REAL( kdaystp )
645	IF ( idayend == 0 ) THEN
646	DO jk = 1, jpk
647	DO jj = 1, jpj
648	DO ji = 1, jpi
649	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
650	& * zdaystp
651	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
652	& * zdaystp
653	END DO
654	END DO
655	END DO
656	ENDIF
657
658	! Get the data for interpolation
659	ALLOCATE( &
660	& igrdi(2,2,ipro), &
661	& igrdj(2,2,ipro), &
662	& zglam(2,2,ipro), &
663	& zgphi(2,2,ipro), &
664	& zmask(2,2,kpk,ipro), &
665	& zintt(2,2,kpk,ipro), &
666	& zints(2,2,kpk,ipro), &
667	& zgdept(2,2,kpk,ipro), &
668	& zgdepw(2,2,kpk,ipro) &
669	& )
670
671	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
672	iobs = jobs - prodatqc%nprofup
673	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
674	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
675	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
676	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
677	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
678	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
679	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
680	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
681	END DO
682
683	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
684	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
685	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
686	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
687	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
688	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdept(:,:,:), &
689	& zgdept )
690	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdepw(:,:,:), &
691	& zgdepw )
692
693	! At the end of the day also get interpolated means
694	IF ( idayend == 0 ) THEN
695
696	ALLOCATE( &
697	& zinmt(2,2,kpk,ipro), &
698	& zinms(2,2,kpk,ipro) &
699	& )
700
701	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
702	& prodatqc%vdmean(:,:,:,1), zinmt )
703	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
704	& prodatqc%vdmean(:,:,:,2), zinms )
705
706	ENDIF
707
708	! Return if no observations to process
709	! Has to be done after comm commands to ensure processors
710	! stay in sync
711	IF ( ipro == 0 ) RETURN
712
713	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
714
715	iobs = jobs - prodatqc%nprofup
716
717	IF ( kt /= prodatqc%mstp(jobs) ) THEN
718
719	IF(lwp) THEN
720	WRITE(numout,*)
721	WRITE(numout,*) ' E R R O R : Observation', &
722	& ' time step is not consistent with the', &
723	& ' model time step'
724	WRITE(numout,*) ' ========='
725	WRITE(numout,*)
726	WRITE(numout,*) ' Record = ', jobs, &
727	& ' kt = ', kt, &
728	& ' mstp = ', prodatqc%mstp(jobs), &
729	& ' ntyp = ', prodatqc%ntyp(jobs)
730	ENDIF
731	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
732	ENDIF
733
734	zlam = prodatqc%rlam(jobs)
735	zphi = prodatqc%rphi(jobs)
736
737	! Horizontal weights
738	! Only calculated once, for both T and S.
739	! Masked values are calculated later.
740
741	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
742	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
743
744	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
745	& zglam(:,:,iobs), zgphi(:,:,iobs), &
746	& zmask(:,:,1,iobs), zweig, zmsk_1 )
747
748	ENDIF
749
750	! IF zmsk_1 = 0; then ob is on land
751	IF (zmsk_1(1) < 0.1) THEN
752	WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'
753
754	ELSE
755
756	! Temperature
757
758	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
759
760	zobsk(:) = obfillflt
761
762	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
763
764	IF ( idayend == 0 ) THEN
765
766	! Daily averaged moored buoy (MRB) data
767
768	! vertically interpolate all 4 corners
769	ista = prodatqc%npvsta(jobs,1)
770	iend = prodatqc%npvend(jobs,1)
771	inum_obs = iend - ista + 1
772	ALLOCATE(interp_corner(2,2,inum_obs),v_indic(inum_obs))
773
774	DO iin=1,2
775	DO ijn=1,2
776
777
778
779	IF ( k1dint == 1 ) THEN
780	CALL obs_int_z1d_spl( kpk, &
781	& zinmt(iin,ijn,:,jobs), &
782	& zobs2k, zgdept(iin,ijn,:,jobs), &
783	& zmask(iin,ijn,:,jobs))
784	ENDIF
785
786	CALL obs_level_search(kpk, &
787	& zgdept(iin,ijn,:,jobs), &
788	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
789	& v_indic)
790	CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, &
791	& prodatqc%var(1)%vdep(ista:iend), &
792	& zinmt(iin,ijn,:,jobs), &
793	& zobs2k, interp_corner(iin,ijn,:), &
794	& zgdept(iin,ijn,:,jobs), &
795	& zmask(iin,ijn,:,jobs))
796
797	ENDDO
798	ENDDO
799
800
801	ELSE
802
803	CALL ctl_stop( ' A nonzero' // &
804	& ' number of profile T BUOY data should' // &
805	& ' only occur at the end of a given day' )
806
807	ENDIF
808
809	ELSE
810
811	! Point data
812
813	! vertically interpolate all 4 corners
814	ista = prodatqc%npvsta(jobs,1)
815	iend = prodatqc%npvend(jobs,1)
816	inum_obs = iend - ista + 1
817	ALLOCATE(interp_corner(2,2,inum_obs), v_indic(inum_obs))
818	DO iin=1,2
819	DO ijn=1,2
820
821
822	IF ( k1dint == 1 ) THEN
823	CALL obs_int_z1d_spl( kpk, &
824	& zintt(iin,ijn,:,jobs),&
825	& zobs2k, zgdept(iin,ijn,:,jobs), &
826	& zmask(iin,ijn,:,jobs))
827
828	ENDIF
829
830	CALL obs_level_search(kpk, &
831	& zgdept(iin,ijn,:,jobs),&
832	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
833	& v_indic)
834	CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, &
835	& prodatqc%var(1)%vdep(ista:iend), &
836	& zintt(iin,ijn,:,jobs), &
837	& zobs2k,interp_corner(iin,ijn,:), &
838	& zgdept(iin,ijn,:,jobs), &
839	& zmask(iin,ijn,:,jobs) )
840
841	ENDDO
842	ENDDO
843
844	ENDIF
845
846	!-------------------------------------------------------------
847	! Compute the horizontal interpolation for every profile level
848	!-------------------------------------------------------------
849
850	DO ikn=1,inum_obs
851	iend=ista+ikn-1
852
853	! This code forces the horizontal weights to be
854	! zero IF the observation is below the bottom of the
855	! corners of the interpolation nodes, Or if it is in
856	! the mask. This is important for observations are near
857	! steep bathymetry
858	DO iin=1,2
859	DO ijn=1,2
860
861	depth_loop1: DO ik=kpk,2,-1
862	IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN
863
864	l_zweig(iin,ijn,1) = &
865	& zweig(iin,ijn,1) * &
866	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) &
867	& - prodatqc%var(1)%vdep(iend)),0._wp)
868
869	EXIT depth_loop1
870	ENDIF
871	ENDDO depth_loop1
872
873	ENDDO
874	ENDDO
875
876	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
877	& prodatqc%var(1)%vmod(iend:iend) )
878
879	ENDDO
880
881
882	DEALLOCATE(interp_corner,v_indic)
883
884	ENDIF
885
886
887	! Salinity
888
889	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
890
891	zobsk(:) = obfillflt
892
893	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
894
895	IF ( idayend == 0 ) THEN
896
897	! Daily averaged moored buoy (MRB) data
898
899	! vertically interpolate all 4 corners
900	ista = prodatqc%npvsta(jobs,2)
901	iend = prodatqc%npvend(jobs,2)
902	inum_obs = iend - ista + 1
903	ALLOCATE(interp_corner(2,2,inum_obs),v_indic(inum_obs))
904
905	DO iin=1,2
906	DO ijn=1,2
907
908
909
910	IF ( k1dint == 1 ) THEN
911	CALL obs_int_z1d_spl( kpk, &
912	& zinms(iin,ijn,:,jobs), &
913	& zobs2k, zgdept(iin,ijn,:,jobs), &
914	& zmask(iin,ijn,:,jobs))
915	ENDIF
916
917	CALL obs_level_search(kpk, &
918	& zgdept(iin,ijn,:,jobs), &
919	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
920	& v_indic)
921	CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, &
922	& prodatqc%var(2)%vdep(ista:iend), &
923	& zinms(iin,ijn,:,jobs), &
924	& zobs2k, interp_corner(iin,ijn,:), &
925	& zgdept(iin,ijn,:,jobs), &
926	& zmask(iin,ijn,:,jobs))
927
928	ENDDO
929	ENDDO
930
931
932	ELSE
933
934	CALL ctl_stop( ' A nonzero' // &
935	& ' number of profile T BUOY data should' // &
936	& ' only occur at the end of a given day' )
937
938	ENDIF
939
940	ELSE
941
942	! Point data
943
944	! vertically interpolate all 4 corners
945	ista = prodatqc%npvsta(jobs,2)
946	iend = prodatqc%npvend(jobs,2)
947	inum_obs = iend - ista + 1
948	ALLOCATE(interp_corner(2,2,inum_obs), v_indic(inum_obs))
949
950	DO iin=1,2
951	DO ijn=1,2
952
953
954	IF ( k1dint == 1 ) THEN
955	CALL obs_int_z1d_spl( kpk, &
956	& zints(iin,ijn,:,jobs),&
957	& zobs2k, zgdept(iin,ijn,:,jobs), &
958	& zmask(iin,ijn,:,jobs))
959
960	ENDIF
961
962	CALL obs_level_search(kpk, &
963	& zgdept(iin,ijn,:,jobs),&
964	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
965	& v_indic)
966	CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, &
967	& prodatqc%var(2)%vdep(ista:iend), &
968	& zints(iin,ijn,:,jobs), &
969	& zobs2k,interp_corner(iin,ijn,:), &
970	& zgdept(iin,ijn,:,jobs), &
971	& zmask(iin,ijn,:,jobs) )
972
973	ENDDO
974	ENDDO
975
976	ENDIF
977
978	!-------------------------------------------------------------
979	! Compute the horizontal interpolation for every profile level
980	!-------------------------------------------------------------
981
982	DO ikn=1,inum_obs
983	iend=ista+ikn-1
984
985	! This code forces the horizontal weights to be
986	! zero IF the observation is below the bottom of the
987	! corners of the interpolation nodes, Or if it is in
988	! the mask. This is important for observations are near
989	! steep bathymetry
990	DO iin=1,2
991	DO ijn=1,2
992
993	depth_loop2: DO ik=kpk,2,-1
994	IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN
995
996	l_zweig(iin,ijn,1) = &
997	& zweig(iin,ijn,1) * &
998	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) &
999	& - prodatqc%var(2)%vdep(iend)),0._wp)
1000
1001	EXIT depth_loop2
1002	ENDIF
1003	ENDDO depth_loop2
1004
1005	ENDDO
1006	ENDDO
1007
1008	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
1009	& prodatqc%var(2)%vmod(iend:iend) )
1010
1011	ENDDO
1012
1013
1014	DEALLOCATE(interp_corner,v_indic)
1015
1016	ENDIF
1017
1018	ENDIF
1019
1020	END DO
1021
1022	! Deallocate the data for interpolation
1023	DEALLOCATE( &
1024	& igrdi, &
1025	& igrdj, &
1026	& zglam, &
1027	& zgphi, &
1028	& zmask, &
1029	& zintt, &
1030	& zints &
1031	& )
1032	! At the end of the day also get interpolated means
1033	IF ( idayend == 0 ) THEN
1034	DEALLOCATE( &
1035	& zinmt, &
1036	& zinms &
1037	& )
1038	ENDIF
1039
1040	prodatqc%nprofup = prodatqc%nprofup + ipro
1041
1042	END SUBROUTINE obs_pro_sco_opt
1043
1044	SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
1045	& psshn, psshmask, k2dint )
1046	!!-----------------------------------------------------------------------
1047	!!
1048	!! * ROUTINE obs_sla_opt *
1049	!!
1050	!! ** Purpose : Compute the model counterpart of sea level anomaly
1051	!! data by interpolating from the model grid to the
1052	!! observation point.
1053	!!
1054	!! ** Method : Linearly interpolate to each observation point using
1055	!! the model values at the corners of the surrounding grid box.
1056	!!
1057	!! The now model SSH is first computed at the obs (lon, lat) point.
1058	!!
1059	!! Several horizontal interpolation schemes are available:
1060	!! - distance-weighted (great circle) (k2dint = 0)
1061	!! - distance-weighted (small angle) (k2dint = 1)
1062	!! - bilinear (geographical grid) (k2dint = 2)
1063	!! - bilinear (quadrilateral grid) (k2dint = 3)
1064	!! - polynomial (quadrilateral grid) (k2dint = 4)
1065	!!
1066	!! The sea level anomaly at the observation points is then computed
1067	!! by removing a mean dynamic topography (defined at the obs. point).
1068	!!
1069	!! ** Action :
1070	!!
1071	!! History :
1072	!! ! 07-03 (A. Weaver)
1073	!!-----------------------------------------------------------------------
1074
1075	!! * Modules used
1076	USE obs_surf_def ! Definition of storage space for surface observations
1077
1078	IMPLICIT NONE
1079
1080	!! * Arguments
1081	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening
1082	INTEGER, INTENT(IN) :: kt ! Time step
1083	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1084	INTEGER, INTENT(IN) :: kpj
1085	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1086	! (kit000-1 = restart time)
1087	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1088	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1089	& psshn, & ! Model SSH field
1090	& psshmask ! Land-sea mask
1091
1092	!! * Local declarations
1093	INTEGER :: ji
1094	INTEGER :: jj
1095	INTEGER :: jobs
1096	INTEGER :: inrc
1097	INTEGER :: isla
1098	INTEGER :: iobs
1099	REAL(KIND=wp) :: zlam
1100	REAL(KIND=wp) :: zphi
1101	REAL(KIND=wp) :: zext(1), zobsmask(1)
1102	REAL(kind=wp), DIMENSION(2,2,1) :: &
1103	& zweig
1104	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1105	& zmask, &
1106	& zsshl, &
1107	& zglam, &
1108	& zgphi
1109	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1110	& igrdi, &
1111	& igrdj
1112
1113	!------------------------------------------------------------------------
1114	! Local initialization
1115	!------------------------------------------------------------------------
1116	! ... Record and data counters
1117	inrc = kt - kit000 + 2
1118	isla = sladatqc%nsstp(inrc)
1119
1120	! Get the data for interpolation
1121
1122	ALLOCATE( &
1123	& igrdi(2,2,isla), &
1124	& igrdj(2,2,isla), &
1125	& zglam(2,2,isla), &
1126	& zgphi(2,2,isla), &
1127	& zmask(2,2,isla), &
1128	& zsshl(2,2,isla) &
1129	& )
1130
1131	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1132	iobs = jobs - sladatqc%nsurfup
1133	igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
1134	igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
1135	igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
1136	igrdj(1,2,iobs) = sladatqc%mj(jobs)
1137	igrdi(2,1,iobs) = sladatqc%mi(jobs)
1138	igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
1139	igrdi(2,2,iobs) = sladatqc%mi(jobs)
1140	igrdj(2,2,iobs) = sladatqc%mj(jobs)
1141	END DO
1142
1143	CALL obs_int_comm_2d( 2, 2, isla, &
1144	& igrdi, igrdj, glamt, zglam )
1145	CALL obs_int_comm_2d( 2, 2, isla, &
1146	& igrdi, igrdj, gphit, zgphi )
1147	CALL obs_int_comm_2d( 2, 2, isla, &
1148	& igrdi, igrdj, psshmask, zmask )
1149	CALL obs_int_comm_2d( 2, 2, isla, &
1150	& igrdi, igrdj, psshn, zsshl )
1151
1152	! Loop over observations
1153
1154	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1155
1156	iobs = jobs - sladatqc%nsurfup
1157
1158	IF ( kt /= sladatqc%mstp(jobs) ) THEN
1159
1160	IF(lwp) THEN
1161	WRITE(numout,*)
1162	WRITE(numout,*) ' E R R O R : Observation', &
1163	& ' time step is not consistent with the', &
1164	& ' model time step'
1165	WRITE(numout,*) ' ========='
1166	WRITE(numout,*)
1167	WRITE(numout,*) ' Record = ', jobs, &
1168	& ' kt = ', kt, &
1169	& ' mstp = ', sladatqc%mstp(jobs), &
1170	& ' ntyp = ', sladatqc%ntyp(jobs)
1171	ENDIF
1172	CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
1173
1174	ENDIF
1175
1176	zlam = sladatqc%rlam(jobs)
1177	zphi = sladatqc%rphi(jobs)
1178
1179	! Get weights to interpolate the model SSH to the observation point
1180	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1181	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1182	& zmask(:,:,iobs), zweig, zobsmask )
1183
1184
1185	! Interpolate the model SSH to the observation point
1186	CALL obs_int_h2d( 1, 1, &
1187	& zweig, zsshl(:,:,iobs), zext )
1188
1189	sladatqc%rext(jobs,1) = zext(1)
1190	! ... Remove the MDT at the observation point
1191	sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
1192
1193	END DO
1194
1195	! Deallocate the data for interpolation
1196	DEALLOCATE( &
1197	& igrdi, &
1198	& igrdj, &
1199	& zglam, &
1200	& zgphi, &
1201	& zmask, &
1202	& zsshl &
1203	& )
1204
1205	sladatqc%nsurfup = sladatqc%nsurfup + isla
1206
1207	END SUBROUTINE obs_sla_opt
1208
1209	SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
1210	& psstn, psstmask, k2dint, ld_nightav )
1211	!!-----------------------------------------------------------------------
1212	!!
1213	!! * ROUTINE obs_sst_opt *
1214	!!
1215	!! ** Purpose : Compute the model counterpart of surface temperature
1216	!! data by interpolating from the model grid to the
1217	!! observation point.
1218	!!
1219	!! ** Method : Linearly interpolate to each observation point using
1220	!! the model values at the corners of the surrounding grid box.
1221	!!
1222	!! The now model SST is first computed at the obs (lon, lat) point.
1223	!!
1224	!! Several horizontal interpolation schemes are available:
1225	!! - distance-weighted (great circle) (k2dint = 0)
1226	!! - distance-weighted (small angle) (k2dint = 1)
1227	!! - bilinear (geographical grid) (k2dint = 2)
1228	!! - bilinear (quadrilateral grid) (k2dint = 3)
1229	!! - polynomial (quadrilateral grid) (k2dint = 4)
1230	!!
1231	!!
1232	!! ** Action :
1233	!!
1234	!! History :
1235	!! ! 07-07 (S. Ricci ) : Original
1236	!!
1237	!!-----------------------------------------------------------------------
1238
1239	!! * Modules used
1240	USE obs_surf_def ! Definition of storage space for surface observations
1241	USE sbcdcy
1242
1243	IMPLICIT NONE
1244
1245	!! * Arguments
1246	TYPE(obs_surf), INTENT(INOUT) :: &
1247	& sstdatqc ! Subset of surface data not failing screening
1248	INTEGER, INTENT(IN) :: kt ! Time step
1249	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1250	INTEGER, INTENT(IN) :: kpj
1251	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1252	! (kit000-1 = restart time)
1253	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1254	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1255	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1256	& psstn, & ! Model SST field
1257	& psstmask ! Land-sea mask
1258
1259	!! * Local declarations
1260	INTEGER :: ji
1261	INTEGER :: jj
1262	INTEGER :: jobs
1263	INTEGER :: inrc
1264	INTEGER :: isst
1265	INTEGER :: iobs
1266	INTEGER :: idayend
1267	REAL(KIND=wp) :: zlam
1268	REAL(KIND=wp) :: zphi
1269	REAL(KIND=wp) :: zext(1), zobsmask(1)
1270	REAL(KIND=wp) :: zdaystp
1271	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1272	& icount_sstnight, &
1273	& imask_night
1274	REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1275	& zintmp, &
1276	& zouttmp, &
1277	& zmeanday ! to compute model sst in region of 24h daylight (pole)
1278	REAL(kind=wp), DIMENSION(2,2,1) :: &
1279	& zweig
1280	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1281	& zmask, &
1282	& zsstl, &
1283	& zsstm, &
1284	& zglam, &
1285	& zgphi
1286	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1287	& igrdi, &
1288	& igrdj
1289	LOGICAL, INTENT(IN) :: ld_nightav
1290
1291	!-----------------------------------------------------------------------
1292	! Local initialization
1293	!-----------------------------------------------------------------------
1294	! ... Record and data counters
1295	inrc = kt - kit000 + 2
1296	isst = sstdatqc%nsstp(inrc)
1297
1298	IF ( ld_nightav ) THEN
1299
1300	! Initialize array for night mean
1301
1302	IF ( kt .EQ. 0 ) THEN
1303	ALLOCATE ( icount_sstnight(kpi,kpj) )
1304	ALLOCATE ( imask_night(kpi,kpj) )
1305	ALLOCATE ( zintmp(kpi,kpj) )
1306	ALLOCATE ( zouttmp(kpi,kpj) )
1307	ALLOCATE ( zmeanday(kpi,kpj) )
1308	nday_qsr = -1 ! initialisation flag for nbc_dcy
1309	ENDIF
1310
1311	! Initialize daily mean for first timestep
1312	idayend = MOD( kt - kit000 + 1, kdaystp )
1313
1314	! Added kt == 0 test to catch restart case
1315	IF ( idayend == 1 .OR. kt == 0) THEN
1316	IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
1317	DO jj = 1, jpj
1318	DO ji = 1, jpi
1319	sstdatqc%vdmean(ji,jj) = 0.0
1320	zmeanday(ji,jj) = 0.0
1321	icount_sstnight(ji,jj) = 0
1322	END DO
1323	END DO
1324	ENDIF
1325
1326	zintmp(:,:) = 0.0
1327	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1328	imask_night(:,:) = INT( zouttmp(:,:) )
1329
1330	DO jj = 1, jpj
1331	DO ji = 1, jpi
1332	! Increment the temperature field for computing night mean and counter
1333	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1334	& + psstn(ji,jj)*imask_night(ji,jj)
1335	zmeanday(ji,jj) = zmeanday(ji,jj) + psstn(ji,jj)
1336	icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
1337	END DO
1338	END DO
1339
1340	! Compute the daily mean at the end of day
1341
1342	zdaystp = 1.0 / REAL( kdaystp )
1343
1344	IF ( idayend == 0 ) THEN
1345	DO jj = 1, jpj
1346	DO ji = 1, jpi
1347	! Test if "no night" point
1348	IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
1349	sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1350	& / icount_sstnight(ji,jj)
1351	ELSE
1352	sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1353	ENDIF
1354	END DO
1355	END DO
1356	ENDIF
1357
1358	ENDIF
1359
1360	! Get the data for interpolation
1361
1362	ALLOCATE( &
1363	& igrdi(2,2,isst), &
1364	& igrdj(2,2,isst), &
1365	& zglam(2,2,isst), &
1366	& zgphi(2,2,isst), &
1367	& zmask(2,2,isst), &
1368	& zsstl(2,2,isst) &
1369	& )
1370
1371	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1372	iobs = jobs - sstdatqc%nsurfup
1373	igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
1374	igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
1375	igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
1376	igrdj(1,2,iobs) = sstdatqc%mj(jobs)
1377	igrdi(2,1,iobs) = sstdatqc%mi(jobs)
1378	igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
1379	igrdi(2,2,iobs) = sstdatqc%mi(jobs)
1380	igrdj(2,2,iobs) = sstdatqc%mj(jobs)
1381	END DO
1382
1383	CALL obs_int_comm_2d( 2, 2, isst, &
1384	& igrdi, igrdj, glamt, zglam )
1385	CALL obs_int_comm_2d( 2, 2, isst, &
1386	& igrdi, igrdj, gphit, zgphi )
1387	CALL obs_int_comm_2d( 2, 2, isst, &
1388	& igrdi, igrdj, psstmask, zmask )
1389	CALL obs_int_comm_2d( 2, 2, isst, &
1390	& igrdi, igrdj, psstn, zsstl )
1391
1392	! At the end of the day get interpolated means
1393	IF ( idayend == 0 .AND. ld_nightav ) THEN
1394
1395	ALLOCATE( &
1396	& zsstm(2,2,isst) &
1397	& )
1398
1399	CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
1400	& sstdatqc%vdmean(:,:), zsstm )
1401
1402	ENDIF
1403
1404	! Loop over observations
1405
1406	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1407
1408	iobs = jobs - sstdatqc%nsurfup
1409
1410	IF ( kt /= sstdatqc%mstp(jobs) ) THEN
1411
1412	IF(lwp) THEN
1413	WRITE(numout,*)
1414	WRITE(numout,*) ' E R R O R : Observation', &
1415	& ' time step is not consistent with the', &
1416	& ' model time step'
1417	WRITE(numout,*) ' ========='
1418	WRITE(numout,*)
1419	WRITE(numout,*) ' Record = ', jobs, &
1420	& ' kt = ', kt, &
1421	& ' mstp = ', sstdatqc%mstp(jobs), &
1422	& ' ntyp = ', sstdatqc%ntyp(jobs)
1423	ENDIF
1424	CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
1425
1426	ENDIF
1427
1428	zlam = sstdatqc%rlam(jobs)
1429	zphi = sstdatqc%rphi(jobs)
1430
1431	! Get weights to interpolate the model SST to the observation point
1432	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1433	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1434	& zmask(:,:,iobs), zweig, zobsmask )
1435
1436	! Interpolate the model SST to the observation point
1437
1438	IF ( ld_nightav ) THEN
1439
1440	IF ( idayend == 0 ) THEN
1441	! Daily averaged/diurnal cycle of SST data
1442	CALL obs_int_h2d( 1, 1, &
1443	& zweig, zsstm(:,:,iobs), zext )
1444	ELSE
1445	CALL ctl_stop( ' ld_nightav is set to true: a nonzero' // &
1446	& ' number of night SST data should' // &
1447	& ' only occur at the end of a given day' )
1448	ENDIF
1449
1450	ELSE
1451
1452	CALL obs_int_h2d( 1, 1, &
1453	& zweig, zsstl(:,:,iobs), zext )
1454
1455	ENDIF
1456	sstdatqc%rmod(jobs,1) = zext(1)
1457
1458	END DO
1459
1460	! Deallocate the data for interpolation
1461	DEALLOCATE( &
1462	& igrdi, &
1463	& igrdj, &
1464	& zglam, &
1465	& zgphi, &
1466	& zmask, &
1467	& zsstl &
1468	& )
1469
1470	! At the end of the day also get interpolated means
1471	IF ( idayend == 0 .AND. ld_nightav ) THEN
1472	DEALLOCATE( &
1473	& zsstm &
1474	& )
1475	ENDIF
1476
1477	sstdatqc%nsurfup = sstdatqc%nsurfup + isst
1478
1479	END SUBROUTINE obs_sst_opt
1480
1481	SUBROUTINE obs_sss_opt
1482	!!-----------------------------------------------------------------------
1483	!!
1484	!! * ROUTINE obs_sss_opt *
1485	!!
1486	!! ** Purpose : Compute the model counterpart of sea surface salinity
1487	!! data by interpolating from the model grid to the
1488	!! observation point.
1489	!!
1490	!! ** Method :
1491	!!
1492	!! ** Action :
1493	!!
1494	!! History :
1495	!! ! ??-??
1496	!!-----------------------------------------------------------------------
1497
1498	IMPLICIT NONE
1499
1500	END SUBROUTINE obs_sss_opt
1501
1502	SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
1503	& pseaicen, pseaicemask, k2dint )
1504
1505	!!-----------------------------------------------------------------------
1506	!!
1507	!! * ROUTINE obs_seaice_opt *
1508	!!
1509	!! ** Purpose : Compute the model counterpart of surface temperature
1510	!! data by interpolating from the model grid to the
1511	!! observation point.
1512	!!
1513	!! ** Method : Linearly interpolate to each observation point using
1514	!! the model values at the corners of the surrounding grid box.
1515	!!
1516	!! The now model sea ice is first computed at the obs (lon, lat) point.
1517	!!
1518	!! Several horizontal interpolation schemes are available:
1519	!! - distance-weighted (great circle) (k2dint = 0)
1520	!! - distance-weighted (small angle) (k2dint = 1)
1521	!! - bilinear (geographical grid) (k2dint = 2)
1522	!! - bilinear (quadrilateral grid) (k2dint = 3)
1523	!! - polynomial (quadrilateral grid) (k2dint = 4)
1524	!!
1525	!!
1526	!! ** Action :
1527	!!
1528	!! History :
1529	!! ! 07-07 (S. Ricci ) : Original
1530	!!
1531	!!-----------------------------------------------------------------------
1532
1533	!! * Modules used
1534	USE obs_surf_def ! Definition of storage space for surface observations
1535
1536	IMPLICIT NONE
1537
1538	!! * Arguments
1539	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening
1540	INTEGER, INTENT(IN) :: kt ! Time step
1541	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1542	INTEGER, INTENT(IN) :: kpj
1543	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1544	! (kit000-1 = restart time)
1545	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1546	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1547	& pseaicen, & ! Model sea ice field
1548	& pseaicemask ! Land-sea mask
1549
1550	!! * Local declarations
1551	INTEGER :: ji
1552	INTEGER :: jj
1553	INTEGER :: jobs
1554	INTEGER :: inrc
1555	INTEGER :: iseaice
1556	INTEGER :: iobs
1557
1558	REAL(KIND=wp) :: zlam
1559	REAL(KIND=wp) :: zphi
1560	REAL(KIND=wp) :: zext(1), zobsmask(1)
1561	REAL(kind=wp), DIMENSION(2,2,1) :: &
1562	& zweig
1563	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1564	& zmask, &
1565	& zseaicel, &
1566	& zglam, &
1567	& zgphi
1568	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1569	& igrdi, &
1570	& igrdj
1571
1572	!------------------------------------------------------------------------
1573	! Local initialization
1574	!------------------------------------------------------------------------
1575	! ... Record and data counters
1576	inrc = kt - kit000 + 2
1577	iseaice = seaicedatqc%nsstp(inrc)
1578
1579	! Get the data for interpolation
1580
1581	ALLOCATE( &
1582	& igrdi(2,2,iseaice), &
1583	& igrdj(2,2,iseaice), &
1584	& zglam(2,2,iseaice), &
1585	& zgphi(2,2,iseaice), &
1586	& zmask(2,2,iseaice), &
1587	& zseaicel(2,2,iseaice) &
1588	& )
1589
1590	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1591	iobs = jobs - seaicedatqc%nsurfup
1592	igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1593	igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1594	igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1595	igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1596	igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1597	igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1598	igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1599	igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1600	END DO
1601
1602	CALL obs_int_comm_2d( 2, 2, iseaice, &
1603	& igrdi, igrdj, glamt, zglam )
1604	CALL obs_int_comm_2d( 2, 2, iseaice, &
1605	& igrdi, igrdj, gphit, zgphi )
1606	CALL obs_int_comm_2d( 2, 2, iseaice, &
1607	& igrdi, igrdj, pseaicemask, zmask )
1608	CALL obs_int_comm_2d( 2, 2, iseaice, &
1609	& igrdi, igrdj, pseaicen, zseaicel )
1610
1611	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1612
1613	iobs = jobs - seaicedatqc%nsurfup
1614
1615	IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1616
1617	IF(lwp) THEN
1618	WRITE(numout,*)
1619	WRITE(numout,*) ' E R R O R : Observation', &
1620	& ' time step is not consistent with the', &
1621	& ' model time step'
1622	WRITE(numout,*) ' ========='
1623	WRITE(numout,*)
1624	WRITE(numout,*) ' Record = ', jobs, &
1625	& ' kt = ', kt, &
1626	& ' mstp = ', seaicedatqc%mstp(jobs), &
1627	& ' ntyp = ', seaicedatqc%ntyp(jobs)
1628	ENDIF
1629	CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1630
1631	ENDIF
1632
1633	zlam = seaicedatqc%rlam(jobs)
1634	zphi = seaicedatqc%rphi(jobs)
1635
1636	! Get weights to interpolate the model sea ice to the observation point
1637	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1638	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1639	& zmask(:,:,iobs), zweig, zobsmask )
1640
1641	! ... Interpolate the model sea ice to the observation point
1642	CALL obs_int_h2d( 1, 1, &
1643	& zweig, zseaicel(:,:,iobs), zext )
1644
1645	seaicedatqc%rmod(jobs,1) = zext(1)
1646
1647	END DO
1648
1649	! Deallocate the data for interpolation
1650	DEALLOCATE( &
1651	& igrdi, &
1652	& igrdj, &
1653	& zglam, &
1654	& zgphi, &
1655	& zmask, &
1656	& zseaicel &
1657	& )
1658
1659	seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1660
1661	END SUBROUTINE obs_seaice_opt
1662
1663	SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1664	& pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1665	& ld_dailyav )
1666	!!-----------------------------------------------------------------------
1667	!!
1668	!! * ROUTINE obs_vel_opt *
1669	!!
1670	!! ** Purpose : Compute the model counterpart of velocity profile
1671	!! data by interpolating from the model grid to the
1672	!! observation point.
1673	!!
1674	!! ** Method : Linearly interpolate zonal and meridional components of velocity
1675	!! to each observation point using the model values at the corners of
1676	!! the surrounding grid box. The model velocity components are on a
1677	!! staggered C- grid.
1678	!!
1679	!! For velocity data from the TAO array, the model equivalent is
1680	!! a daily mean velocity field. So, we first compute
1681	!! the mean, then interpolate only at the end of the day.
1682	!!
1683	!! ** Action :
1684	!!
1685	!! History :
1686	!! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles
1687	!! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1688	!!-----------------------------------------------------------------------
1689
1690	!! * Modules used
1691	USE obs_profiles_def ! Definition of storage space for profile obs.
1692
1693	IMPLICIT NONE
1694
1695	!! * Arguments
1696	TYPE(obs_prof), INTENT(INOUT) :: &
1697	& prodatqc ! Subset of profile data not failing screening
1698	INTEGER, INTENT(IN) :: kt ! Time step
1699	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1700	INTEGER, INTENT(IN) :: kpj
1701	INTEGER, INTENT(IN) :: kpk
1702	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1703	! (kit000-1 = restart time)
1704	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
1705	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1706	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1707	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1708	& pun, & ! Model zonal component of velocity
1709	& pvn, & ! Model meridional component of velocity
1710	& pumask, & ! Land-sea mask
1711	& pvmask ! Land-sea mask
1712	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1713	& pgdept ! Model array of depth levels
1714	LOGICAL, INTENT(IN) :: ld_dailyav
1715
1716	!! * Local declarations
1717	INTEGER :: ji
1718	INTEGER :: jj
1719	INTEGER :: jk
1720	INTEGER :: jobs
1721	INTEGER :: inrc
1722	INTEGER :: ipro
1723	INTEGER :: idayend
1724	INTEGER :: ista
1725	INTEGER :: iend
1726	INTEGER :: iobs
1727	INTEGER, DIMENSION(imaxavtypes) :: &
1728	& idailyavtypes
1729	REAL(KIND=wp) :: zlam
1730	REAL(KIND=wp) :: zphi
1731	REAL(KIND=wp) :: zdaystp
1732	REAL(KIND=wp), DIMENSION(kpk) :: &
1733	& zobsmasku, &
1734	& zobsmaskv, &
1735	& zobsmask, &
1736	& zobsk, &
1737	& zobs2k
1738	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1739	& zweigu,zweigv
1740	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1741	& zumask, zvmask, &
1742	& zintu, &
1743	& zintv, &
1744	& zinmu, &
1745	& zinmv
1746	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1747	& zglamu, zglamv, &
1748	& zgphiu, zgphiv
1749	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1750	& igrdiu, &
1751	& igrdju, &
1752	& igrdiv, &
1753	& igrdjv
1754
1755	!------------------------------------------------------------------------
1756	! Local initialization
1757	!------------------------------------------------------------------------
1758	! ... Record and data counters
1759	inrc = kt - kit000 + 2
1760	ipro = prodatqc%npstp(inrc)
1761
1762	! Initialize daily mean for first timestep
1763	idayend = MOD( kt - kit000 + 1, kdaystp )
1764
1765	! Added kt == 0 test to catch restart case
1766	IF ( idayend == 1 .OR. kt == 0) THEN
1767	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1768	prodatqc%vdmean(:,:,:,1) = 0.0
1769	prodatqc%vdmean(:,:,:,2) = 0.0
1770	ENDIF
1771
1772	! Increment the zonal velocity field for computing daily mean
1773	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1774	! Increment the meridional velocity field for computing daily mean
1775	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1776
1777	! Compute the daily mean at the end of day
1778	zdaystp = 1.0 / REAL( kdaystp )
1779	IF ( idayend == 0 ) THEN
1780	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1781	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1782	ENDIF
1783
1784	! Get the data for interpolation
1785	ALLOCATE( &
1786	& igrdiu(2,2,ipro), &
1787	& igrdju(2,2,ipro), &
1788	& igrdiv(2,2,ipro), &
1789	& igrdjv(2,2,ipro), &
1790	& zglamu(2,2,ipro), zglamv(2,2,ipro), &
1791	& zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1792	& zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1793	& zintu(2,2,kpk,ipro), &
1794	& zintv(2,2,kpk,ipro) &
1795	& )
1796
1797	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1798	iobs = jobs - prodatqc%nprofup
1799	igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1800	igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1801	igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1802	igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1803	igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1804	igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1805	igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1806	igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1807	igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1808	igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1809	igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1810	igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1811	igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1812	igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1813	igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1814	igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1815	END DO
1816
1817	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1818	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1819	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1820	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1821
1822	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1823	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1824	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1825	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1826
1827	! At the end of the day also get interpolated means
1828	IF ( idayend == 0 ) THEN
1829
1830	ALLOCATE( &
1831	& zinmu(2,2,kpk,ipro), &
1832	& zinmv(2,2,kpk,ipro) &
1833	& )
1834
1835	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1836	& prodatqc%vdmean(:,:,:,1), zinmu )
1837	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1838	& prodatqc%vdmean(:,:,:,2), zinmv )
1839
1840	ENDIF
1841
1842	! loop over observations
1843
1844	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1845
1846	iobs = jobs - prodatqc%nprofup
1847
1848	IF ( kt /= prodatqc%mstp(jobs) ) THEN
1849
1850	IF(lwp) THEN
1851	WRITE(numout,*)
1852	WRITE(numout,*) ' E R R O R : Observation', &
1853	& ' time step is not consistent with the', &
1854	& ' model time step'
1855	WRITE(numout,*) ' ========='
1856	WRITE(numout,*)
1857	WRITE(numout,*) ' Record = ', jobs, &
1858	& ' kt = ', kt, &
1859	& ' mstp = ', prodatqc%mstp(jobs), &
1860	& ' ntyp = ', prodatqc%ntyp(jobs)
1861	ENDIF
1862	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1863	ENDIF
1864
1865	zlam = prodatqc%rlam(jobs)
1866	zphi = prodatqc%rphi(jobs)
1867
1868	! Initialize observation masks
1869
1870	zobsmasku(:) = 0.0
1871	zobsmaskv(:) = 0.0
1872
1873	! Horizontal weights and vertical mask
1874
1875	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1876
1877	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1878	& zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1879	& zumask(:,:,:,iobs), zweigu, zobsmasku )
1880
1881	ENDIF
1882
1883
1884	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1885
1886	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
1887	& zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1888	& zvmask(:,:,:,iobs), zweigv, zobsmaskv )
1889
1890	ENDIF
1891
1892	! Ensure that the vertical mask on u and v are consistent.
1893
1894	zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1895
1896	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1897
1898	zobsk(:) = obfillflt
1899
1900	IF ( ld_dailyav ) THEN
1901
1902	IF ( idayend == 0 ) THEN
1903
1904	! Daily averaged data
1905
1906	CALL obs_int_h2d( kpk, kpk, &
1907	& zweigu, zinmu(:,:,:,iobs), zobsk )
1908
1909
1910	ELSE
1911
1912	CALL ctl_stop( ' A nonzero' // &
1913	& ' number of U profile data should' // &
1914	& ' only occur at the end of a given day' )
1915
1916	ENDIF
1917
1918	ELSE
1919
1920	! Point data
1921
1922	CALL obs_int_h2d( kpk, kpk, &
1923	& zweigu, zintu(:,:,:,iobs), zobsk )
1924
1925	ENDIF
1926
1927	!-------------------------------------------------------------
1928	! Compute vertical second-derivative of the interpolating
1929	! polynomial at obs points
1930	!-------------------------------------------------------------
1931
1932	IF ( k1dint == 1 ) THEN
1933	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1934	& pgdept, zobsmask )
1935	ENDIF
1936
1937	!-----------------------------------------------------------------
1938	! Vertical interpolation to the observation point
1939	!-----------------------------------------------------------------
1940	ista = prodatqc%npvsta(jobs,1)
1941	iend = prodatqc%npvend(jobs,1)
1942	CALL obs_int_z1d( kpk, &
1943	& prodatqc%var(1)%mvk(ista:iend), &
1944	& k1dint, iend - ista + 1, &
1945	& prodatqc%var(1)%vdep(ista:iend), &
1946	& zobsk, zobs2k, &
1947	& prodatqc%var(1)%vmod(ista:iend), &
1948	& pgdept, zobsmask )
1949
1950	ENDIF
1951
1952	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1953
1954	zobsk(:) = obfillflt
1955
1956	IF ( ld_dailyav ) THEN
1957
1958	IF ( idayend == 0 ) THEN
1959
1960	! Daily averaged data
1961
1962	CALL obs_int_h2d( kpk, kpk, &
1963	& zweigv, zinmv(:,:,:,iobs), zobsk )
1964
1965	ELSE
1966
1967	CALL ctl_stop( ' A nonzero' // &
1968	& ' number of V profile data should' // &
1969	& ' only occur at the end of a given day' )
1970
1971	ENDIF
1972
1973	ELSE
1974
1975	! Point data
1976
1977	CALL obs_int_h2d( kpk, kpk, &
1978	& zweigv, zintv(:,:,:,iobs), zobsk )
1979
1980	ENDIF
1981
1982
1983	!-------------------------------------------------------------
1984	! Compute vertical second-derivative of the interpolating
1985	! polynomial at obs points
1986	!-------------------------------------------------------------
1987
1988	IF ( k1dint == 1 ) THEN
1989	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1990	& pgdept, zobsmask )
1991	ENDIF
1992
1993	!----------------------------------------------------------------
1994	! Vertical interpolation to the observation point
1995	!----------------------------------------------------------------
1996	ista = prodatqc%npvsta(jobs,2)
1997	iend = prodatqc%npvend(jobs,2)
1998	CALL obs_int_z1d( kpk, &
1999	& prodatqc%var(2)%mvk(ista:iend),&
2000	& k1dint, iend - ista + 1, &
2001	& prodatqc%var(2)%vdep(ista:iend),&
2002	& zobsk, zobs2k, &
2003	& prodatqc%var(2)%vmod(ista:iend),&
2004	& pgdept, zobsmask )
2005
2006	ENDIF
2007
2008	END DO
2009
2010	! Deallocate the data for interpolation
2011	DEALLOCATE( &
2012	& igrdiu, &
2013	& igrdju, &
2014	& igrdiv, &
2015	& igrdjv, &
2016	& zglamu, zglamv, &
2017	& zgphiu, zgphiv, &
2018	& zumask, zvmask, &
2019	& zintu, &
2020	& zintv &
2021	& )
2022	! At the end of the day also get interpolated means
2023	IF ( idayend == 0 ) THEN
2024	DEALLOCATE( &
2025	& zinmu, &
2026	& zinmv &
2027	& )
2028	ENDIF
2029
2030	prodatqc%nprofup = prodatqc%nprofup + ipro
2031
2032	END SUBROUTINE obs_vel_opt
2033
2034	END MODULE obs_oper
2035

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

Download in other formats: