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/UKMO/dev_r5518_obs_oper_update_utils366_fabmv1_v2/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_utils366_fabmv1_v2/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 13487

Last change on this file since 13487 was 13393, checked in by mattmartin, 4 years ago
Include surface velocity observation operator. See ticket https://code.metoffice.gov.uk/trac/utils/ticket/376.
File size: 39.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	!!----------------------------------------------------------------------
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_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint
22	& obs_avg_h2d, &
23	& obs_avg_h2d_init, &
24	& obs_max_fpsize
25	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt
26	& obs_int_z1d, &
27	& obs_int_z1d_spl
28	USE obs_const, ONLY : & ! Obs fill value
29	& obfillflt
30	USE dom_oce, ONLY : &
31	& glamt, glamf, &
32	& gphit, gphif
33	USE lib_mpp, ONLY : & ! Warning and stopping routines
34	& ctl_warn, ctl_stop
35	USE sbcdcy, ONLY : & ! For calculation of where it is night-time
36	& sbc_dcy, nday_qsr
37	USE obs_grid, ONLY : &
38	& obs_level_search
39	#if defined key_cice
40	USE ice_constants, ONLY : & ! For conversion from sea ice freeboard to thickness
41	& rhos, rhoi, rhow
42	#endif
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_surf_opt ! Compute the model counterpart of surface obs
51
52	INTEGER, PARAMETER, PUBLIC :: &
53	& imaxavtypes = 20 ! Max number of daily avgd obs types
54
55	!!----------------------------------------------------------------------
56	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
57	!! $Id$
58	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
59	!!----------------------------------------------------------------------
60
61	!! * Substitutions
62	# include "domzgr_substitute.h90"
63	CONTAINS
64
65
66	SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
67	& kit000, kdaystp, kvar, &
68	& pvar, pclim, &
69	& pgdept, pgdepw, pmask, &
70	& plam, pphi, &
71	& k1dint, k2dint, kdailyavtypes )
72
73	!!-----------------------------------------------------------------------
74	!!
75	!! * ROUTINE obs_pro_opt *
76	!!
77	!! ** Purpose : Compute the model counterpart of profiles
78	!! data by interpolating from the model grid to the
79	!! observation point.
80	!!
81	!! ** Method : Linearly interpolate to each observation point using
82	!! the model values at the corners of the surrounding grid box.
83	!!
84	!! First, a vertical profile of horizontally interpolated model
85	!! now values is computed at the obs (lon, lat) point.
86	!! Several horizontal interpolation schemes are available:
87	!! - distance-weighted (great circle) (k2dint = 0)
88	!! - distance-weighted (small angle) (k2dint = 1)
89	!! - bilinear (geographical grid) (k2dint = 2)
90	!! - bilinear (quadrilateral grid) (k2dint = 3)
91	!! - polynomial (quadrilateral grid) (k2dint = 4)
92	!!
93	!! Next, the vertical profile is interpolated to the
94	!! data depth points. Two vertical interpolation schemes are
95	!! available:
96	!! - linear (k1dint = 0)
97	!! - Cubic spline (k1dint = 1)
98	!!
99	!! For the cubic spline the 2nd derivative of the interpolating
100	!! polynomial is computed before entering the vertical interpolation
101	!! routine.
102	!!
103	!! If the logical is switched on, the model equivalent is
104	!! a daily mean model temperature field. So, we first compute
105	!! the mean, then interpolate only at the end of the day.
106	!!
107	!! Note: in situ temperature observations must be converted
108	!! to potential temperature (the model variable) prior to
109	!! assimilation.
110	!!
111	!! ** Action :
112	!!
113	!! History :
114	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
115	!! ! 06-03 (G. Smith) NEMOVAR migration
116	!! ! 06-10 (A. Weaver) Cleanup
117	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
118	!! ! 07-03 (K. Mogensen) General handling of profiles
119	!! ! 15-02 (M. Martin) Combined routine for all profile types
120	!! ! 17-02 (M. Martin) Include generalised vertical coordinate changes
121	!!-----------------------------------------------------------------------
122
123	!! * Modules used
124	USE obs_profiles_def ! Definition of storage space for profile obs.
125
126	IMPLICIT NONE
127
128	!! * Arguments
129	TYPE(obs_prof), INTENT(INOUT) :: &
130	& prodatqc ! Subset of profile data passing QC
131	INTEGER, INTENT(IN) :: kt ! Time step
132	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
133	INTEGER, INTENT(IN) :: kpj
134	INTEGER, INTENT(IN) :: kpk
135	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
136	! (kit000-1 = restart time)
137	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
138	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
139	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
140	INTEGER, INTENT(IN) :: kvar ! Number of variable in prodatqc
141	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
142	& pvar, & ! Model field for variable
143	& pclim, & ! Climatology field for variable
144	& pmask ! Land-sea mask for variable
145	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
146	& plam, & ! Model longitudes for variable
147	& pphi ! Model latitudes for variable
148	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
149	& pgdept, & ! Model array of depth T levels
150	& pgdepw ! Model array of depth W levels
151	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
152	& kdailyavtypes ! Types for daily averages
153
154	!! * Local declarations
155	INTEGER :: ji
156	INTEGER :: jj
157	INTEGER :: jk
158	INTEGER :: jobs
159	INTEGER :: inrc
160	INTEGER :: ipro
161	INTEGER :: idayend
162	INTEGER :: ista
163	INTEGER :: iend
164	INTEGER :: iobs
165	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
166	INTEGER :: inum_obs
167	INTEGER, DIMENSION(imaxavtypes) :: &
168	& idailyavtypes
169	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
170	& igrdi, &
171	& igrdj
172	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
173
174	REAL(KIND=wp) :: zlam
175	REAL(KIND=wp) :: zphi
176	REAL(KIND=wp) :: zdaystp
177	REAL(KIND=wp), DIMENSION(kpk) :: &
178	& zobsk, &
179	& zobs2k, &
180	& zclm2k
181	REAL(KIND=wp), DIMENSION(2,2,1) :: &
182	& zweig1, &
183	& zweig
184	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
185	& zmask, &
186	& zclim, &
187	& zint, &
188	& zinm, &
189	& zgdept, &
190	& zgdepw
191	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
192	& zglam, &
193	& zgphi
194	REAL(KIND=wp), DIMENSION(1) :: zmsk
195	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
196	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim
197
198	LOGICAL :: ld_dailyav
199
200	!------------------------------------------------------------------------
201	! Local initialization
202	!------------------------------------------------------------------------
203	! Record and data counters
204	inrc = kt - kit000 + 2
205	ipro = prodatqc%npstp(inrc)
206
207	! Daily average types
208	ld_dailyav = .FALSE.
209	IF ( PRESENT(kdailyavtypes) ) THEN
210	idailyavtypes(:) = kdailyavtypes(:)
211	IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
212	ELSE
213	idailyavtypes(:) = -1
214	ENDIF
215
216	! Daily means are calculated for values over timesteps:
217	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
218	idayend = MOD( kt - kit000 + 1, kdaystp )
219
220	IF ( ld_dailyav ) THEN
221
222	! Initialize daily mean for first timestep of the day
223	IF ( idayend == 1 .OR. kt == 0 ) THEN
224	DO jk = 1, jpk
225	DO jj = 1, jpj
226	DO ji = 1, jpi
227	prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
228	END DO
229	END DO
230	END DO
231	ENDIF
232
233	DO jk = 1, jpk
234	DO jj = 1, jpj
235	DO ji = 1, jpi
236	! Increment field for computing daily mean
237	prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
238	& + pvar(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,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
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	& igrdi(2,2,ipro), &
263	& igrdj(2,2,ipro), &
264	& zglam(2,2,ipro), &
265	& zgphi(2,2,ipro), &
266	& zmask(2,2,kpk,ipro), &
267	& zint(2,2,kpk,ipro), &
268	& zgdept(2,2,kpk,ipro), &
269	& zgdepw(2,2,kpk,ipro) &
270	& )
271
272	IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) )
273
274	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
275	iobs = jobs - prodatqc%nprofup
276	igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
277	igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
278	igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
279	igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
280	igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
281	igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
282	igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
283	igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
284	END DO
285
286	! Initialise depth arrays
287	zgdept(:,:,:,:) = 0.0
288	zgdepw(:,:,:,:) = 0.0
289
290	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
291	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
292	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask )
293	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint )
294
295	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )
296	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )
297
298	IF ( prodatqc%lclim ) THEN
299	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )
300	ENDIF
301
302	! At the end of the day also get interpolated means
303	IF ( ld_dailyav .AND. idayend == 0 ) THEN
304
305	ALLOCATE( zinm(2,2,kpk,ipro) )
306
307	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
308	& prodatqc%vdmean(:,:,:,kvar), zinm )
309
310	ENDIF
311
312	! Return if no observations to process
313	! Has to be done after comm commands to ensure processors
314	! stay in sync
315	IF ( ipro == 0 ) RETURN
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
342	! Masked values are calculated later.
343	IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
344
345	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
346	& zglam(:,:,iobs), zgphi(:,:,iobs), &
347	& zmask(:,:,1,iobs), zweig1, zmsk )
348
349	ENDIF
350
351	IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
352
353	zobsk(:) = obfillflt
354
355	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
356
357	IF ( idayend == 0 ) THEN
358	! Daily averaged data
359
360	! vertically interpolate all 4 corners
361	ista = prodatqc%npvsta(jobs,kvar)
362	iend = prodatqc%npvend(jobs,kvar)
363	inum_obs = iend - ista + 1
364	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
365	IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
366
367	DO iin=1,2
368	DO ijn=1,2
369
370	IF ( k1dint == 1 ) THEN
371	CALL obs_int_z1d_spl( kpk, &
372	& zinm(iin,ijn,:,iobs), &
373	& zobs2k, zgdept(iin,ijn,:,iobs), &
374	& zmask(iin,ijn,:,iobs))
375
376	IF ( prodatqc%lclim ) THEN
377	CALL obs_int_z1d_spl( kpk, &
378	& zclim(iin,ijn,:,iobs), &
379	& zclm2k, zgdept(iin,ijn,:,iobs), &
380	& zmask(iin,ijn,:,iobs))
381	ENDIF
382
383	ENDIF
384
385	CALL obs_level_search(kpk, &
386	& zgdept(iin,ijn,:,iobs), &
387	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
388	& iv_indic)
389
390	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
391	& prodatqc%var(kvar)%vdep(ista:iend), &
392	& zinm(iin,ijn,:,iobs), &
393	& zobs2k, interp_corner(iin,ijn,:), &
394	& zgdept(iin,ijn,:,iobs), &
395	& zmask(iin,ijn,:,iobs))
396
397	IF ( prodatqc%lclim ) THEN
398	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
399	& prodatqc%var(kvar)%vdep(ista:iend), &
400	& zclim(iin,ijn,:,iobs), &
401	& zclm2k, interp_corner_clim(iin,ijn,:), &
402	& zgdept(iin,ijn,:,iobs), &
403	& zmask(iin,ijn,:,iobs))
404	ENDIF
405
406	ENDDO
407	ENDDO
408
409	ENDIF !idayend
410
411	ELSE
412
413	! Point data
414
415	! vertically interpolate all 4 corners
416	ista = prodatqc%npvsta(jobs,kvar)
417	iend = prodatqc%npvend(jobs,kvar)
418	inum_obs = iend - ista + 1
419	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
420	IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
421	DO iin=1,2
422	DO ijn=1,2
423
424	IF ( k1dint == 1 ) THEN
425	CALL obs_int_z1d_spl( kpk, &
426	& zint(iin,ijn,:,iobs),&
427	& zobs2k, zgdept(iin,ijn,:,iobs), &
428	& zmask(iin,ijn,:,iobs))
429
430	IF ( prodatqc%lclim ) THEN
431	CALL obs_int_z1d_spl( kpk, &
432	& zclim(iin,ijn,:,iobs),&
433	& zclm2k, zgdept(iin,ijn,:,iobs), &
434	& zmask(iin,ijn,:,iobs))
435	ENDIF
436
437	ENDIF
438
439	CALL obs_level_search(kpk, &
440	& zgdept(iin,ijn,:,iobs),&
441	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
442	& iv_indic)
443
444	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
445	& prodatqc%var(kvar)%vdep(ista:iend), &
446	& zint(iin,ijn,:,iobs), &
447	& zobs2k,interp_corner(iin,ijn,:), &
448	& zgdept(iin,ijn,:,iobs), &
449	& zmask(iin,ijn,:,iobs) )
450
451	IF ( prodatqc%lclim ) THEN
452	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
453	& prodatqc%var(kvar)%vdep(ista:iend), &
454	& zclim(iin,ijn,:,iobs), &
455	& zclm2k,interp_corner_clim(iin,ijn,:), &
456	& zgdept(iin,ijn,:,iobs), &
457	& zmask(iin,ijn,:,iobs) )
458	ENDIF
459
460	ENDDO
461	ENDDO
462
463	ENDIF
464
465	!-------------------------------------------------------------
466	! Compute the horizontal interpolation for every profile level
467	!-------------------------------------------------------------
468
469	DO ikn=1,inum_obs
470	iend=ista+ikn-1
471
472	zweig(:,:,1) = 0._wp
473
474	! This code forces the horizontal weights to be
475	! zero IF the observation is below the bottom of the
476	! corners of the interpolation nodes, Or if it is in
477	! the mask. This is important for observations near
478	! steep bathymetry
479	DO iin=1,2
480	DO ijn=1,2
481
482	depth_loop: DO ik=kpk,2,-1
483	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
484
485	zweig(iin,ijn,1) = &
486	& zweig1(iin,ijn,1) * &
487	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
488	& - prodatqc%var(kvar)%vdep(iend)),0._wp)
489
490	EXIT depth_loop
491
492	ENDIF
493
494	ENDDO depth_loop
495
496	ENDDO
497	ENDDO
498
499	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
500	& prodatqc%var(kvar)%vmod(iend:iend) )
501
502	IF ( prodatqc%lclim ) THEN
503	CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), &
504	& prodatqc%var(kvar)%vclm(iend:iend) )
505	ENDIF
506
507	! Set QC flag for any observations found below the bottom
508	! needed as the check here is more strict than that in obs_prep
509	IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
510
511	ENDDO
512
513	DEALLOCATE(interp_corner,iv_indic)
514	IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )
515
516	ENDIF
517
518	ENDDO
519
520	! Deallocate the data for interpolation
521	DEALLOCATE( &
522	& igrdi, &
523	& igrdj, &
524	& zglam, &
525	& zgphi, &
526	& zmask, &
527	& zint, &
528	& zgdept, &
529	& zgdepw &
530	& )
531
532	IF ( prodatqc%lclim ) DEALLOCATE( zclim )
533
534	! At the end of the day also get interpolated means
535	IF ( ld_dailyav .AND. idayend == 0 ) THEN
536	DEALLOCATE( zinm )
537	ENDIF
538
539	IF ( kvar == prodatqc%nvar ) THEN
540	prodatqc%nprofup = prodatqc%nprofup + ipro
541	ENDIF
542
543	END SUBROUTINE obs_prof_opt
544
545	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
546	& kit000, kdaystp, kvar, &
547	& psurf, pclim, psurfmask, &
548	& plam, pphi, &
549	& k2dint, ldnightav, plamscl, pphiscl, &
550	& lindegrees, kmeanstp )
551
552	!!-----------------------------------------------------------------------
553	!!
554	!! * ROUTINE obs_surf_opt *
555	!!
556	!! ** Purpose : Compute the model counterpart of surface
557	!! data by interpolating from the model grid to the
558	!! observation point.
559	!!
560	!! ** Method : Linearly interpolate to each observation point using
561	!! the model values at the corners of the surrounding grid box.
562	!!
563	!! The new model value is first computed at the obs (lon, lat) point.
564	!!
565	!! Several horizontal interpolation schemes are available:
566	!! - distance-weighted (great circle) (k2dint = 0)
567	!! - distance-weighted (small angle) (k2dint = 1)
568	!! - bilinear (geographical grid) (k2dint = 2)
569	!! - bilinear (quadrilateral grid) (k2dint = 3)
570	!! - polynomial (quadrilateral grid) (k2dint = 4)
571	!!
572	!! Two horizontal averaging schemes are also available:
573	!! - weighted radial footprint (k2dint = 5)
574	!! - weighted rectangular footprint (k2dint = 6)
575	!!
576	!!
577	!! ** Action :
578	!!
579	!! History :
580	!! ! 07-03 (A. Weaver)
581	!! ! 15-02 (M. Martin) Combined routine for surface types
582	!! ! 17-03 (M. Martin) Added horizontal averaging options
583	!! ! 20-08 (M. Martin) Added surface velocity options
584	!!-----------------------------------------------------------------------
585
586	!! * Modules used
587	USE obs_surf_def ! Definition of storage space for surface observations
588
589	IMPLICIT NONE
590
591	!! * Arguments
592	TYPE(obs_surf), INTENT(INOUT) :: &
593	& surfdataqc ! Subset of surface data passing QC
594	INTEGER, INTENT(IN) :: kt ! Time step
595	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
596	INTEGER, INTENT(IN) :: kpj
597	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
598	! (kit000-1 = restart time)
599	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
600	INTEGER, INTENT(IN) :: kvar ! Number of variables in surfdataqc
601	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
602	INTEGER, INTENT(IN), OPTIONAL :: &
603	kmeanstp ! Number of time steps for the time meaning
604	! Averaging is triggered if present and greater than one
605	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
606	& psurf, & ! Model surface field
607	& pclim, & ! Climatological surface field
608	& psurfmask ! Land-sea mask
609	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
610	& plam, & ! Model longitudes for variable
611	& pphi ! Model latitudes for variable
612	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
613	REAL(KIND=wp), INTENT(IN) :: &
614	& plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions
615	& pphiscl ! This is the full width (rather than half-width)
616	LOGICAL, INTENT(IN) :: &
617	& lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
618
619	!! * Local declarations
620	INTEGER :: ji
621	INTEGER :: jj
622	INTEGER :: jobs
623	INTEGER :: inrc
624	INTEGER :: isurf
625	INTEGER :: iobs
626	INTEGER :: imaxifp, imaxjfp
627	INTEGER :: imodi, imodj
628	INTEGER :: idayend
629	INTEGER :: imeanend
630	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
631	& igrdi, &
632	& igrdj, &
633	& igrdip1, &
634	& igrdjp1
635	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
636	& icount_night, &
637	& imask_night
638	REAL(wp) :: zlam
639	REAL(wp) :: zphi
640	REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
641	REAL(wp) :: zdaystp
642	REAL(wp) :: zmeanstp
643	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
644	& zweig, &
645	& zmask, &
646	& zsurf, &
647	& zsurfm, &
648	& zsurftmp, &
649	& zclim, &
650	& zglam, &
651	& zgphi, &
652	& zglamf, &
653	& zgphif
654
655	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
656	& zintmp, &
657	& zouttmp, &
658	& zmeanday ! to compute model sst in region of 24h daylight (pole)
659
660	LOGICAL :: l_timemean
661
662	!------------------------------------------------------------------------
663	! Local initialization
664	!------------------------------------------------------------------------
665	! Record and data counters
666	inrc = kt - kit000 + 2
667	isurf = surfdataqc%nsstp(inrc)
668
669	l_timemean = .FALSE.
670	IF ( PRESENT( kmeanstp ) ) THEN
671	IF ( kmeanstp > 1 ) l_timemean = .TRUE.
672	ENDIF
673
674	! Work out the maximum footprint size for the
675	! interpolation/averaging in model grid-points - has to be even.
676
677	CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
678
679	IF ( l_timemean ) THEN
680	! Initialize time mean for first timestep
681	imeanend = MOD( kt - kit000 + 1, kmeanstp )
682	IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend
683
684	! Added kt == 0 test to catch restart case
685	IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN
686	IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt
687	DO jj = 1, jpj
688	DO ji = 1, jpi
689	surfdataqc%vdmean(ji,jj,kvar) = 0.0
690	END DO
691	END DO
692	ENDIF
693
694	! On each time-step, increment the field for computing time mean
695	IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt
696	DO jj = 1, jpj
697	DO ji = 1, jpi
698	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
699	& + psurf(ji,jj)
700	END DO
701	END DO
702
703	! Compute the time mean at the end of time period
704	IF ( imeanend == 0 ) THEN
705	zmeanstp = 1.0 / REAL( kmeanstp )
706	IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp
707	DO jj = 1, jpj
708	DO ji = 1, jpi
709	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
710	& * zmeanstp
711	END DO
712	END DO
713	ENDIF
714	ENDIF !l_timemean
715
716
717	IF ( ldnightav ) THEN
718
719	! Initialize array for night mean
720	IF ( kt == 0 ) THEN
721	ALLOCATE ( icount_night(kpi,kpj) )
722	ALLOCATE ( imask_night(kpi,kpj) )
723	ALLOCATE ( zintmp(kpi,kpj) )
724	ALLOCATE ( zouttmp(kpi,kpj) )
725	ALLOCATE ( zmeanday(kpi,kpj) )
726	nday_qsr = -1 ! initialisation flag for nbc_dcy
727	ENDIF
728
729	! Night-time means are calculated for night-time values over timesteps:
730	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
731	idayend = MOD( kt - kit000 + 1, kdaystp )
732
733	! Initialize night-time mean for first timestep of the day
734	IF ( idayend == 1 .OR. kt == 0 ) THEN
735	DO jj = 1, jpj
736	DO ji = 1, jpi
737	surfdataqc%vdmean(ji,jj,kvar) = 0.0
738	zmeanday(ji,jj) = 0.0
739	icount_night(ji,jj) = 0
740	END DO
741	END DO
742	ENDIF
743
744	zintmp(:,:) = 0.0
745	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
746	imask_night(:,:) = INT( zouttmp(:,:) )
747
748	DO jj = 1, jpj
749	DO ji = 1, jpi
750	! Increment the temperature field for computing night mean and counter
751	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
752	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
753	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
754	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
755	END DO
756	END DO
757
758	! Compute the night-time mean at the end of the day
759	zdaystp = 1.0 / REAL( kdaystp )
760	IF ( idayend == 0 ) THEN
761	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
762	DO jj = 1, jpj
763	DO ji = 1, jpi
764	! Test if "no night" point
765	IF ( icount_night(ji,jj) > 0 ) THEN
766	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
767	& / REAL( icount_night(ji,jj) )
768	ELSE
769	!At locations where there is no night (e.g. poles),
770	! calculate daily mean instead of night-time mean.
771	surfdataqc%vdmean(ji,jj,kvar) = zmeanday(ji,jj) * zdaystp
772	ENDIF
773	END DO
774	END DO
775	ENDIF
776
777	ENDIF
778
779	! Get the data for interpolation
780	ALLOCATE( &
781	& zweig(imaxifp,imaxjfp,1), &
782	& igrdi(imaxifp,imaxjfp,isurf), &
783	& igrdj(imaxifp,imaxjfp,isurf), &
784	& zglam(imaxifp,imaxjfp,isurf), &
785	& zgphi(imaxifp,imaxjfp,isurf), &
786	& zmask(imaxifp,imaxjfp,isurf), &
787	& zsurf(imaxifp,imaxjfp,isurf), &
788	& zsurftmp(imaxifp,imaxjfp,isurf) &
789	& )
790
791	IF ( k2dint > 4 ) THEN
792	ALLOCATE( &
793	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
794	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
795	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
796	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
797	& )
798	ENDIF
799
800	IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
801
802	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
803	iobs = jobs - surfdataqc%nsurfup
804	DO ji = 0, imaxifp
805	imodi = surfdataqc%mi(jobs,kvar) - int(imaxifp/2) + ji - 1
806
807	!Deal with wrap around in longitude
808	IF ( imodi < 1 ) imodi = imodi + jpiglo
809	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
810
811	DO jj = 0, imaxjfp
812	imodj = surfdataqc%mj(jobs,kvar) - int(imaxjfp/2) + jj - 1
813	!If model values are out of the domain to the north/south then
814	!set them to be the edge of the domain
815	IF ( imodj < 1 ) imodj = 1
816	IF ( imodj > jpjglo ) imodj = jpjglo
817
818	IF ( k2dint > 4 ) THEN
819	igrdip1(ji+1,jj+1,iobs) = imodi
820	igrdjp1(ji+1,jj+1,iobs) = imodj
821	ENDIF
822
823	IF ( ji >= 1 .AND. jj >= 1 ) THEN
824	igrdi(ji,jj,iobs) = imodi
825	igrdj(ji,jj,iobs) = imodj
826	ENDIF
827
828	END DO
829	END DO
830	END DO
831
832	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
833	& igrdi, igrdj, plam, zglam )
834	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
835	& igrdi, igrdj, pphi, zgphi )
836	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
837	& igrdi, igrdj, psurfmask, zmask )
838
839	! At the end of the averaging period get interpolated means
840	IF ( l_timemean ) THEN
841	IF ( imeanend == 0 ) THEN
842	ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )
843	IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt
844	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
845	& igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm )
846	ENDIF
847	ELSE
848	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
849	& igrdi, igrdj, psurf, zsurf )
850	ENDIF
851
852	IF ( k2dint > 4 ) THEN
853	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
854	& igrdip1, igrdjp1, glamf, zglamf )
855	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
856	& igrdip1, igrdjp1, gphif, zgphif )
857	ENDIF
858
859	IF ( surfdataqc%lclim ) THEN
860	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
861	& igrdi, igrdj, pclim, zclim )
862	ENDIF
863
864	! At the end of the day get interpolated means
865	IF ( idayend == 0 .AND. ldnightav ) THEN
866
867	ALLOCATE( &
868	& zsurfm(imaxifp,imaxjfp,isurf) &
869	& )
870
871	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
872	& surfdataqc%vdmean(:,:,kvar), zsurfm )
873
874	ENDIF
875
876	! Loop over observations
877	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
878
879	iobs = jobs - surfdataqc%nsurfup
880
881	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
882
883	IF(lwp) THEN
884	WRITE(numout,*)
885	WRITE(numout,*) ' E R R O R : Observation', &
886	& ' time step is not consistent with the', &
887	& ' model time step'
888	WRITE(numout,*) ' ========='
889	WRITE(numout,*)
890	WRITE(numout,*) ' Record = ', jobs, &
891	& ' kt = ', kt, &
892	& ' mstp = ', surfdataqc%mstp(jobs), &
893	& ' ntyp = ', surfdataqc%ntyp(jobs)
894	ENDIF
895	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
896
897	ENDIF
898
899	zlam = surfdataqc%rlam(jobs)
900	zphi = surfdataqc%rphi(jobs)
901
902	IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN
903	! Night-time or N=kmeanstp timestep averaged data
904	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
905	ELSE
906	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
907	ENDIF
908
909	IF ( ( .NOT. l_timemean ) .OR. &
910	& ( l_timemean .AND. imeanend == 0) ) THEN
911	IF ( k2dint <= 4 ) THEN
912
913	! Get weights to interpolate the model value to the observation point
914	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
915	& zglam(:,:,iobs), zgphi(:,:,iobs), &
916	& zmask(:,:,iobs), zweig, zobsmask )
917
918	! Interpolate the model value to the observation point
919	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
920
921	IF ( surfdataqc%lclim ) THEN
922	CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
923	ENDIF
924
925
926	ELSE
927
928	! Get weights to average the model SLA to the observation footprint
929	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
930	& zglam(:,:,iobs), zgphi(:,:,iobs), &
931	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
932	& zmask(:,:,iobs), plamscl, pphiscl, &
933	& lindegrees, zweig, zobsmask )
934
935	! Average the model SST to the observation footprint
936	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
937	& zweig, zsurftmp(:,:,iobs), zext )
938
939	IF ( surfdataqc%lclim ) THEN
940	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
941	& zweig, zclim(:,:,iobs), zclm )
942	ENDIF
943
944	ENDIF
945
946	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
947	! ... Remove the MDT from the SSH at the observation point to get the SLA
948	surfdataqc%rext(jobs,1) = zext(1)
949	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
950	ELSE
951	surfdataqc%rmod(jobs,kvar) = zext(1)
952	ENDIF
953
954	IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
955
956	IF ( zext(1) == obfillflt ) THEN
957	! If the observation value is a fill value, set QC flag to bad
958	surfdataqc%nqc(jobs) = 4
959	ENDIF
960	ENDIF
961
962	IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
963
964	#if defined key_cice
965	IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN
966	! Convert radar freeboard to true freeboard (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum compared to snow (3.0e8 vs 2.4e8 m/s))
967	surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1)
968	surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)
969	! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
970	IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN
971	surfdataqc%nqc(jobs) = 4
972	ENDIF
973	! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
974	surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)rhow + surfdataqc%rext(jobs,2)rhos)/ &
975	(rhow - rhoi)
976	! Flag any negative ice thickness values as bad
977	IF (surfdataqc%robs(jobs,1) < 0.0) THEN
978	surfdataqc%nqc(jobs) = 4
979	ENDIF
980	ENDIF
981
982	IF ( zext(1) == obfillflt ) THEN
983	! If the observation value is a fill value, set QC flag to bad
984	surfdataqc%nqc(jobs) = 4
985	ENDIF
986	#endif defined key_cice
987
988	END DO
989
990	! Deallocate the data for interpolation
991	DEALLOCATE( &
992	& zweig, &
993	& igrdi, &
994	& igrdj, &
995	& zglam, &
996	& zgphi, &
997	& zmask, &
998	& zsurf, &
999	& zsurftmp &
1000	& )
1001
1002	IF ( k2dint > 4 ) THEN
1003	DEALLOCATE( &
1004	& zglamf, &
1005	& zgphif, &
1006	& igrdip1,&
1007	& igrdjp1 &
1008	& )
1009	ENDIF
1010
1011	IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
1012
1013	! At the end of the day also deallocate night-time mean array
1014	IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN
1015	DEALLOCATE( &
1016	& zsurfm &
1017	& )
1018	ENDIF
1019
1020	IF ( kvar == surfdataqc%nvar ) THEN
1021	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
1022	ENDIF
1023
1024	END SUBROUTINE obs_surf_opt
1025
1026	END MODULE obs_oper

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

Download in other formats: