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/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

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

Last change on this file since 11546 was 11546, checked in by dford, 5 years ago

Add an observation operator for Kd490, and add some error checking so that the observation variable name read in from a feedback file matches a defined name for each observation type. Previously it could read in an incorrect variable without failure. See ​https://code.metoffice.gov.uk/trac/utils/ticket/247.
File size: 30.0 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')
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('SKD490')
628
629         clfiletype = 'skd490fb'
630         cllongname = 'Surface attenuation coefficient of downwelling radiation at 490 nm'
631         clunits    = 'm-1'
632         clgrid     = 'T'
633
634      CASE('SFCO2')
635
636         clfiletype = 'sfco2fb'
637         cllongname = 'Surface fugacity of carbon dioxide'
638         clunits    = 'uatm'
639         clgrid     = 'T'
640
641      CASE('SPCO2')
642
643         clfiletype = 'spco2fb'
644         cllongname = 'Surface partial pressure of carbon dioxide'
645         clunits    = 'uatm'
646         clgrid     = 'T'
647
648      CASE DEFAULT
649
650         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
651
652      END SELECT
653
654      ! SLA needs special treatment because of MDT, so is done above
655      ! Remaining variables treated more generically
656
657      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
658     
659         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
660            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. )
661
662         fbdata%cname(1)      = surfdata%cvars(1)
663         fbdata%coblong(1)    = cllongname
664         fbdata%cobunit(1)    = clunits
665         DO je = 1, iext
666            fbdata%cextname(je) = pext%cdname(je)
667            fbdata%cextlong(je) = pext%cdlong(je,1)
668            fbdata%cextunit(je) = pext%cdunit(je,1)
669         END DO
670         IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
671            fbdata%caddlong(1,1) = 'Model interpolated ICE'
672         ELSE
673            fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1))
674         ENDIF
675         fbdata%caddunit(1,1) = clunits
676         fbdata%cgrid(1)      = clgrid
677         DO ja = 1, iadd
678            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja)
679            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1)
680            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1)
681         END DO
682
683      ENDIF
684     
685      fbdata%caddname(1)   = 'Hx'
686      IF ( indx_std /= -1 ) THEN
687         fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std)
688         fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation'
689         fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1)
690      ENDIF
691     
692      IF ( surfdata%lclim ) THEN
693         fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm)   = 'CLM'
694         fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,1) = 'Climatology'
695         fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,1) = fbdata%cobunit(1)
696      ENDIF
697     
698     
699      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
700
701      IF(lwp) THEN
702         WRITE(numout,*)
703         WRITE(numout,*)'obs_wri_surf :'
704         WRITE(numout,*)'~~~~~~~~~~~~~'
705         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
706      ENDIF
707
708      ! Transform surf data structure into obfbdata structure
709      fbdata%cdjuldref = '19500101000000'
710      DO jo = 1, surfdata%nsurf
711         fbdata%plam(jo)      = surfdata%rlam(jo)
712         fbdata%pphi(jo)      = surfdata%rphi(jo)
713         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
714         fbdata%ivqc(jo,:)    = 0
715         fbdata%ivqcf(:,jo,:) = 0
716         IF ( surfdata%nqc(jo) > 255 ) THEN
717            fbdata%ioqc(jo)    = 4
718            fbdata%ioqcf(1,jo) = 0
719            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
720         ELSE
721            fbdata%ioqc(jo)    = surfdata%nqc(jo)
722            fbdata%ioqcf(:,jo) = 0
723         ENDIF
724         fbdata%ipqc(jo)      = 0
725         fbdata%ipqcf(:,jo)   = 0
726         fbdata%itqc(jo)      = 0
727         fbdata%itqcf(:,jo)   = 0
728         fbdata%cdwmo(jo)     = surfdata%cwmo(jo)
729         fbdata%kindex(jo)    = surfdata%nsfil(jo)
730         IF (ln_grid_global) THEN
731            fbdata%iobsi(jo,1) = surfdata%mi(jo)
732            fbdata%iobsj(jo,1) = surfdata%mj(jo)
733         ELSE
734            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
735            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
736         ENDIF
737         CALL greg2jul( 0, &
738            &           surfdata%nmin(jo), &
739            &           surfdata%nhou(jo), &
740            &           surfdata%nday(jo), &
741            &           surfdata%nmon(jo), &
742            &           surfdata%nyea(jo), &
743            &           fbdata%ptim(jo),   &
744            &           krefdate = 19500101 )
745                   
746         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
747         fbdata%pdep(1,jo)     = 0.0
748         fbdata%idqc(1,jo)     = 0
749         fbdata%idqcf(:,1,jo)  = 0
750         IF ( surfdata%nqc(jo) > 255 ) THEN
751            fbdata%ivqc(jo,1)       = 4
752            fbdata%ivlqc(1,jo,1)    = 4
753            fbdata%ivlqcf(1,1,jo,1) = 0
754            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
755         ELSE
756            fbdata%ivqc(jo,1)       = surfdata%nqc(jo)
757            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo)
758            fbdata%ivlqcf(:,1,jo,1) = 0
759         ENDIF
760         fbdata%iobsk(1,jo,1)  = 0
761 
762         ! Additional variables.
763         ! Hx is always the first additional variable
764         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
765         ! MDT is output as an additional variable if SLA obs type
766         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
767            fbdata%padd(1,jo,1+iadd_mdt,1) = surfdata%rext(jo,1)
768         ENDIF   
769         ! STD is output as an additional variable if available
770         IF ( indx_std /= -1 ) THEN
771            fbdata%padd(1,jo,1+iadd_mdt+iadd_std,1) = surfdata%rext(jo,indx_std)
772         ENDIF
773         ! CLM is output as an additional variable if available
774         IF ( surfdata%lclim ) THEN
775            fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,1) = surfdata%rclm(jo,1)
776         ENDIF
777         ! Then other additional variables are output
778         DO ja = 1, iadd
779            fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,1) = &
780               & surfdata%rext(jo,padd%ipoint(ja))
781         END DO
782         
783         ! Extra variables
784         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)         
785         DO je = 1, iext
786            fbdata%pext(1,jo,1+je) = &
787               & surfdata%rext(jo,pext%ipoint(je))
788         END DO
789      END DO
790
791      ! Write the obfbdata structure
792      CALL write_obfbdata( clfname, fbdata )
793
794      ! Output some basic statistics
795      CALL obs_wri_stats( fbdata )
796
797      CALL dealloc_obfbdata( fbdata )
798
799   END SUBROUTINE obs_wri_surf
800
801   SUBROUTINE obs_wri_stats( fbdata )
802      !!-----------------------------------------------------------------------
803      !!
804      !!                     *** ROUTINE obs_wri_stats  ***
805      !!
806      !! ** Purpose : Output some basic statistics of the data being written out
807      !!
808      !! ** Method  :
809      !!
810      !! ** Action  :
811      !!
812      !!      ! 2014-08  (D. Lea) Initial version
813      !!-----------------------------------------------------------------------
814
815      !! * Arguments
816      TYPE(obfbdata) :: fbdata
817
818      !! * Local declarations
819      INTEGER :: jvar
820      INTEGER :: jo
821      INTEGER :: jk
822      INTEGER :: inumgoodobs
823      INTEGER :: inumgoodobsmpp
824      REAL(wp) :: zsumx
825      REAL(wp) :: zsumx2
826      REAL(wp) :: zomb
827     
828
829      IF (lwp) THEN
830         WRITE(numout,*) ''
831         WRITE(numout,*) 'obs_wri_stats :'
832         WRITE(numout,*) '~~~~~~~~~~~~~~~'
833      ENDIF
834
835      DO jvar = 1, fbdata%nvar
836         zsumx=0.0_wp
837         zsumx2=0.0_wp
838         inumgoodobs=0
839         DO jo = 1, fbdata%nobs
840            DO jk = 1, fbdata%nlev
841               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
842                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
843                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
844
845                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
846                  zsumx=zsumx+zomb
847                  zsumx2=zsumx2+zomb**2
848                  inumgoodobs=inumgoodobs+1
849               ENDIF
850            ENDDO
851         ENDDO
852
853         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
854         CALL mpp_sum(zsumx)
855         CALL mpp_sum(zsumx2)
856
857         IF (lwp) THEN
858            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
859            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
860            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
861            WRITE(numout,*) ''
862         ENDIF
863
864      ENDDO
865
866   END SUBROUTINE obs_wri_stats
867
868END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.