Context Navigation

obs_oper.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 8 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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

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

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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 7646

Download in other formats: