New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

obs_write.F90 in branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 12820

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: