Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obs_oper.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago
Mixed precision version, tested up to 30 years on ORCA2.
Property svn:keywords set to `Id`
File size: 30.6 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	!
22	USE par_kind , ONLY : wp ! Precision variables
23	USE in_out_manager ! I/O manager
24
25	IMPLICIT NONE
26	PRIVATE
27
28	PUBLIC obs_prof_opt !: Compute the model counterpart of profile obs
29	PUBLIC obs_surf_opt !: Compute the model counterpart of surface obs
30
31	INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 !: Max number of daily avgd obs types
32
33	!! * Substitutions
34	# include "do_loop_substitute.h90"
35	# include "single_precision_substitute.h90"
36
37	!!----------------------------------------------------------------------
38	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
39	!! $Id$
40	!! Software governed by the CeCILL license (see ./LICENSE)
41	!!----------------------------------------------------------------------
42	CONTAINS
43
44	SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
45	& kit000, kdaystp, kvar, &
46	& pvar, 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 ! Number of variables in prodatqc
110	REAL(KIND=dp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pvar ! Model field
111	REAL(KIND=dp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pmask ! Land-sea mask
112	REAL(KIND=dp) , INTENT(in ), DIMENSION(kpi,kpj) :: plam ! Model longitude
113	REAL(KIND=dp) , INTENT(in ), DIMENSION(kpi,kpj) :: pphi ! Model latitudes
114	REAL(KIND=dp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pgdept, pgdepw ! depth of T and W levels
115	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: kdailyavtypes ! Types for daily averages
116
117	!! * Local declarations
118	INTEGER :: ji
119	INTEGER :: jj
120	INTEGER :: jk
121	INTEGER :: jobs
122	INTEGER :: inrc
123	INTEGER :: ipro
124	INTEGER :: idayend
125	INTEGER :: ista
126	INTEGER :: iend
127	INTEGER :: iobs
128	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
129	INTEGER :: inum_obs
130	INTEGER, DIMENSION(imaxavtypes) :: &
131	& idailyavtypes
132	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
133	& igrdi, &
134	& igrdj
135	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
136
137	REAL(KIND=dp) :: zlam
138	REAL(KIND=dp) :: zphi
139	REAL(KIND=dp) :: zdaystp
140	REAL(KIND=dp), DIMENSION(kpk) :: &
141	& zobsk, &
142	& zobs2k
143	REAL(KIND=dp), DIMENSION(2,2,1) :: &
144	& zweig1, &
145	& zweig
146	REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
147	& zmask, &
148	& zint, &
149	& zinm, &
150	& zgdept, &
151	& zgdepw
152	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: &
153	& zglam, &
154	& zgphi
155	REAL(KIND=dp), DIMENSION(1) :: zmsk
156	REAL(KIND=dp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
157
158	LOGICAL :: ld_dailyav
159
160	!------------------------------------------------------------------------
161	! Local initialization
162	!------------------------------------------------------------------------
163	! Record and data counters
164	inrc = kt - kit000 + 2
165	ipro = prodatqc%npstp(inrc)
166
167	! Daily average types
168	ld_dailyav = .FALSE.
169	IF ( PRESENT(kdailyavtypes) ) THEN
170	idailyavtypes(:) = kdailyavtypes(:)
171	IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
172	ELSE
173	idailyavtypes(:) = -1
174	ENDIF
175
176	! Daily means are calculated for values over timesteps:
177	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
178	idayend = MOD( kt - kit000 + 1, kdaystp )
179
180	IF ( ld_dailyav ) THEN
181
182	! Initialize daily mean for first timestep of the day
183	IF ( idayend == 1 .OR. kt == 0 ) THEN
184	DO_3D( 1, 1, 1, 1, 1, jpk )
185	prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
186	END_3D
187	ENDIF
188
189	DO_3D( 1, 1, 1, 1, 1, jpk )
190	! Increment field 1 for computing daily mean
191	prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
192	& + pvar(ji,jj,jk)
193	END_3D
194
195	! Compute the daily mean at the end of day
196	zdaystp = 1.0 / REAL( kdaystp )
197	IF ( idayend == 0 ) THEN
198	IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
199	CALL FLUSH(numout)
200	DO_3D( 1, 1, 1, 1, 1, jpk )
201	prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
202	& * zdaystp
203	END_3D
204	ENDIF
205
206	ENDIF
207
208	! Get the data for interpolation
209	ALLOCATE( &
210	& igrdi(2,2,ipro), &
211	& igrdj(2,2,ipro), &
212	& zglam(2,2,ipro), &
213	& zgphi(2,2,ipro), &
214	& zmask(2,2,kpk,ipro), &
215	& zint(2,2,kpk,ipro), &
216	& zgdept(2,2,kpk,ipro), &
217	& zgdepw(2,2,kpk,ipro) &
218	& )
219
220	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
221	iobs = jobs - prodatqc%nprofup
222	igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
223	igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
224	igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
225	igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
226	igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
227	igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
228	igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
229	igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
230	END DO
231
232	! Initialise depth arrays
233	zgdept(:,:,:,:) = 0.0
234	zgdepw(:,:,:,:) = 0.0
235
236	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
237	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
238	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask )
239	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint )
240
241	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )
242	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )
243
244	! At the end of the day also get interpolated means
245	IF ( ld_dailyav .AND. idayend == 0 ) THEN
246
247	ALLOCATE( zinm(2,2,kpk,ipro) )
248
249	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
250	& prodatqc%vdmean(:,:,:,kvar), zinm )
251
252	ENDIF
253
254	! Return if no observations to process
255	! Has to be done after comm commands to ensure processors
256	! stay in sync
257	IF ( ipro == 0 ) RETURN
258
259	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
260
261	iobs = jobs - prodatqc%nprofup
262
263	IF ( kt /= prodatqc%mstp(jobs) ) THEN
264
265	IF(lwp) THEN
266	WRITE(numout,*)
267	WRITE(numout,*) ' E R R O R : Observation', &
268	& ' time step is not consistent with the', &
269	& ' model time step'
270	WRITE(numout,*) ' ========='
271	WRITE(numout,*)
272	WRITE(numout,*) ' Record = ', jobs, &
273	& ' kt = ', kt, &
274	& ' mstp = ', prodatqc%mstp(jobs), &
275	& ' ntyp = ', prodatqc%ntyp(jobs)
276	ENDIF
277	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
278	ENDIF
279
280	zlam = prodatqc%rlam(jobs)
281	zphi = prodatqc%rphi(jobs)
282
283	! Horizontal weights
284	! Masked values are calculated later.
285	IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
286
287	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
288	& zglam(:,:,iobs), zgphi(:,:,iobs), &
289	& zmask(:,:,1,iobs), zweig1, zmsk )
290
291	ENDIF
292
293	IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
294
295	zobsk(:) = obfillflt
296
297	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
298
299	IF ( idayend == 0 ) THEN
300	! Daily averaged data
301
302	! vertically interpolate all 4 corners
303	ista = prodatqc%npvsta(jobs,kvar)
304	iend = prodatqc%npvend(jobs,kvar)
305	inum_obs = iend - ista + 1
306	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
307
308	DO iin=1,2
309	DO ijn=1,2
310
311	IF ( k1dint == 1 ) THEN
312	CALL obs_int_z1d_spl( kpk, &
313	& zinm(iin,ijn,:,iobs), &
314	& zobs2k, zgdept(iin,ijn,:,iobs), &
315	& zmask(iin,ijn,:,iobs))
316	ENDIF
317
318	CALL obs_level_search(kpk, &
319	& zgdept(iin,ijn,:,iobs), &
320	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
321	& iv_indic)
322
323	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
324	& prodatqc%var(kvar)%vdep(ista:iend), &
325	& zinm(iin,ijn,:,iobs), &
326	& zobs2k, interp_corner(iin,ijn,:), &
327	& zgdept(iin,ijn,:,iobs), &
328	& zmask(iin,ijn,:,iobs))
329
330	ENDDO
331	ENDDO
332
333	ENDIF !idayend
334
335	ELSE
336
337	! Point data
338
339	! vertically interpolate all 4 corners
340	ista = prodatqc%npvsta(jobs,kvar)
341	iend = prodatqc%npvend(jobs,kvar)
342	inum_obs = iend - ista + 1
343	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
344	DO iin=1,2
345	DO ijn=1,2
346
347	IF ( k1dint == 1 ) THEN
348	CALL obs_int_z1d_spl( kpk, &
349	& zint(iin,ijn,:,iobs),&
350	& zobs2k, zgdept(iin,ijn,:,iobs), &
351	& zmask(iin,ijn,:,iobs))
352
353	ENDIF
354
355	CALL obs_level_search(kpk, &
356	& zgdept(iin,ijn,:,iobs),&
357	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
358	& iv_indic)
359
360	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
361	& prodatqc%var(kvar)%vdep(ista:iend), &
362	& zint(iin,ijn,:,iobs), &
363	& zobs2k,interp_corner(iin,ijn,:), &
364	& zgdept(iin,ijn,:,iobs), &
365	& zmask(iin,ijn,:,iobs) )
366
367	ENDDO
368	ENDDO
369
370	ENDIF
371
372	!-------------------------------------------------------------
373	! Compute the horizontal interpolation for every profile level
374	!-------------------------------------------------------------
375
376	DO ikn=1,inum_obs
377	iend=ista+ikn-1
378
379	zweig(:,:,1) = 0._wp
380
381	! This code forces the horizontal weights to be
382	! zero IF the observation is below the bottom of the
383	! corners of the interpolation nodes, Or if it is in
384	! the mask. This is important for observations near
385	! steep bathymetry
386	DO iin=1,2
387	DO ijn=1,2
388
389	depth_loop: DO ik=kpk,2,-1
390	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
391
392	zweig(iin,ijn,1) = &
393	& zweig1(iin,ijn,1) * &
394	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
395	& - prodatqc%var(kvar)%vdep(iend)),0._wp)
396
397	EXIT depth_loop
398
399	ENDIF
400
401	ENDDO depth_loop
402
403	ENDDO
404	ENDDO
405
406	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
407	& prodatqc%var(kvar)%vmod(iend:iend) )
408
409	! Set QC flag for any observations found below the bottom
410	! needed as the check here is more strict than that in obs_prep
411	IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
412
413	ENDDO
414
415	DEALLOCATE(interp_corner,iv_indic)
416
417	ENDIF
418
419	ENDDO
420
421	! Deallocate the data for interpolation
422	DEALLOCATE( &
423	& igrdi, &
424	& igrdj, &
425	& zglam, &
426	& zgphi, &
427	& zmask, &
428	& zint, &
429	& zgdept, &
430	& zgdepw &
431	& )
432
433	! At the end of the day also get interpolated means
434	IF ( ld_dailyav .AND. idayend == 0 ) THEN
435	DEALLOCATE( zinm )
436	ENDIF
437
438	IF ( kvar == prodatqc%nvar ) THEN
439	prodatqc%nprofup = prodatqc%nprofup + ipro
440	ENDIF
441
442	END SUBROUTINE obs_prof_opt
443
444	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
445	& kit000, kdaystp, psurf, psurfmask, &
446	& k2dint, ldnightav, plamscl, pphiscl, &
447	& lindegrees )
448
449	!!-----------------------------------------------------------------------
450	!!
451	!! * ROUTINE obs_surf_opt *
452	!!
453	!! ** Purpose : Compute the model counterpart of surface
454	!! data by interpolating from the model grid to the
455	!! observation point.
456	!!
457	!! ** Method : Linearly interpolate to each observation point using
458	!! the model values at the corners of the surrounding grid box.
459	!!
460	!! The new model value is first computed at the obs (lon, lat) point.
461	!!
462	!! Several horizontal interpolation schemes are available:
463	!! - distance-weighted (great circle) (k2dint = 0)
464	!! - distance-weighted (small angle) (k2dint = 1)
465	!! - bilinear (geographical grid) (k2dint = 2)
466	!! - bilinear (quadrilateral grid) (k2dint = 3)
467	!! - polynomial (quadrilateral grid) (k2dint = 4)
468	!!
469	!! Two horizontal averaging schemes are also available:
470	!! - weighted radial footprint (k2dint = 5)
471	!! - weighted rectangular footprint (k2dint = 6)
472	!!
473	!!
474	!! ** Action :
475	!!
476	!! History :
477	!! ! 07-03 (A. Weaver)
478	!! ! 15-02 (M. Martin) Combined routine for surface types
479	!! ! 17-03 (M. Martin) Added horizontal averaging options
480	!!-----------------------------------------------------------------------
481	USE obs_surf_def ! Definition of storage space for surface observations
482
483	IMPLICIT NONE
484
485	TYPE(obs_surf), INTENT(INOUT) :: &
486	& surfdataqc ! Subset of surface data passing QC
487	INTEGER, INTENT(IN) :: kt ! Time step
488	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
489	INTEGER, INTENT(IN) :: kpj
490	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
491	! (kit000-1 = restart time)
492	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
493	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
494	REAL(dp), INTENT(IN), DIMENSION(kpi,kpj) :: &
495	& psurf, & ! Model surface field
496	& psurfmask ! Land-sea mask
497	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
498	REAL(KIND=dp), INTENT(IN) :: &
499	& plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions
500	& pphiscl ! This is the full width (rather than half-width)
501	LOGICAL, INTENT(IN) :: &
502	& lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
503
504	!! * Local declarations
505	INTEGER :: ji
506	INTEGER :: jj
507	INTEGER :: jobs
508	INTEGER :: inrc
509	INTEGER :: isurf
510	INTEGER :: iobs
511	INTEGER :: imaxifp, imaxjfp
512	INTEGER :: imodi, imodj
513	INTEGER :: idayend
514	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
515	& igrdi, &
516	& igrdj, &
517	& igrdip1, &
518	& igrdjp1
519	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
520	& icount_night, &
521	& imask_night
522	REAL(dp) :: zlam
523	REAL(dp) :: zphi
524	REAL(dp), DIMENSION(1) :: zext, zobsmask
525	REAL(dp) :: zdaystp
526	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: &
527	& zweig, &
528	& zmask, &
529	& zsurf, &
530	& zsurfm, &
531	& zsurftmp, &
532	& zglam, &
533	& zgphi, &
534	& zglamf, &
535	& zgphif
536
537	REAL(dp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
538	& zintmp, &
539	& zouttmp, &
540	& zmeanday ! to compute model sst in region of 24h daylight (pole)
541
542	!------------------------------------------------------------------------
543	! Local initialization
544	!------------------------------------------------------------------------
545	! Record and data counters
546	inrc = kt - kit000 + 2
547	isurf = surfdataqc%nsstp(inrc)
548
549	! Work out the maximum footprint size for the
550	! interpolation/averaging in model grid-points - has to be even.
551
552	CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
553
554
555	IF ( ldnightav ) THEN
556
557	! Initialize array for night mean
558	IF ( kt == 0 ) THEN
559	ALLOCATE ( icount_night(kpi,kpj) )
560	ALLOCATE ( imask_night(kpi,kpj) )
561	ALLOCATE ( zintmp(kpi,kpj) )
562	ALLOCATE ( zouttmp(kpi,kpj) )
563	ALLOCATE ( zmeanday(kpi,kpj) )
564	nday_qsr = -1 ! initialisation flag for nbc_dcy
565	ENDIF
566
567	! Night-time means are calculated for night-time values over timesteps:
568	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
569	idayend = MOD( kt - kit000 + 1, kdaystp )
570
571	! Initialize night-time mean for first timestep of the day
572	IF ( idayend == 1 .OR. kt == 0 ) THEN
573	DO_2D( 1, 1, 1, 1 )
574	surfdataqc%vdmean(ji,jj) = 0.0
575	zmeanday(ji,jj) = 0.0
576	icount_night(ji,jj) = 0
577	END_2D
578	ENDIF
579
580	zintmp(:,:) = 0.0
581	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
582	imask_night(:,:) = INT( zouttmp(:,:) )
583
584	DO_2D( 1, 1, 1, 1 )
585	! Increment the temperature field for computing night mean and counter
586	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
587	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
588	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
589	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
590	END_2D
591
592	! Compute the night-time mean at the end of the day
593	zdaystp = 1.0 / REAL( kdaystp )
594	IF ( idayend == 0 ) THEN
595	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
596	DO_2D( 1, 1, 1, 1 )
597	! Test if "no night" point
598	IF ( icount_night(ji,jj) > 0 ) THEN
599	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
600	& / REAL( icount_night(ji,jj) )
601	ELSE
602	!At locations where there is no night (e.g. poles),
603	! calculate daily mean instead of night-time mean.
604	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
605	ENDIF
606	END_2D
607	ENDIF
608
609	ENDIF
610
611	! Get the data for interpolation
612
613	ALLOCATE( &
614	& zweig(imaxifp,imaxjfp,1), &
615	& igrdi(imaxifp,imaxjfp,isurf), &
616	& igrdj(imaxifp,imaxjfp,isurf), &
617	& zglam(imaxifp,imaxjfp,isurf), &
618	& zgphi(imaxifp,imaxjfp,isurf), &
619	& zmask(imaxifp,imaxjfp,isurf), &
620	& zsurf(imaxifp,imaxjfp,isurf), &
621	& zsurftmp(imaxifp,imaxjfp,isurf), &
622	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
623	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
624	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
625	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
626	& )
627
628	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
629	iobs = jobs - surfdataqc%nsurfup
630	DO ji = 0, imaxifp
631	imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
632	!
633	!Deal with wrap around in longitude
634	IF ( imodi < 1 ) imodi = imodi + jpiglo
635	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
636	!
637	DO jj = 0, imaxjfp
638	imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
639	!If model values are out of the domain to the north/south then
640	!set them to be the edge of the domain
641	IF ( imodj < 1 ) imodj = 1
642	IF ( imodj > jpjglo ) imodj = jpjglo
643	!
644	igrdip1(ji+1,jj+1,iobs) = imodi
645	igrdjp1(ji+1,jj+1,iobs) = imodj
646	!
647	IF ( ji >= 1 .AND. jj >= 1 ) THEN
648	igrdi(ji,jj,iobs) = imodi
649	igrdj(ji,jj,iobs) = imodj
650	ENDIF
651	!
652	END DO
653	END DO
654	END DO
655
656	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
657	& igrdi, igrdj, glamt, zglam )
658	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
659	& igrdi, igrdj, gphit, zgphi )
660	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
661	& igrdi, igrdj, psurfmask, zmask )
662	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
663	& igrdi, igrdj, psurf, zsurf )
664	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
665	& igrdip1, igrdjp1, glamf, zglamf )
666	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
667	& igrdip1, igrdjp1, gphif, zgphif )
668
669	! At the end of the day get interpolated means
670	IF ( idayend == 0 .AND. ldnightav ) THEN
671
672	ALLOCATE( &
673	& zsurfm(imaxifp,imaxjfp,isurf) &
674	& )
675
676	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
677	& surfdataqc%vdmean(:,:), zsurfm )
678
679	ENDIF
680
681	! Loop over observations
682	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
683
684	iobs = jobs - surfdataqc%nsurfup
685
686	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
687
688	IF(lwp) THEN
689	WRITE(numout,*)
690	WRITE(numout,*) ' E R R O R : Observation', &
691	& ' time step is not consistent with the', &
692	& ' model time step'
693	WRITE(numout,*) ' ========='
694	WRITE(numout,*)
695	WRITE(numout,*) ' Record = ', jobs, &
696	& ' kt = ', kt, &
697	& ' mstp = ', surfdataqc%mstp(jobs), &
698	& ' ntyp = ', surfdataqc%ntyp(jobs)
699	ENDIF
700	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
701
702	ENDIF
703
704	zlam = surfdataqc%rlam(jobs)
705	zphi = surfdataqc%rphi(jobs)
706
707	IF ( ldnightav .AND. idayend == 0 ) THEN
708	! Night-time averaged data
709	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
710	ELSE
711	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
712	ENDIF
713
714	IF ( k2dint <= 4 ) THEN
715
716	! Get weights to interpolate the model value to the observation point
717	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
718	& zglam(:,:,iobs), zgphi(:,:,iobs), &
719	& zmask(:,:,iobs), zweig, zobsmask )
720
721	! Interpolate the model value to the observation point
722	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
723
724	ELSE
725
726	! Get weights to average the model SLA to the observation footprint
727	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
728	& zglam(:,:,iobs), zgphi(:,:,iobs), &
729	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
730	& zmask(:,:,iobs), plamscl, pphiscl, &
731	& lindegrees, zweig )
732
733	! Average the model SST to the observation footprint
734	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
735	& zweig, zsurftmp(:,:,iobs), zext )
736
737	ENDIF
738
739	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
740	! ... Remove the MDT from the SSH at the observation point to get the SLA
741	surfdataqc%rext(jobs,1) = zext(1)
742	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
743	ELSE
744	surfdataqc%rmod(jobs,1) = zext(1)
745	ENDIF
746
747	IF ( zext(1) == obfillflt ) THEN
748	! If the observation value is a fill value, set QC flag to bad
749	surfdataqc%nqc(jobs) = 4
750	ENDIF
751
752	END DO
753
754	! Deallocate the data for interpolation
755	DEALLOCATE( &
756	& zweig, &
757	& igrdi, &
758	& igrdj, &
759	& zglam, &
760	& zgphi, &
761	& zmask, &
762	& zsurf, &
763	& zsurftmp, &
764	& zglamf, &
765	& zgphif, &
766	& igrdip1,&
767	& igrdjp1 &
768	& )
769
770	! At the end of the day also deallocate night-time mean array
771	IF ( idayend == 0 .AND. ldnightav ) THEN
772	DEALLOCATE( &
773	& zsurfm &
774	& )
775	ENDIF
776	!
777	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
778	!
779	END SUBROUTINE obs_surf_opt
780
781	!!======================================================================
782	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/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obs_oper.F90 @ 15540

Download in other formats: