Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obs_write.F90 @ 11932

Last change on this file since 11932 was 11932, checked in by dcarneir, 4 years ago
Changes in OBS and SBC routines for sea ice thickness data assimilation
File size: 30.2 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'0000 0000 1111 1111')
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) :: cllongname ! Long name of variable
440	CHARACTER(LEN=ilenunit) :: clunits ! Units of variable
441	CHARACTER(LEN=ilengrid) :: clgrid ! Grid of variable
442	INTEGER :: jo
443	INTEGER :: ja
444	INTEGER :: je
445	INTEGER :: iadd
446	INTEGER :: iext
447	INTEGER :: indx_std
448	INTEGER :: iadd_std
449	INTEGER :: iadd_clm
450	INTEGER :: iadd_mdt
451
452	IF ( PRESENT( pext ) ) THEN
453	iext = pext%inum
454	ELSE
455	iext = 0
456	ENDIF
457
458
459	! Set up number of additional variables to be ouput:
460	! Hx, CLM, STD, MDT...
461
462	IF ( PRESENT( padd ) ) THEN
463	iadd = padd%inum
464	ELSE
465	iadd = 0
466	ENDIF
467
468	iadd_std = 0
469	indx_std = -1
470	IF ( surfdata%nextra > 0 ) THEN
471	DO je = 1, surfdata%nextra
472	IF ( TRIM( surfdata%cext(je) ) == 'STD' ) THEN
473	iadd_std = 1
474	indx_std = je
475	ENDIF
476	END DO
477	ENDIF
478
479	iadd_clm = 0
480	IF ( surfdata%lclim ) iadd_clm = 1
481
482	iadd_mdt = 0
483	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1
484
485	CALL init_obfbdata( fbdata )
486
487	SELECT CASE ( TRIM(surfdata%cvars(1)) )
488	CASE('SLA')
489
490	! SLA needs special treatment because of MDT, so is all done here
491	! Other variables are done more generically
492	! No climatology for SLA, MDT is our best estimate of that and is already output.
493
494	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
495	& 1 + iadd_mdt + iadd_std + iadd, &
496	& 1 + iext, .TRUE. )
497
498	clfiletype = 'slafb'
499	fbdata%cname(1) = surfdata%cvars(1)
500	fbdata%coblong(1) = 'Sea level anomaly'
501	fbdata%cobunit(1) = 'Metres'
502	fbdata%cextname(1) = 'MDT'
503	fbdata%cextlong(1) = 'Mean dynamic topography'
504	fbdata%cextunit(1) = 'Metres'
505	DO je = 1, iext
506	fbdata%cextname(je) = pext%cdname(je)
507	fbdata%cextlong(je) = pext%cdlong(je,1)
508	fbdata%cextunit(je) = pext%cdunit(je,1)
509	END DO
510	fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
511	fbdata%caddunit(1,1) = 'Metres'
512	fbdata%caddname(2) = 'SSH'
513	fbdata%caddlong(2,1) = 'Model Sea surface height'
514	fbdata%caddunit(2,1) = 'Metres'
515	fbdata%cgrid(1) = 'T'
516	DO ja = 1, iadd
517	fbdata%caddname(1+iadd_mdt+iadd_std+ja) = padd%cdname(ja)
518	fbdata%caddlong(1+iadd_mdt+iadd_std+ja,1) = padd%cdlong(ja,1)
519	fbdata%caddunit(1+iadd_mdt+iadd_std+ja,1) = padd%cdunit(ja,1)
520	END DO
521
522	CASE('SST')
523
524	clfiletype = 'sstfb'
525	cllongname = 'Sea surface temperature'
526	clunits = 'Degree centigrade'
527	clgrid = 'T'
528
529	CASE('ICECONC')
530
531	clfiletype = 'sicfb'
532	cllongname = 'Sea ice concentration'
533	clunits = 'Fraction'
534	clgrid = 'T'
535
536	CASE('SIT')
537
538	clfiletype = 'sitfb'
539	cllongname = 'Sea ice thickness'
540	clunits = 'm'
541	clgrid = 'T'
542
543	CASE('SSS')
544
545	clfiletype = 'sssfb'
546	cllongname = 'Sea surface salinity'
547	clunits = 'psu'
548	clgrid = 'T'
549
550	CASE('SLCHLTOT')
551
552	clfiletype = 'slchltotfb'
553	cllongname = 'Surface total log10(chlorophyll)'
554	clunits = 'log10(mg/m3)'
555	clgrid = 'T'
556
557	CASE('SLCHLDIA')
558
559	clfiletype = 'slchldiafb'
560	cllongname = 'Surface diatom log10(chlorophyll)'
561	clunits = 'log10(mg/m3)'
562	clgrid = 'T'
563
564	CASE('SLCHLNON')
565
566	clfiletype = 'slchlnonfb'
567	cllongname = 'Surface non-diatom log10(chlorophyll)'
568	clunits = 'log10(mg/m3)'
569	clgrid = 'T'
570
571	CASE('SLCHLDIN')
572
573	clfiletype = 'slchldinfb'
574	cllongname = 'Surface dinoflagellate log10(chlorophyll)'
575	clunits = 'log10(mg/m3)'
576	clgrid = 'T'
577
578	CASE('SLCHLMIC')
579
580	clfiletype = 'slchlmicfb'
581	cllongname = 'Surface microphytoplankton log10(chlorophyll)'
582	clunits = 'log10(mg/m3)'
583	clgrid = 'T'
584
585	CASE('SLCHLNAN')
586
587	clfiletype = 'slchlnanfb'
588	cllongname = 'Surface nanophytoplankton log10(chlorophyll)'
589	clunits = 'log10(mg/m3)'
590	clgrid = 'T'
591
592	CASE('SLCHLPIC')
593
594	clfiletype = 'slchlpicfb'
595	cllongname = 'Surface picophytoplankton log10(chlorophyll)'
596	clunits = 'log10(mg/m3)'
597	clgrid = 'T'
598
599	CASE('SCHLTOT')
600
601	clfiletype = 'schltotfb'
602	cllongname = 'Surface total chlorophyll'
603	clunits = 'mg/m3'
604	clgrid = 'T'
605
606	CASE('SLPHYTOT')
607
608	clfiletype = 'slphytotfb'
609	cllongname = 'Surface total log10(phytoplankton carbon)'
610	clunits = 'log10(mmolC/m3)'
611	clgrid = 'T'
612
613	CASE('SLPHYDIA')
614
615	clfiletype = 'slphydiafb'
616	cllongname = 'Surface diatom log10(phytoplankton carbon)'
617	clunits = 'log10(mmolC/m3)'
618	clgrid = 'T'
619
620	CASE('SLPHYNON')
621
622	clfiletype = 'slphynonfb'
623	cllongname = 'Surface non-diatom log10(phytoplankton carbon)'
624	clunits = 'log10(mmolC/m3)'
625	clgrid = 'T'
626
627	CASE('SSPM')
628
629	clfiletype = 'sspmfb'
630	cllongname = 'Surface suspended particulate matter'
631	clunits = 'g/m3'
632	clgrid = 'T'
633
634	CASE('SKD490')
635
636	clfiletype = 'skd490fb'
637	cllongname = 'Surface attenuation coefficient of downwelling radiation at 490 nm'
638	clunits = 'm-1'
639	clgrid = 'T'
640
641	CASE('SFCO2')
642
643	clfiletype = 'sfco2fb'
644	cllongname = 'Surface fugacity of carbon dioxide'
645	clunits = 'uatm'
646	clgrid = 'T'
647
648	CASE('SPCO2')
649
650	clfiletype = 'spco2fb'
651	cllongname = 'Surface partial pressure of carbon dioxide'
652	clunits = 'uatm'
653	clgrid = 'T'
654
655	CASE DEFAULT
656
657	CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
658
659	END SELECT
660
661	! SLA needs special treatment because of MDT, so is done above
662	! Remaining variables treated more generically
663
664	IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
665
666	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
667	& 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
668
669	fbdata%cname(1) = surfdata%cvars(1)
670	fbdata%coblong(1) = cllongname
671	fbdata%cobunit(1) = clunits
672	DO je = 1, iext
673	fbdata%cextname(je) = pext%cdname(je)
674	fbdata%cextlong(je) = pext%cdlong(je,1)
675	fbdata%cextunit(je) = pext%cdunit(je,1)
676	END DO
677	IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
678	fbdata%caddlong(1,1) = 'Model interpolated ICE'
679	ELSE
680	fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1))
681	ENDIF
682	fbdata%caddunit(1,1) = clunits
683	fbdata%cgrid(1) = clgrid
684	DO ja = 1, iadd
685	fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)
686	fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1)
687	fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1)
688	END DO
689
690	ENDIF
691
692	fbdata%caddname(1) = 'Hx'
693	IF ( indx_std /= -1 ) THEN
694	fbdata%caddname(1+iadd_mdt+iadd_std) = surfdata%cext(indx_std)
695	fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'
696	fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)
697	ENDIF
698
699	IF ( surfdata%lclim ) THEN
700	fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm) = 'CLM'
701	fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,1) = 'Climatology'
702	fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,1) = fbdata%cobunit(1)
703	ENDIF
704
705
706	WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
707
708	IF(lwp) THEN
709	WRITE(numout,*)
710	WRITE(numout,*)'obs_wri_surf :'
711	WRITE(numout,*)'~~~~~~~~~~~~~'
712	WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
713	ENDIF
714
715	! Transform surf data structure into obfbdata structure
716	fbdata%cdjuldref = '19500101000000'
717	DO jo = 1, surfdata%nsurf
718	fbdata%plam(jo) = surfdata%rlam(jo)
719	fbdata%pphi(jo) = surfdata%rphi(jo)
720	WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
721	fbdata%ivqc(jo,:) = 0
722	fbdata%ivqcf(:,jo,:) = 0
723	IF ( surfdata%nqc(jo) > 255 ) THEN
724	fbdata%ioqc(jo) = 4
725	fbdata%ioqcf(1,jo) = 0
726	fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
727	ELSE
728	fbdata%ioqc(jo) = surfdata%nqc(jo)
729	fbdata%ioqcf(:,jo) = 0
730	ENDIF
731	fbdata%ipqc(jo) = 0
732	fbdata%ipqcf(:,jo) = 0
733	fbdata%itqc(jo) = 0
734	fbdata%itqcf(:,jo) = 0
735	fbdata%cdwmo(jo) = surfdata%cwmo(jo)
736	fbdata%kindex(jo) = surfdata%nsfil(jo)
737	IF (ln_grid_global) THEN
738	fbdata%iobsi(jo,1) = surfdata%mi(jo)
739	fbdata%iobsj(jo,1) = surfdata%mj(jo)
740	ELSE
741	fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
742	fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
743	ENDIF
744	CALL greg2jul( 0, &
745	& surfdata%nmin(jo), &
746	& surfdata%nhou(jo), &
747	& surfdata%nday(jo), &
748	& surfdata%nmon(jo), &
749	& surfdata%nyea(jo), &
750	& fbdata%ptim(jo), &
751	& krefdate = 19500101 )
752
753	fbdata%pob(1,jo,1) = surfdata%robs(jo,1)
754	fbdata%pdep(1,jo) = 0.0
755	fbdata%idqc(1,jo) = 0
756	fbdata%idqcf(:,1,jo) = 0
757	IF ( surfdata%nqc(jo) > 255 ) THEN
758	fbdata%ivqc(jo,1) = 4
759	fbdata%ivlqc(1,jo,1) = 4
760	fbdata%ivlqcf(1,1,jo,1) = 0
761	fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
762	ELSE
763	fbdata%ivqc(jo,1) = surfdata%nqc(jo)
764	fbdata%ivlqc(1,jo,1) = surfdata%nqc(jo)
765	fbdata%ivlqcf(:,1,jo,1) = 0
766	ENDIF
767	fbdata%iobsk(1,jo,1) = 0
768
769	! Additional variables.
770	! Hx is always the first additional variable
771	fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
772	! MDT is output as an additional variable if SLA obs type
773	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
774	fbdata%padd(1,jo,1+iadd_mdt,1) = surfdata%rext(jo,1)
775	ENDIF
776	! STD is output as an additional variable if available
777	IF ( indx_std /= -1 ) THEN
778	fbdata%padd(1,jo,1+iadd_mdt+iadd_std,1) = surfdata%rext(jo,indx_std)
779	ENDIF
780	! CLM is output as an additional variable if available
781	IF ( surfdata%lclim ) THEN
782	fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,1) = surfdata%rclm(jo,1)
783	ENDIF
784	! Then other additional variables are output
785	DO ja = 1, iadd
786	fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,1) = &
787	& surfdata%rext(jo,padd%ipoint(ja))
788	END DO
789
790	! Extra variables
791	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
792	DO je = 1, iext
793	fbdata%pext(1,jo,1+je) = &
794	& surfdata%rext(jo,pext%ipoint(je))
795	END DO
796	END DO
797
798	! Write the obfbdata structure
799	CALL write_obfbdata( clfname, fbdata )
800
801	! Output some basic statistics
802	CALL obs_wri_stats( fbdata )
803
804	CALL dealloc_obfbdata( fbdata )
805
806	END SUBROUTINE obs_wri_surf
807
808	SUBROUTINE obs_wri_stats( fbdata )
809	!!-----------------------------------------------------------------------
810	!!
811	!! * ROUTINE obs_wri_stats *
812	!!
813	!! ** Purpose : Output some basic statistics of the data being written out
814	!!
815	!! ** Method :
816	!!
817	!! ** Action :
818	!!
819	!! ! 2014-08 (D. Lea) Initial version
820	!!-----------------------------------------------------------------------
821
822	!! * Arguments
823	TYPE(obfbdata) :: fbdata
824
825	!! * Local declarations
826	INTEGER :: jvar
827	INTEGER :: jo
828	INTEGER :: jk
829	INTEGER :: inumgoodobs
830	INTEGER :: inumgoodobsmpp
831	REAL(wp) :: zsumx
832	REAL(wp) :: zsumx2
833	REAL(wp) :: zomb
834
835
836	IF (lwp) THEN
837	WRITE(numout,*) ''
838	WRITE(numout,*) 'obs_wri_stats :'
839	WRITE(numout,*) '~~~~~~~~~~~~~~~'
840	ENDIF
841
842	DO jvar = 1, fbdata%nvar
843	zsumx=0.0_wp
844	zsumx2=0.0_wp
845	inumgoodobs=0
846	DO jo = 1, fbdata%nobs
847	DO jk = 1, fbdata%nlev
848	IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
849	& ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
850	& ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
851
852	zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
853	zsumx=zsumx+zomb
854	zsumx2=zsumx2+zomb**2
855	inumgoodobs=inumgoodobs+1
856	ENDIF
857	ENDDO
858	ENDDO
859
860	CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
861	CALL mpp_sum(zsumx)
862	CALL mpp_sum(zsumx2)
863
864	IF (lwp) THEN
865	WRITE(numout,*) 'Type: ',fbdata%cname(jvar),' Total number of good observations: ',inumgoodobsmpp
866	WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
867	WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
868	WRITE(numout,*) ''
869	ENDIF
870
871	ENDDO
872
873	END SUBROUTINE obs_wri_stats
874
875	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/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 11932

Download in other formats: