Context Navigation

obs_write.F90 @ 14165

Last change on this file since 14165 was 14165, checked in by dcarneir, 4 years ago
Merging trunk into my branch to keep it updated
File size: 31.8 KB

Line
1	MODULE obs_write
2	!!======================================================================
3	!! * MODULE obs_write *
4	!! Observation diagnosticss: Write observation related diagnostics
5	!!=====================================================================
6
7	!!----------------------------------------------------------------------
8	!! obs_wri_prof : Write profile observations in feedback format
9	!! obs_wri_surf : Write surface observations in feedback format
10	!! obs_wri_stats : Print basic statistics on the data being written out
11	!!----------------------------------------------------------------------
12
13	!! * Modules used
14	USE par_kind, ONLY : & ! Precision variables
15	& wp
16	USE in_out_manager ! I/O manager
17	USE dom_oce ! Ocean space and time domain variables
18	USE obs_types ! Observation type integer to character translation
19	USE julian, ONLY : & ! Julian date routines
20	& greg2jul
21	USE obs_utils, ONLY : & ! Observation operator utility functions
22	& chkerr
23	USE obs_profiles_def ! Type definitions for profiles
24	USE obs_surf_def ! Type defintions for surface observations
25	USE obs_fbm ! Observation feedback I/O
26	USE obs_grid ! Grid tools
27	USE obs_conv ! Conversion between units
28	USE obs_const
29	USE obs_mpp ! MPP support routines for observation diagnostics
30	USE lib_mpp ! MPP routines
31
32	IMPLICIT NONE
33
34	!! * Routine accessibility
35	PRIVATE
36	PUBLIC obs_wri_prof, & ! Write profile observation files
37	& obs_wri_surf, & ! Write surface observation files
38	& obswriinfo
39
40	TYPE obswriinfo
41	INTEGER :: inum
42	INTEGER, POINTER, DIMENSION(:) :: ipoint
43	CHARACTER(len=ilenname), POINTER, DIMENSION(:) :: cdname
44	CHARACTER(len=ilenlong), POINTER, DIMENSION(:,:) :: cdlong
45	CHARACTER(len=ilenunit), POINTER, DIMENSION(:,:) :: cdunit
46	END TYPE obswriinfo
47
48	!!----------------------------------------------------------------------
49	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
50	!! $Id$
51	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
52	!!----------------------------------------------------------------------
53
54	CONTAINS
55
56	SUBROUTINE obs_wri_prof( profdata, padd, pext )
57	!!-----------------------------------------------------------------------
58	!!
59	!! * ROUTINE obs_wri_prof *
60	!!
61	!! ** Purpose : Write profile feedback files
62	!!
63	!! ** Method : NetCDF
64	!!
65	!! ** Action :
66	!!
67	!! History :
68	!! ! 06-04 (A. Vidard) Original
69	!! ! 06-04 (A. Vidard) Reformatted
70	!! ! 06-10 (A. Weaver) Cleanup
71	!! ! 07-01 (K. Mogensen) Use profile data types
72	!! ! 07-03 (K. Mogensen) General handling of profiles
73	!! ! 09-01 (K. Mogensen) New feedback format
74	!! ! 15-02 (M. Martin) Combined routine for writing profiles
75	!!-----------------------------------------------------------------------
76
77	!! * Arguments
78	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
79	TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable
80	TYPE(obswriinfo), OPTIONAL :: pext ! Extra info
81
82	!! * Local declarations
83	TYPE(obfbdata) :: fbdata
84	CHARACTER(LEN=40) :: clfname
85	CHARACTER(LEN=10) :: clfiletype
86	CHARACTER(LEN=ilenlong) :: cllongname ! Long name of variable
87	CHARACTER(LEN=ilenunit) :: clunits ! Units of variable
88	CHARACTER(LEN=ilengrid) :: clgrid ! Grid of variable
89	INTEGER :: ilevel
90	INTEGER :: jvar
91	INTEGER :: jo
92	INTEGER :: jk
93	INTEGER :: ik
94	INTEGER :: ja
95	INTEGER :: je
96	INTEGER :: iadd
97	INTEGER :: iadd_clm ! 1 if climatology present
98	INTEGER :: iext
99	REAL(wp) :: zpres
100
101
102	iadd_clm = 0
103	IF ( profdata%lclim ) iadd_clm = 1
104
105	IF ( PRESENT( padd ) ) THEN
106	iadd = padd%inum
107	ELSE
108	iadd = 0
109	ENDIF
110
111	IF ( PRESENT( pext ) ) THEN
112	iext = pext%inum
113	ELSE
114	iext = 0
115	ENDIF
116
117	CALL init_obfbdata( fbdata )
118
119	! Find maximum level
120	ilevel = 0
121	DO jvar = 1, profdata%nvar
122	ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
123	END DO
124
125	SELECT CASE ( TRIM(profdata%cvars(1)) )
126	CASE('POTM')
127
128	clfiletype='profb'
129	CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
130	& 1 + iadd_clm + iadd, 1 + iext, .TRUE. )
131	fbdata%cname(1) = profdata%cvars(1)
132	fbdata%cname(2) = profdata%cvars(2)
133	fbdata%coblong(1) = 'Potential temperature'
134	fbdata%coblong(2) = 'Practical salinity'
135	fbdata%cobunit(1) = 'Degrees centigrade'
136	fbdata%cobunit(2) = 'PSU'
137	fbdata%cextname(1) = 'TEMP'
138	fbdata%cextlong(1) = 'Insitu temperature'
139	fbdata%cextunit(1) = 'Degrees centigrade'
140	fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
141	fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
142	fbdata%caddunit(1,1) = 'Degrees centigrade'
143	fbdata%caddunit(1,2) = 'PSU'
144	IF ( profdata%lclim ) THEN
145	fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature'
146	fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity'
147	fbdata%caddunit(2,1) = 'Degrees centigrade'
148	fbdata%caddunit(2,2) = 'PSU'
149	ENDIF
150	fbdata%cgrid(:) = 'T'
151	DO je = 1, iext
152	fbdata%cextname(1+je) = pext%cdname(je)
153	fbdata%cextlong(1+je) = pext%cdlong(je,1)
154	fbdata%cextunit(1+je) = pext%cdunit(je,1)
155	END DO
156	DO ja = 1, iadd
157	fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
158	DO jvar = 1, 2
159	fbdata%caddlong(1+iadd_clm+ja,jvar) = padd%cdlong(ja,jvar)
160	fbdata%caddunit(1+iadd_clm+ja,jvar) = padd%cdunit(ja,jvar)
161	END DO
162	END DO
163
164	CASE('UVEL')
165
166	clfiletype='velfb'
167	CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
168	& 1 + iadd_clm + iadd, 0, .TRUE. )
169	fbdata%cname(1) = profdata%cvars(1)
170	fbdata%cname(2) = profdata%cvars(2)
171	fbdata%coblong(1) = 'Zonal velocity'
172	fbdata%coblong(2) = 'Meridional velocity'
173	fbdata%cobunit(1) = 'm/s'
174	fbdata%cobunit(2) = 'm/s'
175	DO je = 1, iext
176	fbdata%cextname(je) = pext%cdname(je)
177	fbdata%cextlong(je) = pext%cdlong(je,1)
178	fbdata%cextunit(je) = pext%cdunit(je,1)
179	END DO
180	fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
181	fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
182	fbdata%caddunit(1,1) = 'm/s'
183	fbdata%caddunit(1,2) = 'm/s'
184	IF ( profdata%lclim ) THEN
185	fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity'
186	fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity'
187	fbdata%caddunit(2,1) = 'm/s'
188	fbdata%caddunit(2,2) = 'm/s'
189	ENDIF
190	fbdata%cgrid(1) = 'U'
191	fbdata%cgrid(2) = 'V'
192	DO ja = 1, iadd
193	fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
194	fbdata%caddlong(1+iadd_clm+ja,1) = padd%cdlong(ja,1)
195	fbdata%caddunit(1+iadd_clm+ja,1) = padd%cdunit(ja,1)
196	END DO
197
198	CASE('PLCHLTOT')
199
200	clfiletype = 'plchltotfb'
201	cllongname = 'log10(chlorophyll concentration)'
202	clunits = 'log10(mg/m3)'
203	clgrid = 'T'
204
205	CASE('PCHLTOT')
206
207	clfiletype = 'pchltotfb'
208	cllongname = 'chlorophyll concentration'
209	clunits = 'mg/m3'
210	clgrid = 'T'
211
212	CASE('PNO3')
213
214	clfiletype = 'pno3fb'
215	cllongname = 'nitrate'
216	clunits = 'mmol/m3'
217	clgrid = 'T'
218
219	CASE('PSI4')
220
221	clfiletype = 'psi4fb'
222	cllongname = 'silicate'
223	clunits = 'mmol/m3'
224	clgrid = 'T'
225
226	CASE('PPO4')
227
228	clfiletype = 'ppo4fb'
229	cllongname = 'phosphate'
230	clunits = 'mmol/m3'
231	clgrid = 'T'
232
233	CASE('PDIC')
234
235	clfiletype = 'pdicfb'
236	cllongname = 'dissolved inorganic carbon'
237	clunits = 'mmol/m3'
238	clgrid = 'T'
239
240	CASE('PALK')
241
242	clfiletype = 'palkfb'
243	cllongname = 'alkalinity'
244	clunits = 'meq/m3'
245	clgrid = 'T'
246
247	CASE('PPH')
248
249	clfiletype = 'pphfb'
250	cllongname = 'pH'
251	clunits = '-'
252	clgrid = 'T'
253
254	CASE('PO2')
255
256	clfiletype = 'po2fb'
257	cllongname = 'dissolved oxygen'
258	clunits = 'mmol/m3'
259	clgrid = 'T'
260
261	END SELECT
262
263	IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. &
264	& ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
265	CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
266	& 1 + iadd_clm + iadd, iext, .TRUE. )
267	fbdata%cname(1) = profdata%cvars(1)
268	fbdata%coblong(1) = cllongname
269	fbdata%cobunit(1) = clunits
270	fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
271	fbdata%caddunit(1,1) = clunits
272	IF ( profdata%lclim ) THEN
273	fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname)
274	fbdata%caddunit(2,1) = clunits
275	ENDIF
276	fbdata%cgrid(:) = clgrid
277	DO je = 1, iext
278	fbdata%cextname(je) = pext%cdname(je)
279	fbdata%cextlong(je) = pext%cdlong(je,1)
280	fbdata%cextunit(je) = pext%cdunit(je,1)
281	END DO
282	DO ja = 1, iadd
283	fbdata%caddname(1+iadd_clm+ja) = padd%cdname(ja)
284	fbdata%caddlong(1+iadd_clm+ja,1) = padd%cdlong(ja,1)
285	fbdata%caddunit(1+iadd_clm+ja,1) = padd%cdunit(ja,1)
286	END DO
287	ENDIF
288
289	fbdata%caddname(1) = 'Hx'
290	IF ( profdata%lclim ) fbdata%caddname(1+iadd_clm) = 'CLM'
291
292	WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
293
294	IF(lwp) THEN
295	WRITE(numout,*)
296	WRITE(numout,*)'obs_wri_prof :'
297	WRITE(numout,*)'~~~~~~~~~~~~~'
298	WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
299	ENDIF
300
301	! Transform obs_prof data structure into obfb data structure
302	fbdata%cdjuldref = '19500101000000'
303	DO jo = 1, profdata%nprof
304	fbdata%plam(jo) = profdata%rlam(jo)
305	fbdata%pphi(jo) = profdata%rphi(jo)
306	WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
307	fbdata%ivqc(jo,:) = profdata%ivqc(jo,:)
308	fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
309	IF ( profdata%nqc(jo) > 255 ) THEN
310	fbdata%ioqc(jo) = IBSET(profdata%nqc(jo),2)
311	fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
312	fbdata%ioqcf(2,jo) = profdata%nqc(jo)
313	ELSE
314	fbdata%ioqc(jo) = profdata%nqc(jo)
315	fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
316	ENDIF
317	fbdata%ipqc(jo) = profdata%ipqc(jo)
318	fbdata%ipqcf(:,jo) = profdata%ipqcf(:,jo)
319	fbdata%itqc(jo) = profdata%itqc(jo)
320	fbdata%itqcf(:,jo) = profdata%itqcf(:,jo)
321	fbdata%cdwmo(jo) = profdata%cwmo(jo)
322	fbdata%kindex(jo) = profdata%npfil(jo)
323	DO jvar = 1, profdata%nvar
324	IF (ln_grid_global) THEN
325	fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
326	fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
327	ELSE
328	fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
329	fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
330	ENDIF
331	END DO
332	CALL greg2jul( 0, &
333	& profdata%nmin(jo), &
334	& profdata%nhou(jo), &
335	& profdata%nday(jo), &
336	& profdata%nmon(jo), &
337	& profdata%nyea(jo), &
338	& fbdata%ptim(jo), &
339	& krefdate = 19500101 )
340	! Reform the profiles arrays for output
341	DO jvar = 1, profdata%nvar
342	DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
343	ik = profdata%var(jvar)%nvlidx(jk)
344	fbdata%pob(ik,jo,jvar) = profdata%var(jvar)%vobs(jk)
345	fbdata%pdep(ik,jo) = profdata%var(jvar)%vdep(jk)
346	fbdata%idqc(ik,jo) = profdata%var(jvar)%idqc(jk)
347	fbdata%idqcf(:,ik,jo) = profdata%var(jvar)%idqcf(:,jk)
348	IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
349	fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
350	fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
351	fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111')
352	ELSE
353	fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
354	fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
355	ENDIF
356	fbdata%iobsk(ik,jo,jvar) = profdata%var(jvar)%mvk(jk)
357
358	fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
359	IF ( profdata%lclim ) THEN
360	fbdata%padd(ik,jo,1+iadd_clm,jvar) = profdata%var(jvar)%vclm(jk)
361	ENDIF
362	DO ja = 1, iadd
363	fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = &
364	& profdata%var(jvar)%vext(jk,padd%ipoint(ja))
365	END DO
366	DO je = 1, iext
367	fbdata%pext(ik,jo,1+je) = &
368	& profdata%var(jvar)%vext(jk,pext%ipoint(je))
369	END DO
370	IF ( ( jvar == 1 ) .AND. &
371	& ( TRIM(profdata%cvars(1)) == 'POTM' ) ) THEN
372	fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
373	ENDIF
374	END DO
375	END DO
376	END DO
377
378	IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
379	! Convert insitu temperature to potential temperature using the model
380	! salinity if no potential temperature
381	DO jo = 1, fbdata%nobs
382	IF ( fbdata%pphi(jo) < 9999.0 ) THEN
383	DO jk = 1, fbdata%nlev
384	IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
385	& ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
386	& ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
387	& ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
388	zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
389	& REAL(fbdata%pphi(jo),wp) )
390	fbdata%pob(jk,jo,1) = potemp( &
391	& REAL(fbdata%padd(jk,jo,1,2), wp), &
392	& REAL(fbdata%pext(jk,jo,1), wp), &
393	& zpres, 0.0_wp )
394	ENDIF
395	END DO
396	ENDIF
397	END DO
398	ENDIF
399
400	! Write the obfbdata structure
401	CALL write_obfbdata( clfname, fbdata )
402
403	! Output some basic statistics
404	CALL obs_wri_stats( fbdata )
405
406	CALL dealloc_obfbdata( fbdata )
407
408	END SUBROUTINE obs_wri_prof
409
410	SUBROUTINE obs_wri_surf( surfdata, padd, pext )
411	!!-----------------------------------------------------------------------
412	!!
413	!! * ROUTINE obs_wri_surf *
414	!!
415	!! ** Purpose : Write surface observation files
416	!!
417	!! ** Method : NetCDF
418	!!
419	!! ** Action :
420	!!
421	!! ! 07-03 (K. Mogensen) Original
422	!! ! 09-01 (K. Mogensen) New feedback format.
423	!! ! 15-02 (M. Martin) Combined surface writing routine.
424	!!-----------------------------------------------------------------------
425
426	!! * Modules used
427	IMPLICIT NONE
428
429	!! * Arguments
430	TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data
431	TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable
432	TYPE(obswriinfo), OPTIONAL :: pext ! Extra info
433
434	!! * Local declarations
435	TYPE(obfbdata) :: fbdata
436	CHARACTER(LEN=40) :: clfname ! netCDF filename
437	CHARACTER(LEN=10) :: clfiletype
438	CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
439	CHARACTER(LEN=ilenlong), DIMENSION(surfdata%nvar) :: cllongname ! Long name of variable
440	CHARACTER(LEN=ilenunit), DIMENSION(surfdata%nvar) :: clunits ! Units of variable
441	CHARACTER(LEN=ilengrid), DIMENSION(surfdata%nvar) :: clgrid ! Grid of variable
442	INTEGER :: jo
443	INTEGER :: ja
444	INTEGER :: je
445	INTEGER :: jv
446	INTEGER :: iadd
447	INTEGER :: iext
448	INTEGER :: indx_std
449	INTEGER :: iadd_std
450	INTEGER :: iadd_clm
451	INTEGER :: iadd_mdt
452
453	IF ( PRESENT( pext ) ) THEN
454	iext = pext%inum
455	ELSE
456	iext = 0
457	ENDIF
458
459
460	! Set up number of additional variables to be ouput:
461	! Hx, CLM, STD, MDT...
462
463	IF ( PRESENT( padd ) ) THEN
464	iadd = padd%inum
465	ELSE
466	iadd = 0
467	ENDIF
468
469	iadd_std = 0
470	indx_std = -1
471	IF ( surfdata%nextra > 0 ) THEN
472	DO je = 1, surfdata%nextra
473	IF ( TRIM( surfdata%cext(je) ) == 'STD' ) THEN
474	iadd_std = 1
475	indx_std = je
476	ENDIF
477	END DO
478	ENDIF
479
480	iadd_clm = 0
481	IF ( surfdata%lclim ) iadd_clm = 1
482
483	iadd_mdt = 0
484	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1
485
486	CALL init_obfbdata( fbdata )
487
488	SELECT CASE ( TRIM(surfdata%cvars(1)) )
489	CASE('SLA')
490
491	! SLA needs special treatment because of MDT, so is all done here
492	! Other variables are done more generically
493	! No climatology for SLA, MDT is our best estimate of that and is already output.
494
495	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
496	& 1 + iadd_mdt + iadd_std + iadd, &
497	& 1 + iext, .TRUE. )
498
499	clfiletype = 'slafb'
500	fbdata%cname(1) = surfdata%cvars(1)
501	fbdata%coblong(1) = 'Sea level anomaly'
502	fbdata%cobunit(1) = 'Metres'
503	fbdata%cextname(1) = 'MDT'
504	fbdata%cextlong(1) = 'Mean dynamic topography'
505	fbdata%cextunit(1) = 'Metres'
506	DO je = 1, iext
507	fbdata%cextname(je) = pext%cdname(je)
508	fbdata%cextlong(je) = pext%cdlong(je,1)
509	fbdata%cextunit(je) = pext%cdunit(je,1)
510	END DO
511	fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
512	fbdata%caddunit(1,1) = 'Metres'
513	fbdata%caddname(2) = 'SSH'
514	fbdata%caddlong(2,1) = 'Model Sea surface height'
515	fbdata%caddunit(2,1) = 'Metres'
516	fbdata%cgrid(1) = 'T'
517	DO ja = 1, iadd
518	fbdata%caddname(1+iadd_mdt+iadd_std+ja) = padd%cdname(ja)
519	fbdata%caddlong(1+iadd_mdt+iadd_std+ja,1) = padd%cdlong(ja,1)
520	fbdata%caddunit(1+iadd_mdt+iadd_std+ja,1) = padd%cdunit(ja,1)
521	END DO
522
523	CASE('SST')
524
525	clfiletype = 'sstfb'
526	cllongname(1) = 'Sea surface temperature'
527	clunits(1) = 'Degree centigrade'
528	clgrid(1) = 'T'
529
530	CASE('ICECONC')
531
532	clfiletype = 'sicfb'
533	cllongname(1) = 'Sea ice concentration'
534	clunits(1) = 'Fraction'
535	clgrid(1) = 'T'
536
537	CASE('SIT')
538
539	clfiletype = 'sitfb'
540	cllongname(1) = 'Sea ice thickness'
541	clunits(1) = 'm'
542	clgrid(1) = 'T'
543
544	CASE('FBD')
545
546	clfiletype = 'fbdfb'
547	! Change label from FBD ("freeboard") to SIT ("Sea Ice Thickness")
548	surfdata%cvars(1) = 'SIT'
549	cllongname = 'Sea ice thickness'
550	clunits = 'm'
551	clgrid = 'T'
552
553	CASE('SSS')
554
555	clfiletype = 'sssfb'
556	cllongname(1) = 'Sea surface salinity'
557	clunits(1) = 'psu'
558	clgrid(1) = 'T'
559
560	CASE('UVEL')
561
562	clfiletype = 'ssvfb'
563	cllongname(1) = 'Eastward sea surface velocity'
564	clunits(1) = 'm s-1'
565	clgrid(1) = 'U'
566	cllongname(2) = 'Northward sea surface velocity'
567	clunits(2) = 'm s-1'
568	clgrid(2) = 'V'
569
570	CASE('SLCHLTOT')
571
572	clfiletype = 'slchltotfb'
573	cllongname(1) = 'Surface total log10(chlorophyll)'
574	clunits(1) = 'log10(mg/m3)'
575	clgrid(1) = 'T'
576
577	CASE('SLCHLDIA')
578
579	clfiletype = 'slchldiafb'
580	cllongname(1) = 'Surface diatom log10(chlorophyll)'
581	clunits(1) = 'log10(mg/m3)'
582	clgrid(1) = 'T'
583
584	CASE('SLCHLNON')
585
586	clfiletype = 'slchlnonfb'
587	cllongname(1) = 'Surface non-diatom log10(chlorophyll)'
588	clunits(1) = 'log10(mg/m3)'
589	clgrid(1) = 'T'
590
591	CASE('SLCHLDIN')
592
593	clfiletype = 'slchldinfb'
594	cllongname(1) = 'Surface dinoflagellate log10(chlorophyll)'
595	clunits(1) = 'log10(mg/m3)'
596	clgrid(1) = 'T'
597
598	CASE('SLCHLMIC')
599
600	clfiletype = 'slchlmicfb'
601	cllongname(1) = 'Surface microphytoplankton log10(chlorophyll)'
602	clunits(1) = 'log10(mg/m3)'
603	clgrid(1) = 'T'
604
605	CASE('SLCHLNAN')
606
607	clfiletype = 'slchlnanfb'
608	cllongname(1) = 'Surface nanophytoplankton log10(chlorophyll)'
609	clunits(1) = 'log10(mg/m3)'
610	clgrid(1) = 'T'
611
612	CASE('SLCHLPIC')
613
614	clfiletype = 'slchlpicfb'
615	cllongname(1) = 'Surface picophytoplankton log10(chlorophyll)'
616	clunits(1) = 'log10(mg/m3)'
617	clgrid(1) = 'T'
618
619	CASE('SCHLTOT')
620
621	clfiletype = 'schltotfb'
622	cllongname = 'Surface total chlorophyll'
623	clunits = 'mg/m3'
624	clgrid = 'T'
625
626	CASE('SLPHYTOT')
627
628	clfiletype = 'slphytotfb'
629	cllongname(1) = 'Surface total log10(phytoplankton carbon)'
630	clunits(1) = 'log10(mmolC/m3)'
631	clgrid(1) = 'T'
632
633	CASE('SLPHYDIA')
634
635	clfiletype = 'slphydiafb'
636	cllongname(1) = 'Surface diatom log10(phytoplankton carbon)'
637	clunits(1) = 'log10(mmolC/m3)'
638	clgrid(1) = 'T'
639
640	CASE('SLPHYNON')
641
642	clfiletype = 'slphynonfb'
643	cllongname(1) = 'Surface non-diatom log10(phytoplankton carbon)'
644	clunits(1) = 'log10(mmolC/m3)'
645	clgrid(1) = 'T'
646
647	CASE('SSPM')
648
649	clfiletype = 'sspmfb'
650	cllongname(1) = 'Surface suspended particulate matter'
651	clunits(1) = 'g/m3'
652	clgrid(1) = 'T'
653
654	CASE('SKD490')
655
656	clfiletype = 'skd490fb'
657	cllongname(1) = 'Surface attenuation coefficient of downwelling radiation at 490 nm'
658	clunits(1) = 'm-1'
659	clgrid(1) = 'T'
660
661	CASE('SFCO2')
662
663	clfiletype = 'sfco2fb'
664	cllongname(1) = 'Surface fugacity of carbon dioxide'
665	clunits(1) = 'uatm'
666	clgrid(1) = 'T'
667
668	CASE('SPCO2')
669
670	clfiletype = 'spco2fb'
671	cllongname(1) = 'Surface partial pressure of carbon dioxide'
672	clunits(1) = 'uatm'
673	clgrid(1) = 'T'
674
675	CASE DEFAULT
676
677	CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
678
679	END SELECT
680
681	! SLA needs special treatment because of MDT, so is done above
682	! Remaining variables treated more generically
683
684	IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
685
686	CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, &
687	& 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
688
689	DO jv = 1, surfdata%nvar
690	fbdata%cname(jv) = surfdata%cvars(jv)
691	fbdata%coblong(jv) = cllongname(jv)
692	fbdata%cobunit(jv) = clunits(jv)
693	END DO
694	DO je = 1, iext
695	fbdata%cextname(je) = pext%cdname(je)
696	fbdata%cextlong(je) = pext%cdlong(je,1)
697	fbdata%cextunit(je) = pext%cdunit(je,1)
698	END DO
699	DO jv = 1, surfdata%nvar
700	IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
701	fbdata%caddlong(1,jv) = 'Model interpolated ICE'
702	ELSE
703	fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv))
704	ENDIF
705	fbdata%caddunit(1,jv) = clunits(jv)
706	fbdata%cgrid(jv) = clgrid(jv)
707	END DO
708	DO ja = 1, iadd
709	fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)
710	DO jv = 1, surfdata%nvar
711	fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv)
712	fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv)
713	END DO
714	END DO
715	ENDIF
716
717	fbdata%caddname(1) = 'Hx'
718	IF ( indx_std /= -1 ) THEN
719	fbdata%caddname(1+iadd_mdt+iadd_std) = surfdata%cext(indx_std)
720	DO jv = 1, surfdata%nvar
721	fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'
722	fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)
723	END DO
724	ENDIF
725
726	IF ( surfdata%lclim ) THEN
727	fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm) = 'CLM'
728	DO jv = 1, surfdata%nvar
729	fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology'
730	fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1)
731	END DO
732	ENDIF
733
734	WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
735
736	IF(lwp) THEN
737	WRITE(numout,*)
738	WRITE(numout,*)'obs_wri_surf :'
739	WRITE(numout,*)'~~~~~~~~~~~~~'
740	WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
741	ENDIF
742
743	! Transform surf data structure into obfbdata structure
744	fbdata%cdjuldref = '19500101000000'
745	DO jo = 1, surfdata%nsurf
746	fbdata%plam(jo) = surfdata%rlam(jo)
747	fbdata%pphi(jo) = surfdata%rphi(jo)
748	WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
749	fbdata%ivqc(jo,:) = 0
750	fbdata%ivqcf(:,jo,:) = 0
751	IF ( surfdata%nqc(jo) > 255 ) THEN
752	fbdata%ioqc(jo) = 4
753	fbdata%ioqcf(1,jo) = 0
754	fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111')
755	ELSE
756	fbdata%ioqc(jo) = surfdata%nqc(jo)
757	fbdata%ioqcf(:,jo) = 0
758	ENDIF
759	fbdata%ipqc(jo) = 0
760	fbdata%ipqcf(:,jo) = 0
761	fbdata%itqc(jo) = 0
762	fbdata%itqcf(:,jo) = 0
763	fbdata%cdwmo(jo) = surfdata%cwmo(jo)
764	fbdata%kindex(jo) = surfdata%nsfil(jo)
765	IF (ln_grid_global) THEN
766	DO jv = 1, surfdata%nvar
767	fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv)
768	fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv)
769	END DO
770	ELSE
771	DO jv = 1, surfdata%nvar
772	fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv))
773	fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv))
774	END DO
775	ENDIF
776	CALL greg2jul( 0, &
777	& surfdata%nmin(jo), &
778	& surfdata%nhou(jo), &
779	& surfdata%nday(jo), &
780	& surfdata%nmon(jo), &
781	& surfdata%nyea(jo), &
782	& fbdata%ptim(jo), &
783	& krefdate = 19500101 )
784
785	fbdata%pdep(1,jo) = 0.0
786	fbdata%idqc(1,jo) = 0
787	fbdata%idqcf(:,1,jo) = 0
788	DO jv = 1, surfdata%nvar
789	fbdata%pob(1,jo,jv) = surfdata%robs(jo,jv)
790
791	IF ( surfdata%nqc(jo) > 255 ) THEN
792	fbdata%ivqc(jo,jv) = 4
793	fbdata%ivlqc(1,jo,jv) = 4
794	fbdata%ivlqcf(1,1,jo,jv) = 0
795	fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000000011111111')
796	ELSE
797	fbdata%ivqc(jo,jv) = surfdata%nqc(jo)
798	fbdata%ivlqc(1,jo,jv) = surfdata%nqc(jo)
799	fbdata%ivlqcf(:,1,jo,jv) = 0
800	ENDIF
801	fbdata%iobsk(1,jo,jv) = 0
802
803
804	! Additional variables.
805	! Hx is always the first additional variable
806	fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv)
807	! MDT is output as an additional variable if SLA obs type
808	IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN
809	fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1)
810	ENDIF
811	! STD is output as an additional variable if available
812	IF ( indx_std /= -1 ) THEN
813	fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std)
814	ENDIF
815	! CLM is output as an additional variable if available
816	IF ( surfdata%lclim ) THEN
817	fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv)
818	ENDIF
819	! Then other additional variables are output
820	DO ja = 1, iadd
821	fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = &
822	& surfdata%rext(jo,padd%ipoint(ja))
823	END DO
824	END DO
825	! Extra variables
826	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
827	DO je = 1, iext
828	fbdata%pext(1,jo,1+je) = &
829	& surfdata%rext(jo,pext%ipoint(je))
830	END DO
831	END DO
832
833	! Write the obfbdata structure
834	CALL write_obfbdata( clfname, fbdata )
835
836	! Output some basic statistics
837	CALL obs_wri_stats( fbdata )
838
839	CALL dealloc_obfbdata( fbdata )
840
841	END SUBROUTINE obs_wri_surf
842
843	SUBROUTINE obs_wri_stats( fbdata )
844	!!-----------------------------------------------------------------------
845	!!
846	!! * ROUTINE obs_wri_stats *
847	!!
848	!! ** Purpose : Output some basic statistics of the data being written out
849	!!
850	!! ** Method :
851	!!
852	!! ** Action :
853	!!
854	!! ! 2014-08 (D. Lea) Initial version
855	!!-----------------------------------------------------------------------
856
857	!! * Arguments
858	TYPE(obfbdata) :: fbdata
859
860	!! * Local declarations
861	INTEGER :: jvar
862	INTEGER :: jo
863	INTEGER :: jk
864	INTEGER :: inumgoodobs
865	INTEGER :: inumgoodobsmpp
866	REAL(wp) :: zsumx
867	REAL(wp) :: zsumx2
868	REAL(wp) :: zomb
869
870
871	IF (lwp) THEN
872	WRITE(numout,*) ''
873	WRITE(numout,*) 'obs_wri_stats :'
874	WRITE(numout,*) '~~~~~~~~~~~~~~~'
875	ENDIF
876
877	DO jvar = 1, fbdata%nvar
878	zsumx=0.0_wp
879	zsumx2=0.0_wp
880	inumgoodobs=0
881	DO jo = 1, fbdata%nobs
882	DO jk = 1, fbdata%nlev
883	IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
884	& ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
885	& ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
886
887	zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
888	zsumx=zsumx+zomb
889	zsumx2=zsumx2+zomb**2
890	inumgoodobs=inumgoodobs+1
891	ENDIF
892	ENDDO
893	ENDDO
894
895	CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
896	CALL mpp_sum(zsumx)
897	CALL mpp_sum(zsumx2)
898
899	IF (lwp) THEN
900	WRITE(numout,*) 'Type: ',fbdata%cname(jvar),' Total number of good observations: ',inumgoodobsmpp
901	WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
902	WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
903	WRITE(numout,*) ''
904	ENDIF
905
906	ENDDO
907
908	END SUBROUTINE obs_wri_stats
909
910	END MODULE obs_write

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_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 14165

Download in other formats: