Context Navigation

obs_oper.F90 @ 11200

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

Download in other formats: