Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

obs_oper.F90

Last change on this file was 15764, checked in by jcastill, 2 years ago
Update with the latest changes in branches/UKMO/dev_r5518_obs_oper_update@15400
File size: 39.5 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	ENDIF
383
384	CALL obs_level_search(kpk, &
385	& zgdept(iin,ijn,:,iobs), &
386	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
387	& iv_indic)
388
389	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
390	& prodatqc%var(kvar)%vdep(ista:iend), &
391	& zinm(iin,ijn,:,iobs), &
392	& zobs2k, interp_corner(iin,ijn,:), &
393	& zgdept(iin,ijn,:,iobs), &
394	& zmask(iin,ijn,:,iobs))
395
396	IF ( prodatqc%lclim ) THEN
397	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
398	& prodatqc%var(kvar)%vdep(ista:iend), &
399	& zclim(iin,ijn,:,iobs), &
400	& zclm2k, interp_corner_clim(iin,ijn,:), &
401	& zgdept(iin,ijn,:,iobs), &
402	& zmask(iin,ijn,:,iobs))
403	ENDIF
404	ENDDO
405	ENDDO
406
407	ENDIF !idayend
408
409	ELSE
410
411	! Point data
412
413	! vertically interpolate all 4 corners
414	ista = prodatqc%npvsta(jobs,kvar)
415	iend = prodatqc%npvend(jobs,kvar)
416	inum_obs = iend - ista + 1
417	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
418	IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
419	DO iin=1,2
420	DO ijn=1,2
421
422	IF ( k1dint == 1 ) THEN
423	CALL obs_int_z1d_spl( kpk, &
424	& zint(iin,ijn,:,iobs),&
425	& zobs2k, zgdept(iin,ijn,:,iobs), &
426	& zmask(iin,ijn,:,iobs))
427
428	IF ( prodatqc%lclim ) THEN
429	CALL obs_int_z1d_spl( kpk, &
430	& zclim(iin,ijn,:,iobs),&
431	& zclm2k, zgdept(iin,ijn,:,iobs), &
432	& zmask(iin,ijn,:,iobs))
433	ENDIF
434	ENDIF
435
436	CALL obs_level_search(kpk, &
437	& zgdept(iin,ijn,:,iobs),&
438	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
439	& iv_indic)
440
441	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
442	& prodatqc%var(kvar)%vdep(ista:iend), &
443	& zint(iin,ijn,:,iobs), &
444	& zobs2k,interp_corner(iin,ijn,:), &
445	& zgdept(iin,ijn,:,iobs), &
446	& zmask(iin,ijn,:,iobs) )
447
448	IF ( prodatqc%lclim ) THEN
449	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
450	& prodatqc%var(kvar)%vdep(ista:iend), &
451	& zclim(iin,ijn,:,iobs), &
452	& zclm2k,interp_corner_clim(iin,ijn,:), &
453	& zgdept(iin,ijn,:,iobs), &
454	& zmask(iin,ijn,:,iobs) )
455	ENDIF
456	ENDDO
457	ENDDO
458
459	ENDIF
460
461	!-------------------------------------------------------------
462	! Compute the horizontal interpolation for every profile level
463	!-------------------------------------------------------------
464
465	DO ikn=1,inum_obs
466	iend=ista+ikn-1
467
468	zweig(:,:,1) = 0._wp
469
470	! This code forces the horizontal weights to be
471	! zero IF the observation is below the bottom of the
472	! corners of the interpolation nodes, Or if it is in
473	! the mask. This is important for observations near
474	! steep bathymetry
475	DO iin=1,2
476	DO ijn=1,2
477
478	depth_loop: DO ik=kpk,2,-1
479	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
480
481	zweig(iin,ijn,1) = &
482	& zweig1(iin,ijn,1) * &
483	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
484	& - prodatqc%var(kvar)%vdep(iend)),0._wp)
485
486	EXIT depth_loop
487
488	ENDIF
489
490	ENDDO depth_loop
491
492	ENDDO
493	ENDDO
494
495	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
496	& prodatqc%var(kvar)%vmod(iend:iend) )
497
498	IF ( prodatqc%lclim ) THEN
499	CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), &
500	& prodatqc%var(kvar)%vclm(iend:iend) )
501	ENDIF
502
503	! Set QC flag for any observations found below the bottom
504	! needed as the check here is more strict than that in obs_prep
505	IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
506
507	ENDDO
508
509	DEALLOCATE(interp_corner,iv_indic)
510	IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )
511
512	ENDIF
513
514	ENDDO
515
516	! Deallocate the data for interpolation
517	DEALLOCATE( &
518	& igrdi, &
519	& igrdj, &
520	& zglam, &
521	& zgphi, &
522	& zmask, &
523	& zint, &
524	& zgdept, &
525	& zgdepw &
526	& )
527
528	IF ( prodatqc%lclim ) DEALLOCATE( zclim )
529
530	! At the end of the day also get interpolated means
531	IF ( ld_dailyav .AND. idayend == 0 ) THEN
532	DEALLOCATE( zinm )
533	ENDIF
534
535	IF ( kvar == prodatqc%nvar ) THEN
536	prodatqc%nprofup = prodatqc%nprofup + ipro
537	ENDIF
538
539	END SUBROUTINE obs_prof_opt
540
541	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
542	& kit000, kdaystp, kvar, &
543	& psurf, pclim, psurfmask, &
544	& plam, pphi, &
545	& k2dint, ldnightav, plamscl, pphiscl, &
546	& lindegrees, kmeanstp )
547
548	!!-----------------------------------------------------------------------
549	!!
550	!! * ROUTINE obs_surf_opt *
551	!!
552	!! ** Purpose : Compute the model counterpart of surface
553	!! data by interpolating from the model grid to the
554	!! observation point.
555	!!
556	!! ** Method : Linearly interpolate to each observation point using
557	!! the model values at the corners of the surrounding grid box.
558	!!
559	!! The new model value is first computed at the obs (lon, lat) point.
560	!!
561	!! Several horizontal interpolation schemes are available:
562	!! - distance-weighted (great circle) (k2dint = 0)
563	!! - distance-weighted (small angle) (k2dint = 1)
564	!! - bilinear (geographical grid) (k2dint = 2)
565	!! - bilinear (quadrilateral grid) (k2dint = 3)
566	!! - polynomial (quadrilateral grid) (k2dint = 4)
567	!!
568	!! Two horizontal averaging schemes are also available:
569	!! - weighted radial footprint (k2dint = 5)
570	!! - weighted rectangular footprint (k2dint = 6)
571	!!
572	!!
573	!! ** Action :
574	!!
575	!! History :
576	!! ! 07-03 (A. Weaver)
577	!! ! 15-02 (M. Martin) Combined routine for surface types
578	!! ! 17-03 (M. Martin) Added horizontal averaging options
579	!! ! 20-08 (M. Martin) Added surface velocity options
580	!!-----------------------------------------------------------------------
581
582	!! * Modules used
583	USE obs_surf_def ! Definition of storage space for surface observations
584
585	IMPLICIT NONE
586
587	!! * Arguments
588	TYPE(obs_surf), INTENT(INOUT) :: &
589	& surfdataqc ! Subset of surface data passing QC
590	INTEGER, INTENT(IN) :: kt ! Time step
591	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
592	INTEGER, INTENT(IN) :: kpj
593	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
594	! (kit000-1 = restart time)
595	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
596	INTEGER, INTENT(IN) :: kvar ! Number of variables in surfdataqc
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	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
606	& plam, & ! Model longitudes for variable
607	& pphi ! Model latitudes for variable
608	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
609	REAL(KIND=wp), INTENT(IN) :: &
610	& plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions
611	& pphiscl ! This is the full width (rather than half-width)
612	LOGICAL, INTENT(IN) :: &
613	& lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
614
615	!! * Local declarations
616	INTEGER :: ji
617	INTEGER :: jj
618	INTEGER :: jobs
619	INTEGER :: inrc
620	INTEGER :: isurf
621	INTEGER :: iobs
622	INTEGER :: imaxifp, imaxjfp
623	INTEGER :: imodi, imodj
624	INTEGER :: idayend
625	INTEGER :: imeanend
626	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
627	& igrdi, &
628	& igrdj, &
629	& igrdip1, &
630	& igrdjp1
631	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
632	& icount_night, &
633	& imask_night
634	REAL(wp) :: zlam
635	REAL(wp) :: zphi
636	REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
637	REAL(wp) :: zdaystp
638	REAL(wp) :: zmeanstp
639	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
640	& zweig, &
641	& zmask, &
642	& zsurf, &
643	& zsurfm, &
644	& zsurftmp, &
645	& zclim, &
646	& zglam, &
647	& zgphi, &
648	& zglamf, &
649	& zgphif
650
651	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
652	& zintmp, &
653	& zouttmp, &
654	& zmeanday ! to compute model sst in region of 24h daylight (pole)
655
656	LOGICAL :: l_timemean
657
658	!------------------------------------------------------------------------
659	! Local initialization
660	!------------------------------------------------------------------------
661	! Record and data counters
662	inrc = kt - kit000 + 2
663	isurf = surfdataqc%nsstp(inrc)
664
665	l_timemean = .FALSE.
666	IF ( PRESENT( kmeanstp ) ) THEN
667	IF ( kmeanstp > 1 ) l_timemean = .TRUE.
668	ENDIF
669
670	! Work out the maximum footprint size for the
671	! interpolation/averaging in model grid-points - has to be even.
672
673	CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
674
675	IF ( l_timemean ) THEN
676	! Initialize time mean for first timestep
677	imeanend = MOD( kt - kit000 + 1, kmeanstp )
678	IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend
679
680	! Added kt == 0 test to catch restart case
681	IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN
682	IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt
683	DO jj = 1, jpj
684	DO ji = 1, jpi
685	surfdataqc%vdmean(ji,jj,kvar) = 0.0
686	END DO
687	END DO
688	ENDIF
689
690	! On each time-step, increment the field for computing time mean
691	IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt
692	DO jj = 1, jpj
693	DO ji = 1, jpi
694	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
695	& + psurf(ji,jj)
696	END DO
697	END DO
698
699	! Compute the time mean at the end of time period
700	IF ( imeanend == 0 ) THEN
701	zmeanstp = 1.0 / REAL( kmeanstp )
702	IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp
703	DO jj = 1, jpj
704	DO ji = 1, jpi
705	surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
706	& * zmeanstp
707	END DO
708	END DO
709	ENDIF
710	ENDIF !l_timemean
711
712	IF ( ldnightav ) THEN
713
714	! Initialize array for night mean
715	IF ( kt == 0 ) THEN
716	ALLOCATE ( icount_night(kpi,kpj) )
717	ALLOCATE ( imask_night(kpi,kpj) )
718	ALLOCATE ( zintmp(kpi,kpj) )
719	ALLOCATE ( zouttmp(kpi,kpj) )
720	ALLOCATE ( zmeanday(kpi,kpj) )
721	nday_qsr = -1 ! initialisation flag for nbc_dcy
722	ENDIF
723
724	! Night-time means are calculated for night-time values over timesteps:
725	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
726	idayend = MOD( kt - kit000 + 1, kdaystp )
727
728	! Initialize night-time mean for first timestep of the day
729	IF ( idayend == 1 .OR. kt == 0 ) THEN
730	DO jj = 1, jpj
731	DO ji = 1, jpi
732	surfdataqc%vdmean(ji,jj) = 0.0
733	zmeanday(ji,jj) = 0.0
734	icount_night(ji,jj) = 0
735	END DO
736	END DO
737	ENDIF
738
739	zintmp(:,:) = 0.0
740	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
741	imask_night(:,:) = INT( zouttmp(:,:) )
742
743	DO jj = 1, jpj
744	DO ji = 1, jpi
745	! Increment the temperature field for computing night mean and counter
746	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
747	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
748	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
749	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
750	END DO
751	END DO
752
753	! Compute the night-time mean at the end of the day
754	zdaystp = 1.0 / REAL( kdaystp )
755	IF ( idayend == 0 ) THEN
756	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
757	DO jj = 1, jpj
758	DO ji = 1, jpi
759	! Test if "no night" point
760	IF ( icount_night(ji,jj) > 0 ) THEN
761	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
762	& / REAL( icount_night(ji,jj) )
763	ELSE
764	!At locations where there is no night (e.g. poles),
765	! calculate daily mean instead of night-time mean.
766	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
767	ENDIF
768	END DO
769	END DO
770	ENDIF
771
772	ENDIF
773
774	! Get the data for interpolation
775
776	ALLOCATE( &
777	& zweig(imaxifp,imaxjfp,1), &
778	& igrdi(imaxifp,imaxjfp,isurf), &
779	& igrdj(imaxifp,imaxjfp,isurf), &
780	& zglam(imaxifp,imaxjfp,isurf), &
781	& zgphi(imaxifp,imaxjfp,isurf), &
782	& zmask(imaxifp,imaxjfp,isurf), &
783	& zsurf(imaxifp,imaxjfp,isurf), &
784	& zsurftmp(imaxifp,imaxjfp,isurf) &
785	& )
786
787	IF ( k2dint > 4 ) THEN
788	ALLOCATE( &
789	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
790	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
791	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
792	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
793	& )
794	ENDIF
795
796	IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
797
798	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
799	iobs = jobs - surfdataqc%nsurfup
800	DO ji = 0, imaxifp
801	imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
802
803	!Deal with wrap around in longitude
804	IF ( imodi < 1 ) imodi = imodi + jpiglo
805	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
806
807	DO jj = 0, imaxjfp
808	imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
809	!If model values are out of the domain to the north/south then
810	!set them to be the edge of the domain
811	IF ( imodj < 1 ) imodj = 1
812	IF ( imodj > jpjglo ) imodj = jpjglo
813
814	IF ( k2dint > 4 ) THEN
815	igrdip1(ji+1,jj+1,iobs) = imodi
816	igrdjp1(ji+1,jj+1,iobs) = imodj
817	ENDIF
818
819	IF ( ji >= 1 .AND. jj >= 1 ) THEN
820	igrdi(ji,jj,iobs) = imodi
821	igrdj(ji,jj,iobs) = imodj
822	ENDIF
823
824	END DO
825	END DO
826	END DO
827
828	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
829	& igrdi, igrdj, glamt, zglam )
830	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
831	& igrdi, igrdj, gphit, zgphi )
832	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
833	& igrdi, igrdj, psurfmask, zmask )
834
835	! At the end of the averaging period get interpolated means
836	IF ( l_timemean ) THEN
837	IF ( imeanend == 0 ) THEN
838	ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )
839	IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt
840	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
841	& igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm )
842	ENDIF
843	ELSE
844	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
845	& igrdi, igrdj, psurf, zsurf )
846	ENDIF
847
848	IF ( k2dint > 4 ) THEN
849	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
850	& igrdip1, igrdjp1, glamf, zglamf )
851	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
852	& igrdip1, igrdjp1, gphif, zgphif )
853	ENDIF
854
855	IF ( surfdataqc%lclim ) THEN
856	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
857	& igrdi, igrdj, pclim, zclim )
858	ENDIF
859
860	! At the end of the day get interpolated means
861	IF ( idayend == 0 .AND. ldnightav ) THEN
862
863	ALLOCATE( &
864	& zsurfm(imaxifp,imaxjfp,isurf) &
865	& )
866
867	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
868	& surfdataqc%vdmean(:,:), zsurfm )
869
870	ENDIF
871
872	! Loop over observations
873	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
874
875	iobs = jobs - surfdataqc%nsurfup
876
877	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
878
879	IF(lwp) THEN
880	WRITE(numout,*)
881	WRITE(numout,*) ' E R R O R : Observation', &
882	& ' time step is not consistent with the', &
883	& ' model time step'
884	WRITE(numout,*) ' ========='
885	WRITE(numout,*)
886	WRITE(numout,*) ' Record = ', jobs, &
887	& ' kt = ', kt, &
888	& ' mstp = ', surfdataqc%mstp(jobs), &
889	& ' ntyp = ', surfdataqc%ntyp(jobs)
890	ENDIF
891	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
892
893	ENDIF
894
895	zlam = surfdataqc%rlam(jobs)
896	zphi = surfdataqc%rphi(jobs)
897
898	IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN
899	! Night-time or N=kmeanstp timestep averaged data
900	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
901	ELSE
902	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
903	ENDIF
904
905	IF ( ( .NOT. l_timemean ) .OR. &
906	& ( l_timemean .AND. imeanend == 0) ) THEN
907	IF ( k2dint <= 4 ) THEN
908
909	! Get weights to interpolate the model value to the observation point
910	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
911	& zglam(:,:,iobs), zgphi(:,:,iobs), &
912	& zmask(:,:,iobs), zweig, zobsmask )
913
914	! Interpolate the model value to the observation point
915	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
916
917	IF ( surfdataqc%lclim ) THEN
918	CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
919	ENDIF
920
921
922	ELSE
923
924	! Get weights to average the model SLA to the observation footprint
925	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
926	& zglam(:,:,iobs), zgphi(:,:,iobs), &
927	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
928	& zmask(:,:,iobs), plamscl, pphiscl, &
929	& lindegrees, zweig, zobsmask )
930
931	! Average the model SST to the observation footprint
932	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
933	& zweig, zsurftmp(:,:,iobs), zext )
934
935	IF ( surfdataqc%lclim ) THEN
936	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
937	& zweig, zclim(:,:,iobs), zclm )
938	ENDIF
939
940	ENDIF
941
942	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
943	! ... Remove the MDT from the SSH at the observation point to get the SLA
944	surfdataqc%rext(jobs,1) = zext(1)
945	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
946	ELSE
947	surfdataqc%rmod(jobs,kvar) = zext(1)
948	ENDIF
949
950	IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
951
952	IF ( zext(1) == obfillflt ) THEN
953	! If the observation value is a fill value, set QC flag to bad
954	surfdataqc%nqc(jobs) = 4
955	ENDIF
956	ENDIF
957
958	IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
959
960	#if defined key_cice
961	IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN
962	! 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))
963	surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1)
964	surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)
965	! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
966	IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN
967	surfdataqc%nqc(jobs) = 4
968	ENDIF
969	! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
970	surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)rhow + surfdataqc%rext(jobs,2)rhos)/ &
971	(rhow - rhoi)
972	ENDIF
973
974	IF ( zext(1) == obfillflt ) THEN
975	! If the observation value is a fill value, set QC flag to bad
976	surfdataqc%nqc(jobs) = 4
977	ENDIF
978	#endif defined key_cice
979
980	END DO
981
982	! Deallocate the data for interpolation
983	DEALLOCATE( &
984	& zweig, &
985	& igrdi, &
986	& igrdj, &
987	& zglam, &
988	& zgphi, &
989	& zmask, &
990	& zsurf, &
991	& zsurftmp &
992	& )
993
994	IF ( k2dint > 4 ) THEN
995	DEALLOCATE( &
996	& zglamf, &
997	& zgphif, &
998	& igrdip1,&
999	& igrdjp1 &
1000	& )
1001	ENDIF
1002
1003	IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
1004
1005	! At the end of the day also deallocate night-time mean array
1006	IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN
1007	DEALLOCATE( &
1008	& zsurfm &
1009	& )
1010	ENDIF
1011
1012	IF ( kvar == surfdataqc%nvar ) THEN
1013	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
1014	ENDIF
1015
1016	END SUBROUTINE obs_surf_opt
1017
1018	END MODULE obs_oper

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

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

Context Navigation

source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

Download in other formats: