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

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

obs_oper.F90 in branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

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

Last change on this file since 5998 was 5998, checked in by timgraham, 8 years ago
Merge in changes from OBS simplification branch (branches/2015/dev_r5072_UKMO2_OBS_simplification)
Property svn:keywords set to `Id`
File size: 50.8 KB

Line
1	MODULE obs_oper
2	!!======================================================================
3	!! * MODULE obs_oper *
4	!! Observation diagnostics: Observation operators for various observation
5	!! types
6	!!======================================================================
7
8	!!----------------------------------------------------------------------
9	!! obs_prof_opt : Compute the model counterpart of profile data
10	!! obs_surf_opt : Compute the model counterpart of surface data
11	!! obs_pro_sco_opt: Compute the model counterpart of temperature and
12	!! salinity observations from profiles in generalised
13	!! vertical coordinates
14	!!----------------------------------------------------------------------
15
16	!! * Modules used
17	USE par_kind, ONLY : & ! Precision variables
18	& wp
19	USE in_out_manager ! I/O manager
20	USE obs_inter_sup ! Interpolation support
21	USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt
22	& obs_int_h2d, &
23	& obs_int_h2d_init
24	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt
25	& obs_int_z1d, &
26	& obs_int_z1d_spl
27	USE obs_const, ONLY : &
28	& obfillflt ! Fillvalue
29	USE dom_oce, ONLY : &
30	& glamt, glamu, glamv, &
31	& gphit, gphiu, gphiv, &
32	#if defined key_vvl
33	& gdept_n
34	#else
35	& gdept_0
36	#endif
37	USE lib_mpp, ONLY : &
38	& ctl_warn, ctl_stop
39	USE obs_grid, ONLY : &
40	& obs_level_search
41	USE sbcdcy, ONLY : & ! For calculation of where it is night-time
42	& sbc_dcy, nday_qsr
43
44	IMPLICIT NONE
45
46	!! * Routine accessibility
47	PRIVATE
48
49	PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs
50	& obs_pro_sco_opt, & ! Compute the model counterpart of profile observations
51	& obs_surf_opt ! Compute the model counterpart of surface obs
52
53	INTEGER, PARAMETER, PUBLIC :: &
54	& imaxavtypes = 20 ! Max number of daily avgd obs types
55
56	!!----------------------------------------------------------------------
57	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
58	!! $Id$
59	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
60	!!----------------------------------------------------------------------
61
62	!! * Substitutions
63	# include "domzgr_substitute.h90"
64	CONTAINS
65
66	SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
67	& kit000, kdaystp, &
68	& pvar1, pvar2, pgdept, pmask1, pmask2, &
69	& plam1, plam2, pphi1, pphi2, &
70	& k1dint, k2dint, kdailyavtypes )
71
72	!!-----------------------------------------------------------------------
73	!!
74	!! * ROUTINE obs_pro_opt *
75	!!
76	!! ** Purpose : Compute the model counterpart of profiles
77	!! data by interpolating from the model grid to the
78	!! observation point.
79	!!
80	!! ** Method : Linearly interpolate to each observation point using
81	!! the model values at the corners of the surrounding grid box.
82	!!
83	!! First, a vertical profile of horizontally interpolated model
84	!! now values is computed at the obs (lon, lat) point.
85	!! Several horizontal interpolation schemes are available:
86	!! - distance-weighted (great circle) (k2dint = 0)
87	!! - distance-weighted (small angle) (k2dint = 1)
88	!! - bilinear (geographical grid) (k2dint = 2)
89	!! - bilinear (quadrilateral grid) (k2dint = 3)
90	!! - polynomial (quadrilateral grid) (k2dint = 4)
91	!!
92	!! Next, the vertical profile is interpolated to the
93	!! data depth points. Two vertical interpolation schemes are
94	!! available:
95	!! - linear (k1dint = 0)
96	!! - Cubic spline (k1dint = 1)
97	!!
98	!! For the cubic spline the 2nd derivative of the interpolating
99	!! polynomial is computed before entering the vertical interpolation
100	!! routine.
101	!!
102	!! If the logical is switched on, the model equivalent is
103	!! a daily mean model temperature field. So, we first compute
104	!! the mean, then interpolate only at the end of the day.
105	!!
106	!! Note: in situ temperature observations must be converted
107	!! to potential temperature (the model variable) prior to
108	!! assimilation.
109	!!
110	!! ** Action :
111	!!
112	!! History :
113	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
114	!! ! 06-03 (G. Smith) NEMOVAR migration
115	!! ! 06-10 (A. Weaver) Cleanup
116	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
117	!! ! 07-03 (K. Mogensen) General handling of profiles
118	!! ! 15-02 (M. Martin) Combined routine for all profile types
119	!!-----------------------------------------------------------------------
120
121	!! * Modules used
122	USE obs_profiles_def ! Definition of storage space for profile obs.
123
124	IMPLICIT NONE
125
126	!! * Arguments
127	TYPE(obs_prof), INTENT(INOUT) :: &
128	& prodatqc ! Subset of profile data passing QC
129	INTEGER, INTENT(IN) :: kt ! Time step
130	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
131	INTEGER, INTENT(IN) :: kpj
132	INTEGER, INTENT(IN) :: kpk
133	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
134	! (kit000-1 = restart time)
135	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
136	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
137	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
138	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
139	& pvar1, & ! Model field 1
140	& pvar2, & ! Model field 2
141	& pmask1, & ! Land-sea mask 1
142	& pmask2 ! Land-sea mask 2
143	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
144	& plam1, & ! Model longitudes for variable 1
145	& plam2, & ! Model longitudes for variable 2
146	& pphi1, & ! Model latitudes for variable 1
147	& pphi2 ! Model latitudes for variable 2
148	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
149	& pgdept ! Model array of depth levels
150	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
151	& kdailyavtypes ! Types for daily averages
152
153	!! * Local declarations
154	INTEGER :: ji
155	INTEGER :: jj
156	INTEGER :: jk
157	INTEGER :: jobs
158	INTEGER :: inrc
159	INTEGER :: ipro
160	INTEGER :: idayend
161	INTEGER :: ista
162	INTEGER :: iend
163	INTEGER :: iobs
164	INTEGER, DIMENSION(imaxavtypes) :: &
165	& idailyavtypes
166	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
167	& igrdi1, &
168	& igrdi2, &
169	& igrdj1, &
170	& igrdj2
171	REAL(KIND=wp) :: zlam
172	REAL(KIND=wp) :: zphi
173	REAL(KIND=wp) :: zdaystp
174	REAL(KIND=wp), DIMENSION(kpk) :: &
175	& zobsmask1, &
176	& zobsmask2, &
177	& zobsk, &
178	& zobs2k
179	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
180	& zweig1, &
181	& zweig2
182	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
183	& zmask1, &
184	& zmask2, &
185	& zint1, &
186	& zint2, &
187	& zinm1, &
188	& zinm2
189	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
190	& zglam1, &
191	& zglam2, &
192	& zgphi1, &
193	& zgphi2
194	LOGICAL :: ld_dailyav
195
196	!------------------------------------------------------------------------
197	! Local initialization
198	!------------------------------------------------------------------------
199	! Record and data counters
200	inrc = kt - kit000 + 2
201	ipro = prodatqc%npstp(inrc)
202
203	! Daily average types
204	ld_dailyav = .FALSE.
205	IF ( PRESENT(kdailyavtypes) ) THEN
206	idailyavtypes(:) = kdailyavtypes(:)
207	IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
208	ELSE
209	idailyavtypes(:) = -1
210	ENDIF
211
212	! Daily means are calculated for values over timesteps:
213	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
214	idayend = MOD( kt - kit000 + 1, kdaystp )
215
216	IF ( ld_dailyav ) THEN
217
218	! Initialize daily mean for first timestep of the day
219	IF ( idayend == 1 .OR. kt == 0 ) THEN
220	DO jk = 1, jpk
221	DO jj = 1, jpj
222	DO ji = 1, jpi
223	prodatqc%vdmean(ji,jj,jk,1) = 0.0
224	prodatqc%vdmean(ji,jj,jk,2) = 0.0
225	END DO
226	END DO
227	END DO
228	ENDIF
229
230	DO jk = 1, jpk
231	DO jj = 1, jpj
232	DO ji = 1, jpi
233	! Increment field 1 for computing daily mean
234	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
235	& + pvar1(ji,jj,jk)
236	! Increment field 2 for computing daily mean
237	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
238	& + pvar2(ji,jj,jk)
239	END DO
240	END DO
241	END DO
242
243	! Compute the daily mean at the end of day
244	zdaystp = 1.0 / REAL( kdaystp )
245	IF ( idayend == 0 ) THEN
246	IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
247	CALL FLUSH(numout)
248	DO jk = 1, jpk
249	DO jj = 1, jpj
250	DO ji = 1, jpi
251	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
252	& * zdaystp
253	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
254	& * zdaystp
255	END DO
256	END DO
257	END DO
258	ENDIF
259
260	ENDIF
261
262	! Get the data for interpolation
263	ALLOCATE( &
264	& igrdi1(2,2,ipro), &
265	& igrdi2(2,2,ipro), &
266	& igrdj1(2,2,ipro), &
267	& igrdj2(2,2,ipro), &
268	& zglam1(2,2,ipro), &
269	& zglam2(2,2,ipro), &
270	& zgphi1(2,2,ipro), &
271	& zgphi2(2,2,ipro), &
272	& zmask1(2,2,kpk,ipro), &
273	& zmask2(2,2,kpk,ipro), &
274	& zint1(2,2,kpk,ipro), &
275	& zint2(2,2,kpk,ipro) &
276	& )
277
278	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
279	iobs = jobs - prodatqc%nprofup
280	igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
281	igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
282	igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
283	igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
284	igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
285	igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
286	igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
287	igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
288	igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
289	igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
290	igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
291	igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
292	igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
293	igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
294	igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
295	igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
296	END DO
297
298	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
299	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
300	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
301	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 )
302
303	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
304	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
305	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
306	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 )
307
308	! At the end of the day also get interpolated means
309	IF ( ld_dailyav .AND. idayend == 0 ) THEN
310
311	ALLOCATE( &
312	& zinm1(2,2,kpk,ipro), &
313	& zinm2(2,2,kpk,ipro) &
314	& )
315
316	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
317	& prodatqc%vdmean(:,:,:,1), zinm1 )
318	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
319	& prodatqc%vdmean(:,:,:,2), zinm2 )
320
321	ENDIF
322
323	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
324
325	iobs = jobs - prodatqc%nprofup
326
327	IF ( kt /= prodatqc%mstp(jobs) ) THEN
328
329	IF(lwp) THEN
330	WRITE(numout,*)
331	WRITE(numout,*) ' E R R O R : Observation', &
332	& ' time step is not consistent with the', &
333	& ' model time step'
334	WRITE(numout,*) ' ========='
335	WRITE(numout,*)
336	WRITE(numout,*) ' Record = ', jobs, &
337	& ' kt = ', kt, &
338	& ' mstp = ', prodatqc%mstp(jobs), &
339	& ' ntyp = ', prodatqc%ntyp(jobs)
340	ENDIF
341	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
342	ENDIF
343
344	zlam = prodatqc%rlam(jobs)
345	zphi = prodatqc%rphi(jobs)
346
347	! Horizontal weights and vertical mask
348
349	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
350
351	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
352	& zglam1(:,:,iobs), zgphi1(:,:,iobs), &
353	& zmask1(:,:,:,iobs), zweig1, zobsmask1 )
354
355	ENDIF
356
357	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
358
359	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
360	& zglam2(:,:,iobs), zgphi2(:,:,iobs), &
361	& zmask2(:,:,:,iobs), zweig2, zobsmask2 )
362
363	ENDIF
364
365	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
366
367	zobsk(:) = obfillflt
368
369	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
370
371	IF ( idayend == 0 ) THEN
372	! Daily averaged data
373	CALL obs_int_h2d( kpk, kpk, &
374	& zweig1, zinm1(:,:,:,iobs), zobsk )
375
376	ENDIF
377
378	ELSE
379
380	! Point data
381	CALL obs_int_h2d( kpk, kpk, &
382	& zweig1, zint1(:,:,:,iobs), zobsk )
383
384	ENDIF
385
386	!-------------------------------------------------------------
387	! Compute vertical second-derivative of the interpolating
388	! polynomial at obs points
389	!-------------------------------------------------------------
390
391	IF ( k1dint == 1 ) THEN
392	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
393	& pgdept, zobsmask1 )
394	ENDIF
395
396	!-----------------------------------------------------------------
397	! Vertical interpolation to the observation point
398	!-----------------------------------------------------------------
399	ista = prodatqc%npvsta(jobs,1)
400	iend = prodatqc%npvend(jobs,1)
401	CALL obs_int_z1d( kpk, &
402	& prodatqc%var(1)%mvk(ista:iend), &
403	& k1dint, iend - ista + 1, &
404	& prodatqc%var(1)%vdep(ista:iend), &
405	& zobsk, zobs2k, &
406	& prodatqc%var(1)%vmod(ista:iend), &
407	& pgdept, zobsmask1 )
408
409	ENDIF
410
411	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
412
413	zobsk(:) = obfillflt
414
415	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
416
417	IF ( idayend == 0 ) THEN
418
419	! Daily averaged data
420	CALL obs_int_h2d( kpk, kpk, &
421	& zweig2, zinm2(:,:,:,iobs), zobsk )
422
423	ENDIF
424
425	ELSE
426
427	! Point data
428	CALL obs_int_h2d( kpk, kpk, &
429	& zweig2, zint2(:,:,:,iobs), zobsk )
430
431	ENDIF
432
433
434	!-------------------------------------------------------------
435	! Compute vertical second-derivative of the interpolating
436	! polynomial at obs points
437	!-------------------------------------------------------------
438
439	IF ( k1dint == 1 ) THEN
440	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
441	& pgdept, zobsmask2 )
442	ENDIF
443
444	!----------------------------------------------------------------
445	! Vertical interpolation to the observation point
446	!----------------------------------------------------------------
447	ista = prodatqc%npvsta(jobs,2)
448	iend = prodatqc%npvend(jobs,2)
449	CALL obs_int_z1d( kpk, &
450	& prodatqc%var(2)%mvk(ista:iend),&
451	& k1dint, iend - ista + 1, &
452	& prodatqc%var(2)%vdep(ista:iend),&
453	& zobsk, zobs2k, &
454	& prodatqc%var(2)%vmod(ista:iend),&
455	& pgdept, zobsmask2 )
456
457	ENDIF
458
459	END DO
460
461	! Deallocate the data for interpolation
462	DEALLOCATE( &
463	& igrdi1, &
464	& igrdi2, &
465	& igrdj1, &
466	& igrdj2, &
467	& zglam1, &
468	& zglam2, &
469	& zgphi1, &
470	& zgphi2, &
471	& zmask1, &
472	& zmask2, &
473	& zint1, &
474	& zint2 &
475	& )
476
477	! At the end of the day also get interpolated means
478	IF ( ld_dailyav .AND. idayend == 0 ) THEN
479	DEALLOCATE( &
480	& zinm1, &
481	& zinm2 &
482	& )
483	ENDIF
484
485	prodatqc%nprofup = prodatqc%nprofup + ipro
486
487	END SUBROUTINE obs_prof_opt
488
489	SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
490	& ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, &
491	& kdailyavtypes )
492	!!-----------------------------------------------------------------------
493	!!
494	!! * ROUTINE obs_pro_opt *
495	!!
496	!! ** Purpose : Compute the model counterpart of profiles
497	!! data by interpolating from the model grid to the
498	!! observation point. Generalised vertical coordinate version
499	!!
500	!! ** Method : Linearly interpolate to each observation point using
501	!! the model values at the corners of the surrounding grid box.
502	!!
503	!! First, model values on the model grid are interpolated vertically to the
504	!! Depths of the profile observations. Two vertical interpolation schemes are
505	!! available:
506	!! - linear (k1dint = 0)
507	!! - Cubic spline (k1dint = 1)
508	!!
509	!!
510	!! Secondly the interpolated values are interpolated horizontally to the
511	!! obs (lon, lat) point.
512	!! Several horizontal interpolation schemes are available:
513	!! - distance-weighted (great circle) (k2dint = 0)
514	!! - distance-weighted (small angle) (k2dint = 1)
515	!! - bilinear (geographical grid) (k2dint = 2)
516	!! - bilinear (quadrilateral grid) (k2dint = 3)
517	!! - polynomial (quadrilateral grid) (k2dint = 4)
518	!!
519	!! For the cubic spline the 2nd derivative of the interpolating
520	!! polynomial is computed before entering the vertical interpolation
521	!! routine.
522	!!
523	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
524	!! a daily mean model temperature field. So, we first compute
525	!! the mean, then interpolate only at the end of the day.
526	!!
527	!! This is the procedure to be used with generalised vertical model
528	!! coordinates (ie s-coordinates. It is ~4x slower than the equivalent
529	!! horizontal then vertical interpolation algorithm, but can deal with situations
530	!! where the model levels are not flat.
531	!! ONLY PERFORMED if ln_sco=.TRUE.
532	!!
533	!! Note: the in situ temperature observations must be converted
534	!! to potential temperature (the model variable) prior to
535	!! assimilation.
536	!!??????????????????????????????????????????????????????????????
537	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
538	!!??????????????????????????????????????????????????????????????
539	!!
540	!! ** Action :
541	!!
542	!! History :
543	!! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised
544	!! vertical coordinates
545	!!-----------------------------------------------------------------------
546
547	!! * Modules used
548	USE obs_profiles_def ! Definition of storage space for profile obs.
549	USE dom_oce, ONLY : &
550	#if defined key_vvl
551	& gdepw_n
552	#else
553	& gdepw_0
554	#endif
555
556	IMPLICIT NONE
557
558	!! * Arguments
559	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
560	INTEGER, INTENT(IN) :: kt ! Time step
561	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
562	INTEGER, INTENT(IN) :: kpj
563	INTEGER, INTENT(IN) :: kpk
564	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
565	! (kit000-1 = restart time)
566	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
567	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
568	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
569	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
570	& ptn, & ! Model temperature field
571	& psn, & ! Model salinity field
572	& ptmask ! Land-sea mask
573	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,jpj,kpk) :: &
574	& pgdept, & ! Model array of depth T levels
575	& pgdepw ! Model array of depth W levels
576	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
577	& kdailyavtypes ! Types for daily averages
578
579	!! * Local declarations
580	INTEGER :: ji
581	INTEGER :: jj
582	INTEGER :: jk
583	INTEGER :: iico, ijco
584	INTEGER :: jobs
585	INTEGER :: inrc
586	INTEGER :: ipro
587	INTEGER :: idayend
588	INTEGER :: ista
589	INTEGER :: iend
590	INTEGER :: iobs
591	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
592	INTEGER, DIMENSION(imaxavtypes) :: &
593	& idailyavtypes
594	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
595	& igrdi, &
596	& igrdj
597	INTEGER :: &
598	& inum_obs
599	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
600	REAL(KIND=wp) :: zlam
601	REAL(KIND=wp) :: zphi
602	REAL(KIND=wp) :: zdaystp
603	REAL(KIND=wp), DIMENSION(kpk) :: &
604	& zobsmask, &
605	& zobsk, &
606	& zobs2k
607	REAL(KIND=wp), DIMENSION(2,2,1) :: &
608	& zweig, &
609	& l_zweig
610	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
611	& zmask, &
612	& zintt, &
613	& zints, &
614	& zinmt, &
615	& zgdept,&
616	& zgdepw,&
617	& zinms
618	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
619	& zglam, &
620	& zgphi
621	REAL(KIND=wp), DIMENSION(1) :: zmsk_1
622	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
623
624	!------------------------------------------------------------------------
625	! Local initialization
626	!------------------------------------------------------------------------
627	! ... Record and data counters
628	inrc = kt - kit000 + 2
629	ipro = prodatqc%npstp(inrc)
630
631	! Daily average types
632	IF ( PRESENT(kdailyavtypes) ) THEN
633	idailyavtypes(:) = kdailyavtypes(:)
634	ELSE
635	idailyavtypes(:) = -1
636	ENDIF
637
638	! Initialize daily mean for first time-step
639	idayend = MOD( kt - kit000 + 1, kdaystp )
640
641	! Added kt == 0 test to catch restart case
642	IF ( idayend == 1 .OR. kt == 0) THEN
643
644	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
645	DO jk = 1, jpk
646	DO jj = 1, jpj
647	DO ji = 1, jpi
648	prodatqc%vdmean(ji,jj,jk,1) = 0.0
649	prodatqc%vdmean(ji,jj,jk,2) = 0.0
650	END DO
651	END DO
652	END DO
653
654	ENDIF
655
656	DO jk = 1, jpk
657	DO jj = 1, jpj
658	DO ji = 1, jpi
659	! Increment the temperature field for computing daily mean
660	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
661	& + ptn(ji,jj,jk)
662	! Increment the salinity field for computing daily mean
663	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
664	& + psn(ji,jj,jk)
665	END DO
666	END DO
667	END DO
668
669	! Compute the daily mean at the end of day
670	zdaystp = 1.0 / REAL( kdaystp )
671	IF ( idayend == 0 ) THEN
672	DO jk = 1, jpk
673	DO jj = 1, jpj
674	DO ji = 1, jpi
675	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
676	& * zdaystp
677	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
678	& * zdaystp
679	END DO
680	END DO
681	END DO
682	ENDIF
683
684	! Get the data for interpolation
685	ALLOCATE( &
686	& igrdi(2,2,ipro), &
687	& igrdj(2,2,ipro), &
688	& zglam(2,2,ipro), &
689	& zgphi(2,2,ipro), &
690	& zmask(2,2,kpk,ipro), &
691	& zintt(2,2,kpk,ipro), &
692	& zints(2,2,kpk,ipro), &
693	& zgdept(2,2,kpk,ipro), &
694	& zgdepw(2,2,kpk,ipro) &
695	& )
696
697	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
698	iobs = jobs - prodatqc%nprofup
699	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
700	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
701	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
702	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
703	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
704	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
705	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
706	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
707	END DO
708
709	! Initialise depth arrays
710	zgdept = 0.0
711	zgdepw = 0.0
712
713	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
714	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
715	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
716	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
717	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
718	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdept(:,:,:), &
719	& zgdept )
720	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdepw(:,:,:), &
721	& zgdepw )
722
723	! At the end of the day also get interpolated means
724	IF ( idayend == 0 ) THEN
725
726	ALLOCATE( &
727	& zinmt(2,2,kpk,ipro), &
728	& zinms(2,2,kpk,ipro) &
729	& )
730
731	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
732	& prodatqc%vdmean(:,:,:,1), zinmt )
733	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
734	& prodatqc%vdmean(:,:,:,2), zinms )
735
736	ENDIF
737
738	! Return if no observations to process
739	! Has to be done after comm commands to ensure processors
740	! stay in sync
741	IF ( ipro == 0 ) RETURN
742
743	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
744
745	iobs = jobs - prodatqc%nprofup
746
747	IF ( kt /= prodatqc%mstp(jobs) ) THEN
748
749	IF(lwp) THEN
750	WRITE(numout,*)
751	WRITE(numout,*) ' E R R O R : Observation', &
752	& ' time step is not consistent with the', &
753	& ' model time step'
754	WRITE(numout,*) ' ========='
755	WRITE(numout,*)
756	WRITE(numout,*) ' Record = ', jobs, &
757	& ' kt = ', kt, &
758	& ' mstp = ', prodatqc%mstp(jobs), &
759	& ' ntyp = ', prodatqc%ntyp(jobs)
760	ENDIF
761	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
762	ENDIF
763
764	zlam = prodatqc%rlam(jobs)
765	zphi = prodatqc%rphi(jobs)
766
767	! Horizontal weights
768	! Only calculated once, for both T and S.
769	! Masked values are calculated later.
770
771	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
772	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
773
774	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
775	& zglam(:,:,iobs), zgphi(:,:,iobs), &
776	& zmask(:,:,1,iobs), zweig, zmsk_1 )
777
778	ENDIF
779
780	! IF zmsk_1 = 0; then ob is on land
781	IF (zmsk_1(1) < 0.1) THEN
782	WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'
783
784	ELSE
785
786	! Temperature
787
788	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
789
790	zobsk(:) = obfillflt
791
792	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
793
794	IF ( idayend == 0 ) THEN
795
796	! Daily averaged moored buoy (MRB) data
797
798	! vertically interpolate all 4 corners
799	ista = prodatqc%npvsta(jobs,1)
800	iend = prodatqc%npvend(jobs,1)
801	inum_obs = iend - ista + 1
802	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
803
804	DO iin=1,2
805	DO ijn=1,2
806
807
808
809	IF ( k1dint == 1 ) THEN
810	CALL obs_int_z1d_spl( kpk, &
811	& zinmt(iin,ijn,:,iobs), &
812	& zobs2k, zgdept(iin,ijn,:,iobs), &
813	& zmask(iin,ijn,:,iobs))
814	ENDIF
815
816	CALL obs_level_search(kpk, &
817	& zgdept(iin,ijn,:,iobs), &
818	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
819	& iv_indic)
820	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
821	& prodatqc%var(1)%vdep(ista:iend), &
822	& zinmt(iin,ijn,:,iobs), &
823	& zobs2k, interp_corner(iin,ijn,:), &
824	& zgdept(iin,ijn,:,iobs), &
825	& zmask(iin,ijn,:,iobs))
826
827	ENDDO
828	ENDDO
829
830
831	ELSE
832
833	CALL ctl_stop( ' A nonzero' // &
834	& ' number of profile T BUOY data should' // &
835	& ' only occur at the end of a given day' )
836
837	ENDIF
838
839	ELSE
840
841	! Point data
842
843	! vertically interpolate all 4 corners
844	ista = prodatqc%npvsta(jobs,1)
845	iend = prodatqc%npvend(jobs,1)
846	inum_obs = iend - ista + 1
847	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
848	DO iin=1,2
849	DO ijn=1,2
850
851
852	IF ( k1dint == 1 ) THEN
853	CALL obs_int_z1d_spl( kpk, &
854	& zintt(iin,ijn,:,iobs),&
855	& zobs2k, zgdept(iin,ijn,:,iobs), &
856	& zmask(iin,ijn,:,iobs))
857
858	ENDIF
859
860	CALL obs_level_search(kpk, &
861	& zgdept(iin,ijn,:,iobs),&
862	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
863	& iv_indic)
864	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
865	& prodatqc%var(1)%vdep(ista:iend), &
866	& zintt(iin,ijn,:,iobs), &
867	& zobs2k,interp_corner(iin,ijn,:), &
868	& zgdept(iin,ijn,:,iobs), &
869	& zmask(iin,ijn,:,iobs) )
870
871	ENDDO
872	ENDDO
873
874	ENDIF
875
876	!-------------------------------------------------------------
877	! Compute the horizontal interpolation for every profile level
878	!-------------------------------------------------------------
879
880	DO ikn=1,inum_obs
881	iend=ista+ikn-1
882
883	l_zweig(:,:,1) = 0._wp
884
885	! This code forces the horizontal weights to be
886	! zero IF the observation is below the bottom of the
887	! corners of the interpolation nodes, Or if it is in
888	! the mask. This is important for observations are near
889	! steep bathymetry
890	DO iin=1,2
891	DO ijn=1,2
892
893	depth_loop1: DO ik=kpk,2,-1
894	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
895
896	l_zweig(iin,ijn,1) = &
897	& zweig(iin,ijn,1) * &
898	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
899	& - prodatqc%var(1)%vdep(iend)),0._wp)
900
901	EXIT depth_loop1
902	ENDIF
903	ENDDO depth_loop1
904
905	ENDDO
906	ENDDO
907
908	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
909	& prodatqc%var(1)%vmod(iend:iend) )
910
911	ENDDO
912
913
914	DEALLOCATE(interp_corner,iv_indic)
915
916	ENDIF
917
918
919	! Salinity
920
921	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
922
923	zobsk(:) = obfillflt
924
925	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
926
927	IF ( idayend == 0 ) THEN
928
929	! Daily averaged moored buoy (MRB) data
930
931	! vertically interpolate all 4 corners
932	ista = prodatqc%npvsta(iobs,2)
933	iend = prodatqc%npvend(iobs,2)
934	inum_obs = iend - ista + 1
935	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
936
937	DO iin=1,2
938	DO ijn=1,2
939
940
941
942	IF ( k1dint == 1 ) THEN
943	CALL obs_int_z1d_spl( kpk, &
944	& zinms(iin,ijn,:,iobs), &
945	& zobs2k, zgdept(iin,ijn,:,iobs), &
946	& zmask(iin,ijn,:,iobs))
947	ENDIF
948
949	CALL obs_level_search(kpk, &
950	& zgdept(iin,ijn,:,iobs), &
951	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
952	& iv_indic)
953	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
954	& prodatqc%var(2)%vdep(ista:iend), &
955	& zinms(iin,ijn,:,iobs), &
956	& zobs2k, interp_corner(iin,ijn,:), &
957	& zgdept(iin,ijn,:,iobs), &
958	& zmask(iin,ijn,:,iobs))
959
960	ENDDO
961	ENDDO
962
963
964	ELSE
965
966	CALL ctl_stop( ' A nonzero' // &
967	& ' number of profile T BUOY data should' // &
968	& ' only occur at the end of a given day' )
969
970	ENDIF
971
972	ELSE
973
974	! Point data
975
976	! vertically interpolate all 4 corners
977	ista = prodatqc%npvsta(jobs,2)
978	iend = prodatqc%npvend(jobs,2)
979	inum_obs = iend - ista + 1
980	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
981
982	DO iin=1,2
983	DO ijn=1,2
984
985
986	IF ( k1dint == 1 ) THEN
987	CALL obs_int_z1d_spl( kpk, &
988	& zints(iin,ijn,:,iobs),&
989	& zobs2k, zgdept(iin,ijn,:,iobs), &
990	& zmask(iin,ijn,:,iobs))
991
992	ENDIF
993
994	CALL obs_level_search(kpk, &
995	& zgdept(iin,ijn,:,iobs),&
996	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
997	& iv_indic)
998	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
999	& prodatqc%var(2)%vdep(ista:iend), &
1000	& zints(iin,ijn,:,iobs), &
1001	& zobs2k,interp_corner(iin,ijn,:), &
1002	& zgdept(iin,ijn,:,iobs), &
1003	& zmask(iin,ijn,:,iobs) )
1004
1005	ENDDO
1006	ENDDO
1007
1008	ENDIF
1009
1010	!-------------------------------------------------------------
1011	! Compute the horizontal interpolation for every profile level
1012	!-------------------------------------------------------------
1013
1014	DO ikn=1,inum_obs
1015	iend=ista+ikn-1
1016
1017	l_zweig(:,:,1) = 0._wp
1018
1019	! This code forces the horizontal weights to be
1020	! zero IF the observation is below the bottom of the
1021	! corners of the interpolation nodes, Or if it is in
1022	! the mask. This is important for observations are near
1023	! steep bathymetry
1024	DO iin=1,2
1025	DO ijn=1,2
1026
1027	depth_loop2: DO ik=kpk,2,-1
1028	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
1029
1030	l_zweig(iin,ijn,1) = &
1031	& zweig(iin,ijn,1) * &
1032	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
1033	& - prodatqc%var(2)%vdep(iend)),0._wp)
1034
1035	EXIT depth_loop2
1036	ENDIF
1037	ENDDO depth_loop2
1038
1039	ENDDO
1040	ENDDO
1041
1042	CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &
1043	& prodatqc%var(2)%vmod(iend:iend) )
1044
1045	ENDDO
1046
1047
1048	DEALLOCATE(interp_corner,iv_indic)
1049
1050	ENDIF
1051
1052	ENDIF
1053
1054	END DO
1055
1056	! Deallocate the data for interpolation
1057	DEALLOCATE( &
1058	& igrdi, &
1059	& igrdj, &
1060	& zglam, &
1061	& zgphi, &
1062	& zmask, &
1063	& zintt, &
1064	& zints, &
1065	& zgdept,&
1066	& zgdepw &
1067	& )
1068	! At the end of the day also get interpolated means
1069	IF ( idayend == 0 ) THEN
1070	DEALLOCATE( &
1071	& zinmt, &
1072	& zinms &
1073	& )
1074	ENDIF
1075
1076	prodatqc%nprofup = prodatqc%nprofup + ipro
1077
1078	END SUBROUTINE obs_pro_sco_opt
1079
1080	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
1081	& kit000, kdaystp, psurf, psurfmask, &
1082	& k2dint, ldnightav )
1083
1084	!!-----------------------------------------------------------------------
1085	!!
1086	!! * ROUTINE obs_surf_opt *
1087	!!
1088	!! ** Purpose : Compute the model counterpart of surface
1089	!! data by interpolating from the model grid to the
1090	!! observation point.
1091	!!
1092	!! ** Method : Linearly interpolate to each observation point using
1093	!! the model values at the corners of the surrounding grid box.
1094	!!
1095	!! The new model value is first computed at the obs (lon, lat) point.
1096	!!
1097	!! Several horizontal interpolation schemes are available:
1098	!! - distance-weighted (great circle) (k2dint = 0)
1099	!! - distance-weighted (small angle) (k2dint = 1)
1100	!! - bilinear (geographical grid) (k2dint = 2)
1101	!! - bilinear (quadrilateral grid) (k2dint = 3)
1102	!! - polynomial (quadrilateral grid) (k2dint = 4)
1103	!!
1104	!!
1105	!! ** Action :
1106	!!
1107	!! History :
1108	!! ! 07-03 (A. Weaver)
1109	!! ! 15-02 (M. Martin) Combined routine for surface types
1110	!!-----------------------------------------------------------------------
1111
1112	!! * Modules used
1113	USE obs_surf_def ! Definition of storage space for surface observations
1114
1115	IMPLICIT NONE
1116
1117	!! * Arguments
1118	TYPE(obs_surf), INTENT(INOUT) :: &
1119	& surfdataqc ! Subset of surface data passing QC
1120	INTEGER, INTENT(IN) :: kt ! Time step
1121	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
1122	INTEGER, INTENT(IN) :: kpj
1123	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
1124	! (kit000-1 = restart time)
1125	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
1126	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
1127	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1128	& psurf, & ! Model surface field
1129	& psurfmask ! Land-sea mask
1130	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
1131
1132	!! * Local declarations
1133	INTEGER :: ji
1134	INTEGER :: jj
1135	INTEGER :: jobs
1136	INTEGER :: inrc
1137	INTEGER :: isurf
1138	INTEGER :: iobs
1139	INTEGER :: idayend
1140	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1141	& igrdi, &
1142	& igrdj
1143	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1144	& icount_night, &
1145	& imask_night
1146	REAL(wp) :: zlam
1147	REAL(wp) :: zphi
1148	REAL(wp), DIMENSION(1) :: zext, zobsmask
1149	REAL(wp) :: zdaystp
1150	REAL(wp), DIMENSION(2,2,1) :: &
1151	& zweig
1152	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1153	& zmask, &
1154	& zsurf, &
1155	& zsurfm, &
1156	& zglam, &
1157	& zgphi
1158	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1159	& zintmp, &
1160	& zouttmp, &
1161	& zmeanday ! to compute model sst in region of 24h daylight (pole)
1162
1163	!------------------------------------------------------------------------
1164	! Local initialization
1165	!------------------------------------------------------------------------
1166	! Record and data counters
1167	inrc = kt - kit000 + 2
1168	isurf = surfdataqc%nsstp(inrc)
1169
1170	IF ( ldnightav ) THEN
1171
1172	! Initialize array for night mean
1173	IF ( kt == 0 ) THEN
1174	ALLOCATE ( icount_night(kpi,kpj) )
1175	ALLOCATE ( imask_night(kpi,kpj) )
1176	ALLOCATE ( zintmp(kpi,kpj) )
1177	ALLOCATE ( zouttmp(kpi,kpj) )
1178	ALLOCATE ( zmeanday(kpi,kpj) )
1179	nday_qsr = -1 ! initialisation flag for nbc_dcy
1180	ENDIF
1181
1182	! Night-time means are calculated for night-time values over timesteps:
1183	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
1184	idayend = MOD( kt - kit000 + 1, kdaystp )
1185
1186	! Initialize night-time mean for first timestep of the day
1187	IF ( idayend == 1 .OR. kt == 0 ) THEN
1188	DO jj = 1, jpj
1189	DO ji = 1, jpi
1190	surfdataqc%vdmean(ji,jj) = 0.0
1191	zmeanday(ji,jj) = 0.0
1192	icount_night(ji,jj) = 0
1193	END DO
1194	END DO
1195	ENDIF
1196
1197	zintmp(:,:) = 0.0
1198	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1199	imask_night(:,:) = INT( zouttmp(:,:) )
1200
1201	DO jj = 1, jpj
1202	DO ji = 1, jpi
1203	! Increment the temperature field for computing night mean and counter
1204	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
1205	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
1206	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
1207	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
1208	END DO
1209	END DO
1210
1211	! Compute the night-time mean at the end of the day
1212	zdaystp = 1.0 / REAL( kdaystp )
1213	IF ( idayend == 0 ) THEN
1214	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
1215	DO jj = 1, jpj
1216	DO ji = 1, jpi
1217	! Test if "no night" point
1218	IF ( icount_night(ji,jj) > 0 ) THEN
1219	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
1220	& / REAL( icount_night(ji,jj) )
1221	ELSE
1222	!At locations where there is no night (e.g. poles),
1223	! calculate daily mean instead of night-time mean.
1224	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1225	ENDIF
1226	END DO
1227	END DO
1228	ENDIF
1229
1230	ENDIF
1231
1232	! Get the data for interpolation
1233
1234	ALLOCATE( &
1235	& igrdi(2,2,isurf), &
1236	& igrdj(2,2,isurf), &
1237	& zglam(2,2,isurf), &
1238	& zgphi(2,2,isurf), &
1239	& zmask(2,2,isurf), &
1240	& zsurf(2,2,isurf) &
1241	& )
1242
1243	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
1244	iobs = jobs - surfdataqc%nsurfup
1245	igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1
1246	igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1
1247	igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1
1248	igrdj(1,2,iobs) = surfdataqc%mj(jobs)
1249	igrdi(2,1,iobs) = surfdataqc%mi(jobs)
1250	igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1
1251	igrdi(2,2,iobs) = surfdataqc%mi(jobs)
1252	igrdj(2,2,iobs) = surfdataqc%mj(jobs)
1253	END DO
1254
1255	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
1256	& igrdi, igrdj, glamt, zglam )
1257	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
1258	& igrdi, igrdj, gphit, zgphi )
1259	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
1260	& igrdi, igrdj, psurfmask, zmask )
1261	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
1262	& igrdi, igrdj, psurf, zsurf )
1263
1264	! At the end of the day get interpolated means
1265	IF ( idayend == 0 .AND. ldnightav ) THEN
1266
1267	ALLOCATE( &
1268	& zsurfm(2,2,isurf) &
1269	& )
1270
1271	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, &
1272	& surfdataqc%vdmean(:,:), zsurfm )
1273
1274	ENDIF
1275
1276	! Loop over observations
1277	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
1278
1279	iobs = jobs - surfdataqc%nsurfup
1280
1281	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
1282
1283	IF(lwp) THEN
1284	WRITE(numout,*)
1285	WRITE(numout,*) ' E R R O R : Observation', &
1286	& ' time step is not consistent with the', &
1287	& ' model time step'
1288	WRITE(numout,*) ' ========='
1289	WRITE(numout,*)
1290	WRITE(numout,*) ' Record = ', jobs, &
1291	& ' kt = ', kt, &
1292	& ' mstp = ', surfdataqc%mstp(jobs), &
1293	& ' ntyp = ', surfdataqc%ntyp(jobs)
1294	ENDIF
1295	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
1296
1297	ENDIF
1298
1299	zlam = surfdataqc%rlam(jobs)
1300	zphi = surfdataqc%rphi(jobs)
1301
1302	! Get weights to interpolate the model value to the observation point
1303	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
1304	& zglam(:,:,iobs), zgphi(:,:,iobs), &
1305	& zmask(:,:,iobs), zweig, zobsmask )
1306
1307	! Interpolate the model field to the observation point
1308	IF ( ldnightav .AND. idayend == 0 ) THEN
1309	! Night-time averaged data
1310	CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext )
1311	ELSE
1312	CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext )
1313	ENDIF
1314
1315	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
1316	! ... Remove the MDT from the SSH at the observation point to get the SLA
1317	surfdataqc%rext(jobs,1) = zext(1)
1318	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
1319	ELSE
1320	surfdataqc%rmod(jobs,1) = zext(1)
1321	ENDIF
1322
1323	END DO
1324
1325	! Deallocate the data for interpolation
1326	DEALLOCATE( &
1327	& igrdi, &
1328	& igrdj, &
1329	& zglam, &
1330	& zgphi, &
1331	& zmask, &
1332	& zsurf &
1333	& )
1334
1335	! At the end of the day also deallocate night-time mean array
1336	IF ( idayend == 0 .AND. ldnightav ) THEN
1337	DEALLOCATE( &
1338	& zsurfm &
1339	& )
1340	ENDIF
1341
1342	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
1343
1344	END SUBROUTINE obs_surf_opt
1345
1346	END MODULE obs_oper

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

Download in other formats: