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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 12380

Last change on this file since 12380 was 12380, checked in by dcarneir, 4 years ago
Updating obs_oper with obs_oper version in the apps trunk
File size: 38.9 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, psurf, pclim, psurfmask, &
547	& k2dint, ldnightav, plamscl, pphiscl, &
548	& lindegrees, kmeanstp )
549
550	!!-----------------------------------------------------------------------
551	!!
552	!! * ROUTINE obs_surf_opt *
553	!!
554	!! ** Purpose : Compute the model counterpart of surface
555	!! data by interpolating from the model grid to the
556	!! observation point.
557	!!
558	!! ** Method : Linearly interpolate to each observation point using
559	!! the model values at the corners of the surrounding grid box.
560	!!
561	!! The new model value is first computed at the obs (lon, lat) point.
562	!!
563	!! Several horizontal interpolation schemes are available:
564	!! - distance-weighted (great circle) (k2dint = 0)
565	!! - distance-weighted (small angle) (k2dint = 1)
566	!! - bilinear (geographical grid) (k2dint = 2)
567	!! - bilinear (quadrilateral grid) (k2dint = 3)
568	!! - polynomial (quadrilateral grid) (k2dint = 4)
569	!!
570	!! Two horizontal averaging schemes are also available:
571	!! - weighted radial footprint (k2dint = 5)
572	!! - weighted rectangular footprint (k2dint = 6)
573	!!
574	!!
575	!! ** Action :
576	!!
577	!! History :
578	!! ! 07-03 (A. Weaver)
579	!! ! 15-02 (M. Martin) Combined routine for surface types
580	!! ! 17-03 (M. Martin) Added horizontal averaging options
581	!!-----------------------------------------------------------------------
582
583	!! * Modules used
584	USE obs_surf_def ! Definition of storage space for surface observations
585
586	IMPLICIT NONE
587
588	!! * Arguments
589	TYPE(obs_surf), INTENT(INOUT) :: &
590	& surfdataqc ! Subset of surface data passing QC
591	INTEGER, INTENT(IN) :: kt ! Time step
592	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
593	INTEGER, INTENT(IN) :: kpj
594	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
595	! (kit000-1 = restart time)
596	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
597	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
598	INTEGER, INTENT(IN), OPTIONAL :: &
599	kmeanstp ! Number of time steps for the time meaning
600	! Averaging is triggered if present and greater than one
601	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
602	& psurf, & ! Model surface field
603	& pclim, & ! Climatological surface field
604	& psurfmask ! Land-sea mask
605	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
606	REAL(KIND=wp), INTENT(IN) :: &
607	& plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions
608	& pphiscl ! This is the full width (rather than half-width)
609	LOGICAL, INTENT(IN) :: &
610	& lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
611
612	!! * Local declarations
613	INTEGER :: ji
614	INTEGER :: jj
615	INTEGER :: jobs
616	INTEGER :: inrc
617	INTEGER :: isurf
618	INTEGER :: iobs
619	INTEGER :: imaxifp, imaxjfp
620	INTEGER :: imodi, imodj
621	INTEGER :: idayend
622	INTEGER :: imeanend
623	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
624	& igrdi, &
625	& igrdj, &
626	& igrdip1, &
627	& igrdjp1
628	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
629	& icount_night, &
630	& imask_night
631	REAL(wp) :: zlam
632	REAL(wp) :: zphi
633	REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
634	REAL(wp) :: zdaystp
635	REAL(wp) :: zmeanstp
636	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
637	& zweig, &
638	& zmask, &
639	& zsurf, &
640	& zsurfm, &
641	& zsurftmp, &
642	& zclim, &
643	& zglam, &
644	& zgphi, &
645	& zglamf, &
646	& zgphif
647
648	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
649	& zintmp, &
650	& zouttmp, &
651	& zmeanday ! to compute model sst in region of 24h daylight (pole)
652
653	LOGICAL :: l_timemean
654
655	!------------------------------------------------------------------------
656	! Local initialization
657	!------------------------------------------------------------------------
658	! Record and data counters
659	inrc = kt - kit000 + 2
660	isurf = surfdataqc%nsstp(inrc)
661
662	l_timemean = .FALSE.
663	IF ( PRESENT( kmeanstp ) ) THEN
664	IF ( kmeanstp > 1 ) l_timemean = .TRUE.
665	ENDIF
666
667	! Work out the maximum footprint size for the
668	! interpolation/averaging in model grid-points - has to be even.
669
670	CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
671
672
673	IF ( l_timemean ) THEN
674	! Initialize time mean for first timestep
675	imeanend = MOD( kt - kit000 + 1, kmeanstp )
676	IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend
677
678	! Added kt == 0 test to catch restart case
679	IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN
680	IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt
681	DO jj = 1, jpj
682	DO ji = 1, jpi
683	surfdataqc%vdmean(ji,jj) = 0.0
684	END DO
685	END DO
686	ENDIF
687
688	! On each time-step, increment the field for computing time mean
689	IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt
690	DO jj = 1, jpj
691	DO ji = 1, jpi
692	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
693	& + psurf(ji,jj)
694	END DO
695	END DO
696
697	! Compute the time mean at the end of time period
698	IF ( imeanend == 0 ) THEN
699	zmeanstp = 1.0 / REAL( kmeanstp )
700	IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp
701	DO jj = 1, jpj
702	DO ji = 1, jpi
703	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
704	& * zmeanstp
705	END DO
706	END DO
707	ENDIF
708	ENDIF !l_timemean
709
710
711	IF ( ldnightav ) THEN
712
713	! Initialize array for night mean
714	IF ( kt == 0 ) THEN
715	ALLOCATE ( icount_night(kpi,kpj) )
716	ALLOCATE ( imask_night(kpi,kpj) )
717	ALLOCATE ( zintmp(kpi,kpj) )
718	ALLOCATE ( zouttmp(kpi,kpj) )
719	ALLOCATE ( zmeanday(kpi,kpj) )
720	nday_qsr = -1 ! initialisation flag for nbc_dcy
721	ENDIF
722
723	! Night-time means are calculated for night-time values over timesteps:
724	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
725	idayend = MOD( kt - kit000 + 1, kdaystp )
726
727	! Initialize night-time mean for first timestep of the day
728	IF ( idayend == 1 .OR. kt == 0 ) THEN
729	DO jj = 1, jpj
730	DO ji = 1, jpi
731	surfdataqc%vdmean(ji,jj) = 0.0
732	zmeanday(ji,jj) = 0.0
733	icount_night(ji,jj) = 0
734	END DO
735	END DO
736	ENDIF
737
738	zintmp(:,:) = 0.0
739	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
740	imask_night(:,:) = INT( zouttmp(:,:) )
741
742	DO jj = 1, jpj
743	DO ji = 1, jpi
744	! Increment the temperature field for computing night mean and counter
745	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
746	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
747	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
748	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
749	END DO
750	END DO
751
752	! Compute the night-time mean at the end of the day
753	zdaystp = 1.0 / REAL( kdaystp )
754	IF ( idayend == 0 ) THEN
755	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
756	DO jj = 1, jpj
757	DO ji = 1, jpi
758	! Test if "no night" point
759	IF ( icount_night(ji,jj) > 0 ) THEN
760	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
761	& / REAL( icount_night(ji,jj) )
762	ELSE
763	!At locations where there is no night (e.g. poles),
764	! calculate daily mean instead of night-time mean.
765	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
766	ENDIF
767	END DO
768	END DO
769	ENDIF
770
771	ENDIF
772
773	! Get the data for interpolation
774
775	ALLOCATE( &
776	& zweig(imaxifp,imaxjfp,1), &
777	& igrdi(imaxifp,imaxjfp,isurf), &
778	& igrdj(imaxifp,imaxjfp,isurf), &
779	& zglam(imaxifp,imaxjfp,isurf), &
780	& zgphi(imaxifp,imaxjfp,isurf), &
781	& zmask(imaxifp,imaxjfp,isurf), &
782	& zsurf(imaxifp,imaxjfp,isurf), &
783	& zsurftmp(imaxifp,imaxjfp,isurf), &
784	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
785	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
786	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
787	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
788	& )
789
790	IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
791
792	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
793	iobs = jobs - surfdataqc%nsurfup
794	DO ji = 0, imaxifp
795	imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
796
797	!Deal with wrap around in longitude
798	IF ( imodi < 1 ) imodi = imodi + jpiglo
799	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
800
801	DO jj = 0, imaxjfp
802	imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
803	!If model values are out of the domain to the north/south then
804	!set them to be the edge of the domain
805	IF ( imodj < 1 ) imodj = 1
806	IF ( imodj > jpjglo ) imodj = jpjglo
807
808	igrdip1(ji+1,jj+1,iobs) = imodi
809	igrdjp1(ji+1,jj+1,iobs) = imodj
810
811	IF ( ji >= 1 .AND. jj >= 1 ) THEN
812	igrdi(ji,jj,iobs) = imodi
813	igrdj(ji,jj,iobs) = imodj
814	ENDIF
815
816	END DO
817	END DO
818	END DO
819
820	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
821	& igrdi, igrdj, glamt, zglam )
822	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
823	& igrdi, igrdj, gphit, zgphi )
824	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
825	& igrdi, igrdj, psurfmask, zmask )
826
827	! At the end of the averaging period get interpolated means
828	IF ( l_timemean ) THEN
829	IF ( imeanend == 0 ) THEN
830	ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )
831	IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt
832	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
833	& igrdi, igrdj, surfdataqc%vdmean(:,:), zsurfm )
834	ENDIF
835	ELSE
836	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
837	& igrdi, igrdj, psurf, zsurf )
838	ENDIF
839
840	IF ( k2dint > 4 ) THEN
841	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
842	& igrdip1, igrdjp1, glamf, zglamf )
843	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
844	& igrdip1, igrdjp1, gphif, zgphif )
845	ENDIF
846
847	IF ( surfdataqc%lclim ) THEN
848	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
849	& igrdi, igrdj, pclim, zclim )
850	ENDIF
851
852	! At the end of the day get interpolated means
853	IF ( idayend == 0 .AND. ldnightav ) THEN
854
855	ALLOCATE( &
856	& zsurfm(imaxifp,imaxjfp,isurf) &
857	& )
858
859	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
860	& surfdataqc%vdmean(:,:), zsurfm )
861
862	ENDIF
863
864	! Loop over observations
865	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
866
867	iobs = jobs - surfdataqc%nsurfup
868
869	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
870
871	IF(lwp) THEN
872	WRITE(numout,*)
873	WRITE(numout,*) ' E R R O R : Observation', &
874	& ' time step is not consistent with the', &
875	& ' model time step'
876	WRITE(numout,*) ' ========='
877	WRITE(numout,*)
878	WRITE(numout,*) ' Record = ', jobs, &
879	& ' kt = ', kt, &
880	& ' mstp = ', surfdataqc%mstp(jobs), &
881	& ' ntyp = ', surfdataqc%ntyp(jobs)
882	ENDIF
883	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
884
885	ENDIF
886
887	zlam = surfdataqc%rlam(jobs)
888	zphi = surfdataqc%rphi(jobs)
889
890	IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN
891	! Night-time or N=kmeanstp timestep averaged data
892	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
893	ELSE
894	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
895	ENDIF
896
897	IF ( ( .NOT. l_timemean ) .OR. &
898	& ( l_timemean .AND. imeanend == 0) ) THEN
899	IF ( k2dint <= 4 ) THEN
900
901	! Get weights to interpolate the model value to the observation point
902	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
903	& zglam(:,:,iobs), zgphi(:,:,iobs), &
904	& zmask(:,:,iobs), zweig, zobsmask )
905
906	! Interpolate the model value to the observation point
907	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
908
909	IF ( surfdataqc%lclim ) THEN
910	CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
911	ENDIF
912
913
914	ELSE
915
916	! Get weights to average the model SLA to the observation footprint
917	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
918	& zglam(:,:,iobs), zgphi(:,:,iobs), &
919	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
920	& zmask(:,:,iobs), plamscl, pphiscl, &
921	& lindegrees, zweig, zobsmask )
922
923	! Average the model SST to the observation footprint
924	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
925	& zweig, zsurftmp(:,:,iobs), zext )
926
927	IF ( surfdataqc%lclim ) THEN
928	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
929	& zweig, zclim(:,:,iobs), zclm )
930	ENDIF
931
932	ENDIF
933
934	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
935	! ... Remove the MDT from the SSH at the observation point to get the SLA
936	surfdataqc%rext(jobs,1) = zext(1)
937	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
938	ELSE
939	surfdataqc%rmod(jobs,1) = zext(1)
940	ENDIF
941
942	IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
943
944	IF ( zext(1) == obfillflt ) THEN
945	! If the observation value is a fill value, set QC flag to bad
946	surfdataqc%nqc(jobs) = 4
947	ENDIF
948	ENDIF
949
950	IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
951
952	IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN
953	! 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))
954	surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1)
955	surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)
956	! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
957	IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN
958	surfdataqc%nqc(jobs) = 4
959	ENDIF
960	! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
961	surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)rhow + surfdataqc%rext(jobs,2)rhos)/ &
962	(rhow - rhoi)
963	! Flag any negative ice thickness values as bad
964	IF (surfdataqc%robs(jobs,1) < 0.0) THEN
965	surfdataqc%nqc(jobs) = 4
966	ENDIF
967	ENDIF
968
969	IF ( zext(1) == obfillflt ) THEN
970	! If the observation value is a fill value, set QC flag to bad
971	surfdataqc%nqc(jobs) = 4
972	ENDIF
973
974	END DO
975
976	! Deallocate the data for interpolation
977	DEALLOCATE( &
978	& zweig, &
979	& igrdi, &
980	& igrdj, &
981	& zglam, &
982	& zgphi, &
983	& zmask, &
984	& zsurf, &
985	& zsurftmp, &
986	& zglamf, &
987	& zgphif, &
988	& igrdip1,&
989	& igrdjp1 &
990	& )
991
992	IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
993
994	! At the end of the day also deallocate night-time mean array
995	IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN
996	DEALLOCATE( &
997	& zsurfm &
998	& )
999	ENDIF
1000
1001	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
1002
1003	END SUBROUTINE obs_surf_opt
1004
1005	END MODULE obs_oper

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

Download in other formats: