source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 11468

Last change on this file since 11468 was 11468, checked in by mattmartin, 17 months ago

Merged changes to allow writing of climatological information to feedback files.

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