Context Navigation

obs_oper.F90 @ 15487

Last change on this file since 15487 was 15487, checked in by dford, 3 years ago
Remove some comments and add some spaces.
File size: 40.0 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: ', &
681	& kt, ' with weight: ', zmeanstp
682	DO jj = 1, jpj
683	DO ji = 1, jpi
684	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
685	& * zmeanstp
686	END DO
687	END DO
688	ENDIF
689	ENDIF
690
691	IF ( ldnightav ) THEN
692
693	! Initialize array for night mean
694	IF ( kt == 0 ) THEN
695	ALLOCATE ( icount_night(kpi,kpj) )
696	ALLOCATE ( imask_night(kpi,kpj) )
697	ALLOCATE ( zintmp(kpi,kpj) )
698	ALLOCATE ( zouttmp(kpi,kpj) )
699	ALLOCATE ( zmeanday(kpi,kpj) )
700	nday_qsr = -1 ! initialisation flag for nbc_dcy
701	ENDIF
702
703	! Night-time means are calculated for night-time values over timesteps:
704	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
705	idayend = MOD( kt - kit000 + 1, kdaystp )
706
707	! Initialize night-time mean for first timestep of the day
708	IF ( idayend == 1 .OR. kt == 0 ) THEN
709	DO jj = 1, jpj
710	DO ji = 1, jpi
711	surfdataqc%vdmean(ji,jj,:) = 0.0
712	zmeanday(ji,jj) = 0.0
713	icount_night(ji,jj) = 0
714	END DO
715	END DO
716	ENDIF
717
718	zintmp(:,:) = 0.0
719	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
720	imask_night(:,:) = INT( zouttmp(:,:) )
721
722	DO jj = 1, jpj
723	DO ji = 1, jpi
724	! Increment the model field for computing night mean and counter
725	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
726	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
727	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
728	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
729	END DO
730	END DO
731
732	! Compute the night-time mean at the end of the day
733	zdaystp = 1.0 / REAL( kdaystp )
734	IF ( idayend == 0 ) THEN
735	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ', kt
736	DO jj = 1, jpj
737	DO ji = 1, jpi
738	! Test if "no night" point
739	IF ( icount_night(ji,jj) > 0 ) THEN
740	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
741	& / REAL( icount_night(ji,jj) )
742	ELSE
743	!At locations where there is no night (e.g. poles),
744	! calculate daily mean instead of night-time mean.
745	surfdataqc%vdmean(ji,jj,kvar) = zmeanday(ji,jj) * zdaystp
746	ENDIF
747	END DO
748	END DO
749	ENDIF
750
751	ENDIF
752
753	! Get the data for interpolation
754
755	ALLOCATE( &
756	& zweig(imaxifp,imaxjfp,1), &
757	& igrdi(imaxifp,imaxjfp,isurf), &
758	& igrdj(imaxifp,imaxjfp,isurf), &
759	& zglam(imaxifp,imaxjfp,isurf), &
760	& zgphi(imaxifp,imaxjfp,isurf), &
761	& zmask(imaxifp,imaxjfp,isurf), &
762	& zsurf(imaxifp,imaxjfp,isurf), &
763	& zsurftmp(imaxifp,imaxjfp,isurf) &
764	& )
765
766	IF ( k2dint > 4 ) THEN
767	ALLOCATE( &
768	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
769	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
770	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
771	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
772	& )
773	ENDIF
774
775	IF ( ldclim ) THEN
776	ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
777	ENDIF
778
779	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
780	iobs = jobs - surfdataqc%nsurfup
781	DO ji = 0, imaxifp
782	imodi = surfdataqc%mi(jobs,kvar) - int(imaxifp/2) + ji - 1
783	!
784	!Deal with wrap around in longitude
785	IF ( imodi < 1 ) imodi = imodi + jpiglo
786	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
787	!
788	DO jj = 0, imaxjfp
789	imodj = surfdataqc%mj(jobs,kvar) - int(imaxjfp/2) + jj - 1
790	!If model values are out of the domain to the north/south then
791	!set them to be the edge of the domain
792	IF ( imodj < 1 ) imodj = 1
793	IF ( imodj > jpjglo ) imodj = jpjglo
794	!
795	IF ( k2dint > 4 ) THEN
796	igrdip1(ji+1,jj+1,iobs) = imodi
797	igrdjp1(ji+1,jj+1,iobs) = imodj
798	ENDIF
799	!
800	IF ( ji >= 1 .AND. jj >= 1 ) THEN
801	igrdi(ji,jj,iobs) = imodi
802	igrdj(ji,jj,iobs) = imodj
803	ENDIF
804	!
805	END DO
806	END DO
807	END DO
808
809	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
810	& igrdi, igrdj, glamt, zglam )
811	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
812	& igrdi, igrdj, gphit, zgphi )
813	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
814	& igrdi, igrdj, psurfmask, zmask )
815
816	! At the end of the averaging period get interpolated means
817	IF ( ldtime_mean ) THEN
818	IF ( imeanend == 0 ) THEN
819	ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )
820	IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ', kt
821	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
822	& igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm )
823	ENDIF
824	ELSE
825	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
826	& igrdi, igrdj, psurf, zsurf )
827	ENDIF
828
829	IF ( k2dint > 4 ) THEN
830	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
831	& igrdip1, igrdjp1, glamf, zglamf )
832	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
833	& igrdip1, igrdjp1, gphif, zgphif )
834	ENDIF
835
836	IF ( ldclim ) THEN
837	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
838	& igrdi, igrdj, pclim, zclim )
839	ENDIF
840
841	! At the end of the day get interpolated means
842	IF ( idayend == 0 .AND. ldnightav ) THEN
843
844	ALLOCATE( &
845	& zsurfm(imaxifp,imaxjfp,isurf) &
846	& )
847
848	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
849	& surfdataqc%vdmean(:,:,kvar), zsurfm )
850
851	ENDIF
852
853	! Loop over observations
854	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
855
856	iobs = jobs - surfdataqc%nsurfup
857
858	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
859
860	IF(lwp) THEN
861	WRITE(numout,*)
862	WRITE(numout,*) ' E R R O R : Observation', &
863	& ' time step is not consistent with the', &
864	& ' model time step'
865	WRITE(numout,*) ' ========='
866	WRITE(numout,*)
867	WRITE(numout,*) ' Record = ', jobs, &
868	& ' kt = ', kt, &
869	& ' mstp = ', surfdataqc%mstp(jobs), &
870	& ' ntyp = ', surfdataqc%ntyp(jobs)
871	ENDIF
872	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
873
874	ENDIF
875
876	zlam = surfdataqc%rlam(jobs)
877	zphi = surfdataqc%rphi(jobs)
878
879	IF ( ( ldnightav .AND. idayend == 0 ) .OR. (ldtime_mean .AND. imeanend == 0) ) THEN
880	! Night-time or N=kmeanstp timestep averaged data
881	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
882	ELSE
883	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
884	ENDIF
885
886	IF ( ( .NOT. ldtime_mean ) .OR. ( ldtime_mean .AND. imeanend == 0) ) THEN
887
888	IF ( k2dint <= 4 ) THEN
889
890	! Get weights to interpolate the model value to the observation point
891	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
892	& zglam(:,:,iobs), zgphi(:,:,iobs), &
893	& zmask(:,:,iobs), zweig, zobsmask )
894
895	! Interpolate the model value to the observation point
896	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
897
898	IF ( ldclim ) THEN
899	CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
900	ENDIF
901
902	ELSE
903
904	! Get weights to average the model field to the observation footprint
905	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
906	& zglam(:,:,iobs), zgphi(:,:,iobs), &
907	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
908	& zmask(:,:,iobs), plamscl, pphiscl, &
909	& lindegrees, zweig, zobsmask )
910
911	! Average the model field to the observation footprint
912	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
913	& zweig, zsurftmp(:,:,iobs), zext )
914
915	IF ( ldclim ) THEN
916	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
917	& zweig, zclim(:,:,iobs), zclm )
918	ENDIF
919
920	ENDIF
921
922	IF ( TRIM(surfdataqc%cvars(kvar)) == cobsname_sla .AND. PRESENT(kssh) .AND. PRESENT(kmdt) ) THEN
923	! ... Remove the MDT from the SSH at the observation point to get the SLA
924	surfdataqc%radd(jobs,kssh,kvar) = zext(1)
925	surfdataqc%rmod(jobs,kvar) = surfdataqc%radd(jobs,kssh,kvar) - surfdataqc%rext(jobs,kmdt)
926	#if defined key_si3 \|\| defined key_cice
927	ELSE IF ( TRIM(surfdataqc%cvars(kvar)) == cobsname_fbd .AND. PRESENT(kfbd) .AND. PRESENT(ksnow) ) THEN
928	surfdataqc%rmod(jobs,kvar) = zext(1)
929	! Convert radar freeboard to true freeboard
930	! (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum
931	! compared to snow (3.0e8 vs 2.4e8 m/s))
932	surfdataqc%radd(jobs,kfbd,kvar) = surfdataqc%robs(jobs,kvar)
933	surfdataqc%robs(jobs,kvar) = surfdataqc%radd(jobs,kfbd,kvar) + 0.25*surfdataqc%rext(jobs,ksnow)
934	! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
935	IF ((surfdataqc%robs(jobs,kvar) < -0.3) .OR. (surfdataqc%robs(jobs,kvar) > 3.0)) THEN
936	surfdataqc%nqc(jobs) = 4
937	ENDIF
938	! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
939	surfdataqc%robs(jobs,kvar) = (surfdataqc%robs(jobs,kvar)rhow + surfdataqc%rext(jobs,ksnow)rhos)/ &
940	& (rhow - rhoi)
941	#endif
942	ELSE
943	surfdataqc%rmod(jobs,kvar) = zext(1)
944	ENDIF
945
946	IF ( ldclim ) THEN
947	surfdataqc%radd(jobs,kclim,kvar) = zclm(1)
948	ENDIF
949
950	IF ( zext(1) == obfillflt ) THEN
951	! If the observation value is a fill value, set QC flag to bad
952	surfdataqc%nqc(jobs) = 4
953	ENDIF
954
955	ENDIF
956
957	END DO
958
959	! Deallocate the data for interpolation
960	DEALLOCATE( &
961	& zweig, &
962	& igrdi, &
963	& igrdj, &
964	& zglam, &
965	& zgphi, &
966	& zmask, &
967	& zsurf, &
968	& zsurftmp &
969	& )
970
971	IF ( k2dint > 4 ) THEN
972	DEALLOCATE( &
973	& zglamf, &
974	& zgphif, &
975	& igrdip1,&
976	& igrdjp1 &
977	& )
978	ENDIF
979
980	IF ( ldclim ) THEN
981	DEALLOCATE( zclim )
982	ENDIF
983
984	! At the end of the day also deallocate time mean array
985	IF ( ( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. ldtime_mean ) ) THEN
986	DEALLOCATE( &
987	& zsurfm &
988	& )
989	ENDIF
990	!
991	IF ( kvar == surfdataqc%nvar ) THEN
992	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
993	ENDIF
994	!
995	END SUBROUTINE obs_surf_opt
996
997	!!======================================================================
998	END MODULE obs_oper

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

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

Context Navigation

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

Download in other formats: