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 @ 15227

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

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

Download in other formats: