Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obs_write.F90 @ 13384

Last change on this file since 13384 was 13384, checked in by mattmartin, 4 years ago
First working version of surface velocity observation operator code.
File size: 31.5 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), 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('SSS')
545
546	clfiletype = 'sssfb'
547	cllongname(1) = 'Sea surface salinity'
548	clunits(1) = 'psu'
549	clgrid(1) = 'T'
550
551	CASE('UVEL')
552
553	clfiletype = 'ssvfb'
554	cllongname(1) = 'Eastward sea surface velocity'
555	clunits(1) = 'm s-1'
556	clgrid(1) = 'U'
557	cllongname(2) = 'Northward sea surface velocity'
558	clunits(2) = 'm s-1'
559	clgrid(2) = 'V'
560
561	CASE('SLCHLTOT')
562
563	clfiletype = 'slchltotfb'
564	cllongname(1) = 'Surface total log10(chlorophyll)'
565	clunits(1) = 'log10(mg/m3)'
566	clgrid(1) = 'T'
567
568	CASE('SLCHLDIA')
569
570	clfiletype = 'slchldiafb'
571	cllongname(1) = 'Surface diatom log10(chlorophyll)'
572	clunits(1) = 'log10(mg/m3)'
573	clgrid(1) = 'T'
574
575	CASE('SLCHLNON')
576
577	clfiletype = 'slchlnonfb'
578	cllongname(1) = 'Surface non-diatom log10(chlorophyll)'
579	clunits(1) = 'log10(mg/m3)'
580	clgrid(1) = 'T'
581
582	CASE('SLCHLDIN')
583
584	clfiletype = 'slchldinfb'
585	cllongname(1) = 'Surface dinoflagellate log10(chlorophyll)'
586	clunits(1) = 'log10(mg/m3)'
587	clgrid(1) = 'T'
588
589	CASE('SLCHLMIC')
590
591	clfiletype = 'slchlmicfb'
592	cllongname(1) = 'Surface microphytoplankton log10(chlorophyll)'
593	clunits(1) = 'log10(mg/m3)'
594	clgrid(1) = 'T'
595
596	CASE('SLCHLNAN')
597
598	clfiletype = 'slchlnanfb'
599	cllongname(1) = 'Surface nanophytoplankton log10(chlorophyll)'
600	clunits(1) = 'log10(mg/m3)'
601	clgrid(1) = 'T'
602
603	CASE('SLCHLPIC')
604
605	clfiletype = 'slchlpicfb'
606	cllongname(1) = 'Surface picophytoplankton log10(chlorophyll)'
607	clunits(1) = 'log10(mg/m3)'
608	clgrid(1) = 'T'
609
610	CASE('SCHLTOT')
611
612	clfiletype = 'schltotfb'
613	cllongname = 'Surface total chlorophyll'
614	clunits = 'mg/m3'
615	clgrid = 'T'
616
617	CASE('SLPHYTOT')
618
619	clfiletype = 'slphytotfb'
620	cllongname(1) = 'Surface total log10(phytoplankton carbon)'
621	clunits(1) = 'log10(mmolC/m3)'
622	clgrid(1) = 'T'
623
624	CASE('SLPHYDIA')
625
626	clfiletype = 'slphydiafb'
627	cllongname(1) = 'Surface diatom log10(phytoplankton carbon)'
628	clunits(1) = 'log10(mmolC/m3)'
629	clgrid(1) = 'T'
630
631	CASE('SLPHYNON')
632
633	clfiletype = 'slphynonfb'
634	cllongname(1) = 'Surface non-diatom log10(phytoplankton carbon)'
635	clunits(1) = 'log10(mmolC/m3)'
636	clgrid(1) = 'T'
637
638	CASE('SSPM')
639
640	clfiletype = 'sspmfb'
641	cllongname(1) = 'Surface suspended particulate matter'
642	clunits(1) = 'g/m3'
643	clgrid(1) = 'T'
644
645	CASE('SKD490')
646
647	clfiletype = 'skd490fb'
648	cllongname(1) = 'Surface attenuation coefficient of downwelling radiation at 490 nm'
649	clunits(1) = 'm-1'
650	clgrid(1) = 'T'
651
652	CASE('SFCO2')
653
654	clfiletype = 'sfco2fb'
655	cllongname(1) = 'Surface fugacity of carbon dioxide'
656	clunits(1) = 'uatm'
657	clgrid(1) = 'T'
658
659	CASE('SPCO2')
660
661	clfiletype = 'spco2fb'
662	cllongname(1) = 'Surface partial pressure of carbon dioxide'
663	clunits(1) = 'uatm'
664	clgrid(1) = 'T'
665
666	CASE DEFAULT
667
668	CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
669
670	END SELECT
671
672	! SLA needs special treatment because of MDT, so is done above
673	! Remaining variables treated more generically
674
675	IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
676
677	CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, &
678	& 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
679
680	DO jv = 1, surfdata%nvar
681	fbdata%cname(jv) = surfdata%cvars(jv)
682	fbdata%coblong(jv) = cllongname(jv)
683	fbdata%cobunit(jv) = clunits(jv)
684	END DO
685	DO je = 1, iext
686	fbdata%cextname(je) = pext%cdname(je)
687	fbdata%cextlong(je) = pext%cdlong(je,1)
688	fbdata%cextunit(je) = pext%cdunit(je,1)
689	END DO
690	DO jv = 1, surfdata%nvar
691	IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
692	fbdata%caddlong(1,jv) = 'Model interpolated ICE'
693	ELSE
694	fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv))
695	ENDIF
696	fbdata%caddunit(1,jv) = clunits(jv)
697	fbdata%cgrid(jv) = clgrid(jv)
698	END DO
699	DO ja = 1, iadd
700	fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)
701	DO jv = 1, surfdata%nvar
702	fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv)
703	fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv)
704	END DO
705	END DO
706	ENDIF
707
708	fbdata%caddname(1) = 'Hx'
709	IF ( indx_std /= -1 ) THEN
710	fbdata%caddname(1+iadd_mdt+iadd_std) = surfdata%cext(indx_std)
711	DO jv = 1, surfdata%nvar
712	fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'
713	fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)
714	END DO
715	ENDIF
716
717	IF ( surfdata%lclim ) THEN
718	fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm) = 'CLM'
719	DO jv = 1, surfdata%nvar
720	fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology'
721	fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1)
722	END DO
723	ENDIF
724
725	WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
726
727	IF(lwp) THEN
728	WRITE(numout,*)
729	WRITE(numout,*)'obs_wri_surf :'
730	WRITE(numout,*)'~~~~~~~~~~~~~'
731	WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
732	ENDIF
733
734	! Transform surf data structure into obfbdata structure
735	fbdata%cdjuldref = '19500101000000'
736	DO jo = 1, surfdata%nsurf
737	fbdata%plam(jo) = surfdata%rlam(jo)
738	fbdata%pphi(jo) = surfdata%rphi(jo)
739	WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
740	fbdata%ivqc(jo,:) = 0
741	fbdata%ivqcf(:,jo,:) = 0
742	IF ( surfdata%nqc(jo) > 255 ) THEN
743	fbdata%ioqc(jo) = 4
744	fbdata%ioqcf(1,jo) = 0
745	fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
746	ELSE
747	fbdata%ioqc(jo) = surfdata%nqc(jo)
748	fbdata%ioqcf(:,jo) = 0
749	ENDIF
750	fbdata%ipqc(jo) = 0
751	fbdata%ipqcf(:,jo) = 0
752	fbdata%itqc(jo) = 0
753	fbdata%itqcf(:,jo) = 0
754	fbdata%cdwmo(jo) = surfdata%cwmo(jo)
755	fbdata%kindex(jo) = surfdata%nsfil(jo)
756	IF (ln_grid_global) THEN
757	DO jv = 1, surfdata%nvar
758	fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv)
759	fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv)
760	END DO
761	ELSE
762	DO jv = 1, surfdata%nvar
763	fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv))
764	fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv))
765	END DO
766	ENDIF
767	CALL greg2jul( 0, &
768	& surfdata%nmin(jo), &
769	& surfdata%nhou(jo), &
770	& surfdata%nday(jo), &
771	& surfdata%nmon(jo), &
772	& surfdata%nyea(jo), &
773	& fbdata%ptim(jo), &
774	& krefdate = 19500101 )
775
776	fbdata%pdep(1,jo) = 0.0
777	fbdata%idqc(1,jo) = 0
778	fbdata%idqcf(:,1,jo) = 0
779	DO jv = 1, surfdata%nvar
780	fbdata%pob(1,jo,jv) = surfdata%robs(jo,jv)
781
782	IF ( surfdata%nqc(jo) > 255 ) THEN
783	fbdata%ivqc(jo,jv) = 4
784	fbdata%ivlqc(1,jo,jv) = 4
785	fbdata%ivlqcf(1,1,jo,jv) = 0
786	fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
787	ELSE
788	fbdata%ivqc(jo,jv) = surfdata%nqc(jo)
789	fbdata%ivlqc(1,jo,jv) = surfdata%nqc(jo)
790	fbdata%ivlqcf(:,1,jo,jv) = 0
791	ENDIF
792	fbdata%iobsk(1,jo,jv) = 0
793
794
795	! Additional variables.
796	! Hx is always the first additional variable
797	fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv)
798	! MDT is output as an additional variable if SLA obs type
799	IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN
800	fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1)
801	ENDIF
802	! STD is output as an additional variable if available
803	IF ( indx_std /= -1 ) THEN
804	fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std)
805	ENDIF
806	! CLM is output as an additional variable if available
807	IF ( surfdata%lclim ) THEN
808	fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv)
809	ENDIF
810	! Then other additional variables are output
811	DO ja = 1, iadd
812	fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = &
813	& surfdata%rext(jo,padd%ipoint(ja))
814	END DO
815	END DO
816	! Extra variables
817	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
818	DO je = 1, iext
819	fbdata%pext(1,jo,1+je) = &
820	& surfdata%rext(jo,pext%ipoint(je))
821	END DO
822	END DO
823
824	! Write the obfbdata structure
825	CALL write_obfbdata( clfname, fbdata )
826
827	! Output some basic statistics
828	CALL obs_wri_stats( fbdata )
829
830	CALL dealloc_obfbdata( fbdata )
831
832	END SUBROUTINE obs_wri_surf
833
834	SUBROUTINE obs_wri_stats( fbdata )
835	!!-----------------------------------------------------------------------
836	!!
837	!! * ROUTINE obs_wri_stats *
838	!!
839	!! ** Purpose : Output some basic statistics of the data being written out
840	!!
841	!! ** Method :
842	!!
843	!! ** Action :
844	!!
845	!! ! 2014-08 (D. Lea) Initial version
846	!!-----------------------------------------------------------------------
847
848	!! * Arguments
849	TYPE(obfbdata) :: fbdata
850
851	!! * Local declarations
852	INTEGER :: jvar
853	INTEGER :: jo
854	INTEGER :: jk
855	INTEGER :: inumgoodobs
856	INTEGER :: inumgoodobsmpp
857	REAL(wp) :: zsumx
858	REAL(wp) :: zsumx2
859	REAL(wp) :: zomb
860
861
862	IF (lwp) THEN
863	WRITE(numout,*) ''
864	WRITE(numout,*) 'obs_wri_stats :'
865	WRITE(numout,*) '~~~~~~~~~~~~~~~'
866	ENDIF
867
868	DO jvar = 1, fbdata%nvar
869	zsumx=0.0_wp
870	zsumx2=0.0_wp
871	inumgoodobs=0
872	DO jo = 1, fbdata%nobs
873	DO jk = 1, fbdata%nlev
874	IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
875	& ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
876	& ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
877
878	zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
879	zsumx=zsumx+zomb
880	zsumx2=zsumx2+zomb**2
881	inumgoodobs=inumgoodobs+1
882	ENDIF
883	ENDDO
884	ENDDO
885
886	CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
887	CALL mpp_sum(zsumx)
888	CALL mpp_sum(zsumx2)
889
890	IF (lwp) THEN
891	WRITE(numout,*) 'Type: ',fbdata%cname(jvar),' Total number of good observations: ',inumgoodobsmpp
892	WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
893	WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
894	WRITE(numout,*) ''
895	ENDIF
896
897	ENDDO
898
899	END SUBROUTINE obs_wri_stats
900
901	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_vel_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 13384

Download in other formats: