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 @ 5950

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

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

Download in other formats: