Context Navigation

obs_oper.F90 @ 7773

Last change on this file since 7773 was 7773, checked in by mattmartin, 7 years ago

Committing updates after doing the following:

merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
- generalised vertical coordinate for profile obs. This was done so that is now the default option.
- sst bias correction implemented with the new simplified obs code.
- included the biogeochemical obs types int he new simplified obs code.
- included the changes to exclude obs in the boundary for limited area models
- included other changes for the efficiency of the obs operator to remove global arrays.

File size: 34.7 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	!!----------------------------------------------------------------------
12
13	!! * Modules used
14	USE par_kind, ONLY : & ! Precision variables
15	& wp
16	USE in_out_manager ! I/O manager
17	USE obs_inter_sup ! Interpolation support
18	USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt
19	& obs_int_h2d, &
20	& obs_int_h2d_init
21	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt
22	& obs_int_z1d, &
23	& obs_int_z1d_spl
24	USE obs_const, ONLY : & ! Obs fill value
25	& obfillflt
26	USE dom_oce, ONLY : & ! Model lats/lons
27	& glamt, &
28	& gphit
29	USE lib_mpp, ONLY : & ! Warning and stopping routines
30	& ctl_warn, ctl_stop
31	USE sbcdcy, ONLY : & ! For calculation of where it is night-time
32	& sbc_dcy, nday_qsr
33
34	IMPLICIT NONE
35
36	!! * Routine accessibility
37	PRIVATE
38
39	PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs
40	& obs_surf_opt ! Compute the model counterpart of surface obs
41
42	INTEGER, PARAMETER, PUBLIC :: &
43	& imaxavtypes = 20 ! Max number of daily avgd obs types
44
45	!!----------------------------------------------------------------------
46	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
47	!! $Id$
48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49	!!----------------------------------------------------------------------
50
51	!! * Substitutions
52	# include "domzgr_substitute.h90"
53	CONTAINS
54
55	SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
56	& kit000, kdaystp, &
57	& pvar1, pvar2, pgdept, pgdepw,
58	& pmask1, pmask2, &
59	& plam1, plam2, pphi1, pphi2, &
60	& k1dint, k2dint, kdailyavtypes )
61
62	!!-----------------------------------------------------------------------
63	!!
64	!! * ROUTINE obs_pro_opt *
65	!!
66	!! ** Purpose : Compute the model counterpart of profiles
67	!! data by interpolating from the model grid to the
68	!! observation point.
69	!!
70	!! ** Method : Linearly interpolate to each observation point using
71	!! the model values at the corners of the surrounding grid box.
72	!!
73	!! First, a vertical profile of horizontally interpolated model
74	!! now values is computed at the obs (lon, lat) point.
75	!! Several horizontal interpolation schemes are available:
76	!! - distance-weighted (great circle) (k2dint = 0)
77	!! - distance-weighted (small angle) (k2dint = 1)
78	!! - bilinear (geographical grid) (k2dint = 2)
79	!! - bilinear (quadrilateral grid) (k2dint = 3)
80	!! - polynomial (quadrilateral grid) (k2dint = 4)
81	!!
82	!! Next, the vertical profile is interpolated to the
83	!! data depth points. Two vertical interpolation schemes are
84	!! available:
85	!! - linear (k1dint = 0)
86	!! - Cubic spline (k1dint = 1)
87	!!
88	!! For the cubic spline the 2nd derivative of the interpolating
89	!! polynomial is computed before entering the vertical interpolation
90	!! routine.
91	!!
92	!! If the logical is switched on, the model equivalent is
93	!! a daily mean model temperature field. So, we first compute
94	!! the mean, then interpolate only at the end of the day.
95	!!
96	!! Note: in situ temperature observations must be converted
97	!! to potential temperature (the model variable) prior to
98	!! assimilation.
99	!!
100	!! ** Action :
101	!!
102	!! History :
103	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
104	!! ! 06-03 (G. Smith) NEMOVAR migration
105	!! ! 06-10 (A. Weaver) Cleanup
106	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
107	!! ! 07-03 (K. Mogensen) General handling of profiles
108	!! ! 15-02 (M. Martin) Combined routine for all profile types
109	!! ! 17-02 (M. Martin) Include generalised vertical coordinate changes
110	!!-----------------------------------------------------------------------
111
112	!! * Modules used
113	USE obs_profiles_def ! Definition of storage space for profile obs.
114
115	IMPLICIT NONE
116
117	!! * Arguments
118	TYPE(obs_prof), INTENT(INOUT) :: &
119	& prodatqc ! Subset of profile data passing QC
120	INTEGER, INTENT(IN) :: kt ! Time step
121	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
122	INTEGER, INTENT(IN) :: kpj
123	INTEGER, INTENT(IN) :: kpk
124	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
125	! (kit000-1 = restart time)
126	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
127	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
128	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
129	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
130	& pvar1, & ! Model field 1
131	& pvar2, & ! Model field 2
132	& pmask1, & ! Land-sea mask 1
133	& pmask2 ! Land-sea mask 2
134	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
135	& plam1, & ! Model longitudes for variable 1
136	& plam2, & ! Model longitudes for variable 2
137	& pphi1, & ! Model latitudes for variable 1
138	& pphi2 ! Model latitudes for variable 2
139	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
140	& pgdept, & ! Model array of depth T levels
141	& pgdepw ! Model array of depth W levels
142	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
143	& kdailyavtypes ! Types for daily averages
144
145	!! * Local declarations
146	INTEGER :: ji
147	INTEGER :: jj
148	INTEGER :: jk
149	INTEGER :: jobs
150	INTEGER :: inrc
151	INTEGER :: ipro
152	INTEGER :: idayend
153	INTEGER :: ista
154	INTEGER :: iend
155	INTEGER :: iobs
156	INTEGER, DIMENSION(imaxavtypes) :: &
157	& idailyavtypes
158	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
159	& igrdi1, &
160	& igrdi2, &
161	& igrdj1, &
162	& igrdj2
163	REAL(KIND=wp) :: zlam
164	REAL(KIND=wp) :: zphi
165	REAL(KIND=wp) :: zdaystp
166	REAL(KIND=wp), DIMENSION(kpk) :: &
167	& zobsmask1, &
168	& zobsmask2, &
169	& zobsk, &
170	& zobs2k
171	REAL(KIND=wp), DIMENSION(2,2,1) :: &
172	& zweig1, &
173	& zweig2, &
174	& zweig
175	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
176	& zmask1, &
177	& zmask2, &
178	& zint1, &
179	& zint2, &
180	& zinm1, &
181	& zinm2, &
182	& zgdept, &
183	& zgdepw
184	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
185	& zglam1, &
186	& zglam2, &
187	& zgphi1, &
188	& zgphi2
189	REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2
190	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
191
192	LOGICAL :: ld_dailyav
193
194	!------------------------------------------------------------------------
195	! Local initialization
196	!------------------------------------------------------------------------
197	! Record and data counters
198	inrc = kt - kit000 + 2
199	ipro = prodatqc%npstp(inrc)
200
201	! Daily average types
202	ld_dailyav = .FALSE.
203	IF ( PRESENT(kdailyavtypes) ) THEN
204	idailyavtypes(:) = kdailyavtypes(:)
205	IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
206	ELSE
207	idailyavtypes(:) = -1
208	ENDIF
209
210	! Daily means are calculated for values over timesteps:
211	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
212	idayend = MOD( kt - kit000 + 1, kdaystp )
213
214	IF ( ld_dailyav ) THEN
215
216	! Initialize daily mean for first timestep of the day
217	IF ( idayend == 1 .OR. kt == 0 ) THEN
218	DO jk = 1, jpk
219	DO jj = 1, jpj
220	DO ji = 1, jpi
221	prodatqc%vdmean(ji,jj,jk,1) = 0.0
222	prodatqc%vdmean(ji,jj,jk,2) = 0.0
223	END DO
224	END DO
225	END DO
226	ENDIF
227
228	DO jk = 1, jpk
229	DO jj = 1, jpj
230	DO ji = 1, jpi
231	! Increment field 1 for computing daily mean
232	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
233	& + pvar1(ji,jj,jk)
234	! Increment field 2 for computing daily mean
235	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
236	& + pvar2(ji,jj,jk)
237	END DO
238	END DO
239	END DO
240
241	! Compute the daily mean at the end of day
242	zdaystp = 1.0 / REAL( kdaystp )
243	IF ( idayend == 0 ) THEN
244	IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
245	CALL FLUSH(numout)
246	DO jk = 1, jpk
247	DO jj = 1, jpj
248	DO ji = 1, jpi
249	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
250	& * zdaystp
251	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
252	& * zdaystp
253	END DO
254	END DO
255	END DO
256	ENDIF
257
258	ENDIF
259
260	! Get the data for interpolation
261	ALLOCATE( &
262	& igrdi1(2,2,ipro), &
263	& igrdi2(2,2,ipro), &
264	& igrdj1(2,2,ipro), &
265	& igrdj2(2,2,ipro), &
266	& zglam1(2,2,ipro), &
267	& zglam2(2,2,ipro), &
268	& zgphi1(2,2,ipro), &
269	& zgphi2(2,2,ipro), &
270	& zmask1(2,2,kpk,ipro), &
271	& zmask2(2,2,kpk,ipro), &
272	& zint1(2,2,kpk,ipro), &
273	& zint2(2,2,kpk,ipro), &
274	& zgdept(2,2,kpk,ipro), &
275	& zgdepw(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	! Initialise depth arrays
299	zgdept(:,:,:,:) = 0.0
300	zgdepw(:,:,:,:) = 0.0
301
302	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
303	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
304	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
305	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 )
306
307	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
308	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
309	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
310	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 )
311
312	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )
313	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )
314
315	! At the end of the day also get interpolated means
316	IF ( ld_dailyav .AND. idayend == 0 ) THEN
317
318	ALLOCATE( &
319	& zinm1(2,2,kpk,ipro), &
320	& zinm2(2,2,kpk,ipro) &
321	& )
322
323	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
324	& prodatqc%vdmean(:,:,:,1), zinm1 )
325	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
326	& prodatqc%vdmean(:,:,:,2), zinm2 )
327
328	ENDIF
329
330	! Return if no observations to process
331	! Has to be done after comm commands to ensure processors
332	! stay in sync
333	IF ( ipro == 0 ) RETURN
334
335	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
336
337	iobs = jobs - prodatqc%nprofup
338
339	IF ( kt /= prodatqc%mstp(jobs) ) THEN
340
341	IF(lwp) THEN
342	WRITE(numout,*)
343	WRITE(numout,*) ' E R R O R : Observation', &
344	& ' time step is not consistent with the', &
345	& ' model time step'
346	WRITE(numout,*) ' ========='
347	WRITE(numout,*)
348	WRITE(numout,*) ' Record = ', jobs, &
349	& ' kt = ', kt, &
350	& ' mstp = ', prodatqc%mstp(jobs), &
351	& ' ntyp = ', prodatqc%ntyp(jobs)
352	ENDIF
353	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
354	ENDIF
355
356	zlam = prodatqc%rlam(jobs)
357	zphi = prodatqc%rphi(jobs)
358
359	! Horizontal weights
360	! Masked values are calculated later.
361	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
362
363	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
364	& zglam1(:,:,iobs), zgphi1(:,:,iobs), &
365	& zmask1(:,:,1,iobs), zweig1, zmsk_1 )
366
367	ENDIF
368
369	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
370
371	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
372	& zglam2(:,:,iobs), zgphi2(:,:,iobs), &
373	& zmask2(:,:,1,iobs), zweig2, zmsk_2 )
374
375	ENDIF
376
377	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
378
379	zobsk(:) = obfillflt
380
381	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
382
383	IF ( idayend == 0 ) THEN
384	! Daily averaged data
385
386	! vertically interpolate all 4 corners
387	ista = prodatqc%npvsta(jobs,1)
388	iend = prodatqc%npvend(jobs,1)
389	inum_obs = iend - ista + 1
390	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
391
392	DO iin=1,2
393	DO ijn=1,2
394
395	IF ( k1dint == 1 ) THEN
396	CALL obs_int_z1d_spl( kpk, &
397	& zinm1(iin,ijn,:,iobs), &
398	& zobs2k, zgdept(iin,ijn,:,iobs), &
399	& zmask1(iin,ijn,:,iobs))
400	ENDIF
401
402	CALL obs_level_search(kpk, &
403	& zgdept(iin,ijn,:,iobs), &
404	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
405	& iv_indic)
406
407	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
408	& prodatqc%var(1)%vdep(ista:iend), &
409	& zinm1(iin,ijn,:,iobs), &
410	& zobs2k, interp_corner(iin,ijn,:), &
411	& zgdept(iin,ijn,:,iobs), &
412	& zmask1(iin,ijn,:,iobs))
413
414	ENDDO
415	ENDDO
416
417	ENDIF !idayend
418
419	ELSE
420
421	! Point data
422
423	! vertically interpolate all 4 corners
424	ista = prodatqc%npvsta(jobs,1)
425	iend = prodatqc%npvend(jobs,1)
426	inum_obs = iend - ista + 1
427	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
428	DO iin=1,2
429	DO ijn=1,2
430
431	IF ( k1dint == 1 ) THEN
432	CALL obs_int_z1d_spl( kpk, &
433	& zint1(iin,ijn,:,iobs),&
434	& zobs2k, zgdept(iin,ijn,:,iobs), &
435	& zmask1(iin,ijn,:,iobs))
436
437	ENDIF
438
439	CALL obs_level_search(kpk, &
440	& zgdept(iin,ijn,:,iobs),&
441	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
442	& iv_indic)
443
444	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
445	& prodatqc%var(1)%vdep(ista:iend), &
446	& zint1(iin,ijn,:,iobs), &
447	& zobs2k,interp_corner(iin,ijn,:), &
448	& zgdept(iin,ijn,:,iobs), &
449	& zmask1(iin,ijn,:,iobs) )
450
451	ENDDO
452	ENDDO
453
454	ENDIF
455
456	!-------------------------------------------------------------
457	! Compute the horizontal interpolation for every profile level
458	!-------------------------------------------------------------
459
460	DO ikn=1,inum_obs
461	iend=ista+ikn-1
462
463	zweig(:,:,1) = 0._wp
464
465	! This code forces the horizontal weights to be
466	! zero IF the observation is below the bottom of the
467	! corners of the interpolation nodes, Or if it is in
468	! the mask. This is important for observations near
469	! steep bathymetry
470	DO iin=1,2
471	DO ijn=1,2
472
473	depth_loop1: DO ik=kpk,2,-1
474	IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN
475
476	zweig(iin,ijn,1) = &
477	& zweig1(iin,ijn,1) * &
478	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
479	& - prodatqc%var(1)%vdep(iend)),0._wp)
480
481	EXIT depth_loop1
482
483	ENDIF
484
485	ENDDO depth_loop1
486
487	ENDDO
488	ENDDO
489
490	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
491	& prodatqc%var(1)%vmod(iend:iend) )
492
493	! Set QC flag for any observations found below the bottom
494	! needed as the check here is more strict than that in obs_prep
495	IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
496
497	ENDDO
498
499	DEALLOCATE(interp_corner,iv_indic)
500
501	ENDIF
502
503	! For the second variable
504	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
505
506	zobsk(:) = obfillflt
507
508	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
509
510	IF ( idayend == 0 ) THEN
511	! Daily averaged data
512
513	! vertically interpolate all 4 corners
514	ista = prodatqc%npvsta(jobs,2)
515	iend = prodatqc%npvend(jobs,2)
516	inum_obs = iend - ista + 1
517	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
518
519	DO iin=1,2
520	DO ijn=1,2
521
522	IF ( k1dint == 1 ) THEN
523	CALL obs_int_z1d_spl( kpk, &
524	& zinm2(iin,ijn,:,iobs), &
525	& zobs2k, zgdept(iin,ijn,:,iobs), &
526	& zmask2(iin,ijn,:,iobs))
527	ENDIF
528
529	CALL obs_level_search(kpk, &
530	& zgdept(iin,ijn,:,iobs), &
531	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
532	& iv_indic)
533
534	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
535	& prodatqc%var(2)%vdep(ista:iend), &
536	& zinm2(iin,ijn,:,iobs), &
537	& zobs2k, interp_corner(iin,ijn,:), &
538	& zgdept(iin,ijn,:,iobs), &
539	& zmask2(iin,ijn,:,iobs))
540
541	ENDDO
542	ENDDO
543
544	ENDIF !idayend
545
546	ELSE
547
548	! Point data
549
550	! vertically interpolate all 4 corners
551	ista = prodatqc%npvsta(jobs,2)
552	iend = prodatqc%npvend(jobs,2)
553	inum_obs = iend - ista + 1
554	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
555	DO iin=1,2
556	DO ijn=1,2
557
558	IF ( k1dint == 1 ) THEN
559	CALL obs_int_z1d_spl( kpk, &
560	& zint2(iin,ijn,:,iobs),&
561	& zobs2k, zgdept(iin,ijn,:,iobs), &
562	& zmask2(iin,ijn,:,iobs))
563
564	ENDIF
565
566	CALL obs_level_search(kpk, &
567	& zgdept(iin,ijn,:,iobs),&
568	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
569	& iv_indic)
570
571	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
572	& prodatqc%var(2)%vdep(ista:iend), &
573	& zint2(iin,ijn,:,iobs), &
574	& zobs2k,interp_corner(iin,ijn,:), &
575	& zgdept(iin,ijn,:,iobs), &
576	& zmask2(iin,ijn,:,iobs) )
577
578	ENDDO
579	ENDDO
580
581	ENDIF
582
583	!-------------------------------------------------------------
584	! Compute the horizontal interpolation for every profile level
585	!-------------------------------------------------------------
586
587	DO ikn=1,inum_obs
588	iend=ista+ikn-1
589
590	zweig(:,:,1) = 0._wp
591
592	! This code forces the horizontal weights to be
593	! zero IF the observation is below the bottom of the
594	! corners of the interpolation nodes, Or if it is in
595	! the mask. This is important for observations near
596	! steep bathymetry
597	DO iin=1,2
598	DO ijn=1,2
599
600	depth_loop2: DO ik=kpk,2,-1
601	IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN
602
603	zweig(iin,ijn,1) = &
604	& zweig2(iin,ijn,1) * &
605	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
606	& - prodatqc%var(2)%vdep(iend)),0._wp)
607
608	EXIT depth_loop2
609
610	ENDIF
611
612	ENDDO depth_loop2
613
614	ENDDO
615	ENDDO
616
617	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
618	& prodatqc%var(2)%vmod(iend:iend) )
619
620	! Set QC flag for any observations found below the bottom
621	! needed as the check here is more strict than that in obs_prep
622	IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
623
624	ENDDO
625
626	DEALLOCATE(interp_corner,iv_indic)
627
628	ENDIF
629
630	! Deallocate the data for interpolation
631	DEALLOCATE( &
632	& igrdi1, &
633	& igrdi2, &
634	& igrdj1, &
635	& igrdj2, &
636	& zglam1, &
637	& zglam2, &
638	& zgphi1, &
639	& zgphi2, &
640	& zmask1, &
641	& zmask2, &
642	& zint1, &
643	& zint2, &
644	& zgdept, &
645	& zgdepw &
646	& )
647
648	! At the end of the day also get interpolated means
649	IF ( ld_dailyav .AND. idayend == 0 ) THEN
650	DEALLOCATE( &
651	& zinm1, &
652	& zinm2 &
653	& )
654	ENDIF
655
656	prodatqc%nprofup = prodatqc%nprofup + ipro
657
658	END SUBROUTINE obs_prof_opt
659
660	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
661	& kit000, kdaystp, psurf, psurfmask, &
662	& k2dint, ldnightav )
663
664	!!-----------------------------------------------------------------------
665	!!
666	!! * ROUTINE obs_surf_opt *
667	!!
668	!! ** Purpose : Compute the model counterpart of surface
669	!! data by interpolating from the model grid to the
670	!! observation point.
671	!!
672	!! ** Method : Linearly interpolate to each observation point using
673	!! the model values at the corners of the surrounding grid box.
674	!!
675	!! The new model value is first computed at the obs (lon, lat) point.
676	!!
677	!! Several horizontal interpolation schemes are available:
678	!! - distance-weighted (great circle) (k2dint = 0)
679	!! - distance-weighted (small angle) (k2dint = 1)
680	!! - bilinear (geographical grid) (k2dint = 2)
681	!! - bilinear (quadrilateral grid) (k2dint = 3)
682	!! - polynomial (quadrilateral grid) (k2dint = 4)
683	!!
684	!!
685	!! ** Action :
686	!!
687	!! History :
688	!! ! 07-03 (A. Weaver)
689	!! ! 15-02 (M. Martin) Combined routine for surface types
690	!!-----------------------------------------------------------------------
691
692	!! * Modules used
693	USE obs_surf_def ! Definition of storage space for surface observations
694
695	IMPLICIT NONE
696
697	!! * Arguments
698	TYPE(obs_surf), INTENT(INOUT) :: &
699	& surfdataqc ! Subset of surface data passing QC
700	INTEGER, INTENT(IN) :: kt ! Time step
701	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
702	INTEGER, INTENT(IN) :: kpj
703	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
704	! (kit000-1 = restart time)
705	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
706	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
707	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
708	& psurf, & ! Model surface field
709	& psurfmask ! Land-sea mask
710	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
711
712	!! * Local declarations
713	INTEGER :: ji
714	INTEGER :: jj
715	INTEGER :: jobs
716	INTEGER :: inrc
717	INTEGER :: isurf
718	INTEGER :: iobs
719	INTEGER :: idayend
720	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
721	& igrdi, &
722	& igrdj
723	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
724	& icount_night, &
725	& imask_night
726	REAL(wp) :: zlam
727	REAL(wp) :: zphi
728	REAL(wp), DIMENSION(1) :: zext, zobsmask
729	REAL(wp) :: zdaystp
730	REAL(wp), DIMENSION(2,2,1) :: &
731	& zweig
732	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
733	& zmask, &
734	& zsurf, &
735	& zsurfm, &
736	& zglam, &
737	& zgphi
738	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
739	& zintmp, &
740	& zouttmp, &
741	& zmeanday ! to compute model sst in region of 24h daylight (pole)
742
743	!------------------------------------------------------------------------
744	! Local initialization
745	!------------------------------------------------------------------------
746	! Record and data counters
747	inrc = kt - kit000 + 2
748	isurf = surfdataqc%nsstp(inrc)
749
750	IF ( ldnightav ) THEN
751
752	! Initialize array for night mean
753	IF ( kt == 0 ) THEN
754	ALLOCATE ( icount_night(kpi,kpj) )
755	ALLOCATE ( imask_night(kpi,kpj) )
756	ALLOCATE ( zintmp(kpi,kpj) )
757	ALLOCATE ( zouttmp(kpi,kpj) )
758	ALLOCATE ( zmeanday(kpi,kpj) )
759	nday_qsr = -1 ! initialisation flag for nbc_dcy
760	ENDIF
761
762	! Night-time means are calculated for night-time values over timesteps:
763	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
764	idayend = MOD( kt - kit000 + 1, kdaystp )
765
766	! Initialize night-time mean for first timestep of the day
767	IF ( idayend == 1 .OR. kt == 0 ) THEN
768	DO jj = 1, jpj
769	DO ji = 1, jpi
770	surfdataqc%vdmean(ji,jj) = 0.0
771	zmeanday(ji,jj) = 0.0
772	icount_night(ji,jj) = 0
773	END DO
774	END DO
775	ENDIF
776
777	zintmp(:,:) = 0.0
778	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
779	imask_night(:,:) = INT( zouttmp(:,:) )
780
781	DO jj = 1, jpj
782	DO ji = 1, jpi
783	! Increment the temperature field for computing night mean and counter
784	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
785	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
786	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
787	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
788	END DO
789	END DO
790
791	! Compute the night-time mean at the end of the day
792	zdaystp = 1.0 / REAL( kdaystp )
793	IF ( idayend == 0 ) THEN
794	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
795	DO jj = 1, jpj
796	DO ji = 1, jpi
797	! Test if "no night" point
798	IF ( icount_night(ji,jj) > 0 ) THEN
799	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
800	& / REAL( icount_night(ji,jj) )
801	ELSE
802	!At locations where there is no night (e.g. poles),
803	! calculate daily mean instead of night-time mean.
804	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
805	ENDIF
806	END DO
807	END DO
808	ENDIF
809
810	ENDIF
811
812	! Get the data for interpolation
813
814	ALLOCATE( &
815	& igrdi(2,2,isurf), &
816	& igrdj(2,2,isurf), &
817	& zglam(2,2,isurf), &
818	& zgphi(2,2,isurf), &
819	& zmask(2,2,isurf), &
820	& zsurf(2,2,isurf) &
821	& )
822
823	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
824	iobs = jobs - surfdataqc%nsurfup
825	igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1
826	igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1
827	igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1
828	igrdj(1,2,iobs) = surfdataqc%mj(jobs)
829	igrdi(2,1,iobs) = surfdataqc%mi(jobs)
830	igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1
831	igrdi(2,2,iobs) = surfdataqc%mi(jobs)
832	igrdj(2,2,iobs) = surfdataqc%mj(jobs)
833	END DO
834
835	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
836	& igrdi, igrdj, glamt, zglam )
837	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
838	& igrdi, igrdj, gphit, zgphi )
839	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
840	& igrdi, igrdj, psurfmask, zmask )
841	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &
842	& igrdi, igrdj, psurf, zsurf )
843
844	! At the end of the day get interpolated means
845	IF ( idayend == 0 .AND. ldnightav ) THEN
846
847	ALLOCATE( &
848	& zsurfm(2,2,isurf) &
849	& )
850
851	CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, &
852	& surfdataqc%vdmean(:,:), zsurfm )
853
854	ENDIF
855
856	! Loop over observations
857	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
858
859	iobs = jobs - surfdataqc%nsurfup
860
861	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
862
863	IF(lwp) THEN
864	WRITE(numout,*)
865	WRITE(numout,*) ' E R R O R : Observation', &
866	& ' time step is not consistent with the', &
867	& ' model time step'
868	WRITE(numout,*) ' ========='
869	WRITE(numout,*)
870	WRITE(numout,*) ' Record = ', jobs, &
871	& ' kt = ', kt, &
872	& ' mstp = ', surfdataqc%mstp(jobs), &
873	& ' ntyp = ', surfdataqc%ntyp(jobs)
874	ENDIF
875	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
876
877	ENDIF
878
879	zlam = surfdataqc%rlam(jobs)
880	zphi = surfdataqc%rphi(jobs)
881
882	! Get weights to interpolate the model value to the observation point
883	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
884	& zglam(:,:,iobs), zgphi(:,:,iobs), &
885	& zmask(:,:,iobs), zweig, zobsmask )
886
887	! Interpolate the model field to the observation point
888	IF ( ldnightav .AND. idayend == 0 ) THEN
889	! Night-time averaged data
890	CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext )
891	ELSE
892	CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext )
893	ENDIF
894
895	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
896	! ... Remove the MDT from the SSH at the observation point to get the SLA
897	surfdataqc%rext(jobs,1) = zext(1)
898	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
899	ELSE
900	surfdataqc%rmod(jobs,1) = zext(1)
901	ENDIF
902
903	END DO
904
905	! Deallocate the data for interpolation
906	DEALLOCATE( &
907	& igrdi, &
908	& igrdj, &
909	& zglam, &
910	& zgphi, &
911	& zmask, &
912	& zsurf &
913	& )
914
915	! At the end of the day also deallocate night-time mean array
916	IF ( idayend == 0 .AND. ldnightav ) THEN
917	DEALLOCATE( &
918	& zsurfm &
919	& )
920	ENDIF
921
922	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
923
924	END SUBROUTINE obs_surf_opt
925
926	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: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 7773

Download in other formats: