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 NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/obs_oper.F90 @ 15248

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

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

Download in other formats: