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

source: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 11449

Last change on this file since 11449 was 11449, checked in by mattmartin, 5 years ago

Committed first version of changes to output climatology values at obs locations in the fdbk files.
File size: 29.6 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   USE diaobs, ONLY: &
32      & ln_output_clim
33
34   IMPLICIT NONE
35
36   !! * Routine accessibility
37   PRIVATE
38   PUBLIC obs_wri_prof, &    ! Write profile observation files
39      &   obs_wri_surf, &    ! Write surface observation files
40      &   obswriinfo
41   
42   TYPE obswriinfo
43      INTEGER :: inum
44      INTEGER, POINTER, DIMENSION(:) :: ipoint
45      CHARACTER(len=ilenname), POINTER, DIMENSION(:) :: cdname
46      CHARACTER(len=ilenlong), POINTER, DIMENSION(:,:) :: cdlong
47      CHARACTER(len=ilenunit), POINTER, DIMENSION(:,:) :: cdunit
48   END TYPE obswriinfo
49
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE obs_wri_prof( profdata, padd, pext )
59      !!-----------------------------------------------------------------------
60      !!
61      !!                     *** ROUTINE obs_wri_prof  ***
62      !!
63      !! ** Purpose : Write profile feedback files
64      !!
65      !! ** Method  : NetCDF
66      !!
67      !! ** Action  :
68      !!
69      !! History :
70      !!      ! 06-04  (A. Vidard) Original
71      !!      ! 06-04  (A. Vidard) Reformatted
72      !!      ! 06-10  (A. Weaver) Cleanup
73      !!      ! 07-01  (K. Mogensen) Use profile data types
74      !!      ! 07-03  (K. Mogensen) General handling of profiles
75      !!      ! 09-01  (K. Mogensen) New feedback format
76      !!      ! 15-02  (M. Martin) Combined routine for writing profiles
77      !!-----------------------------------------------------------------------
78
79      !! * Arguments
80      TYPE(obs_prof), INTENT(INOUT) :: profdata      ! Full set of profile data
81      TYPE(obswriinfo), OPTIONAL :: padd             ! Additional info for each variable
82      TYPE(obswriinfo), OPTIONAL :: pext             ! Extra info
83
84      !! * Local declarations
85      TYPE(obfbdata) :: fbdata
86      CHARACTER(LEN=40) :: clfname
87      CHARACTER(LEN=10) :: clfiletype
88      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
89      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
90      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
91      INTEGER :: ilevel
92      INTEGER :: jvar
93      INTEGER :: jo
94      INTEGER :: jk
95      INTEGER :: ik
96      INTEGER :: ja
97      INTEGER :: je
98      INTEGER :: iadd
99      INTEGER :: iadd_exp  ! expected additional variables
100      INTEGER :: iext
101      REAL(wp) :: zpres
102
103
104      ! Set up number of additional variables to be ouput:
105      ! Hx, CLIM, ...
106      iadd_exp = 1   ! Hx
107      IF ( ln_output_clim ) iadd_exp = iadd_exp + 1
108     
109      IF ( PRESENT( padd ) ) THEN
110         iadd = padd%inum
111      ELSE
112         iadd = 0
113      ENDIF
114
115      IF ( PRESENT( pext ) ) THEN
116         iext = pext%inum
117      ELSE
118         iext = 0
119      ENDIF
120
121      CALL init_obfbdata( fbdata )
122
123      ! Find maximum level
124      ilevel = 0
125      DO jvar = 1, profdata%nvar
126         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
127      END DO
128
129      SELECT CASE ( TRIM(profdata%cvars(1)) )
130      CASE('POTM')
131
132         clfiletype='profb'
133         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
134            &                 iadd_exp + iadd, 1 + iext, .TRUE. )
135         fbdata%cname(1)      = profdata%cvars(1)
136         fbdata%cname(2)      = profdata%cvars(2)
137         fbdata%coblong(1)    = 'Potential temperature'
138         fbdata%coblong(2)    = 'Practical salinity'
139         fbdata%cobunit(1)    = 'Degrees centigrade'
140         fbdata%cobunit(2)    = 'PSU'
141         
142         fbdata%cextname(1)   = 'TEMP'
143         fbdata%cextlong(1)   = 'Insitu temperature'
144         fbdata%cextunit(1)   = 'Degrees centigrade'
145         
146         fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
147         fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
148         fbdata%caddunit(1,1) = 'Degrees centigrade'
149         fbdata%caddunit(1,2) = 'PSU'
150         IF ( ln_output_clim ) THEN
151            fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature'
152            fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity'
153            fbdata%caddunit(2,1) = 'Degrees centigrade'
154            fbdata%caddunit(2,2) = 'PSU'
155         ENDIF
156         fbdata%cgrid(:)      = 'T'
157         DO je = 1, iext
158            fbdata%cextname(1+je) = pext%cdname(je)
159            fbdata%cextlong(1+je) = pext%cdlong(je,1)
160            fbdata%cextunit(1+je) = pext%cdunit(je,1)
161         END DO
162         DO ja = 1, iadd
163            fbdata%caddname(iadd_exp+ja) = padd%cdname(ja)
164            DO jvar = 1, 2
165               fbdata%caddlong(iadd_exp+ja,jvar) = padd%cdlong(ja,jvar)
166               fbdata%caddunit(iadd_exp+ja,jvar) = padd%cdunit(ja,jvar)
167            END DO
168         END DO
169
170      CASE('UVEL')
171
172         clfiletype='velfb'
173         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
174            &                 iadd_exp + iadd, 0, .TRUE. )
175         fbdata%cname(1)      = profdata%cvars(1)
176         fbdata%cname(2)      = profdata%cvars(2)
177         fbdata%coblong(1)    = 'Zonal velocity'
178         fbdata%coblong(2)    = 'Meridional velocity'
179         fbdata%cobunit(1)    = 'm/s'
180         fbdata%cobunit(2)    = 'm/s'
181         DO je = 1, iext
182            fbdata%cextname(je) = pext%cdname(je)
183            fbdata%cextlong(je) = pext%cdlong(je,1)
184            fbdata%cextunit(je) = pext%cdunit(je,1)
185         END DO
186         fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
187         fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
188         fbdata%caddunit(1,1) = 'm/s'
189         fbdata%caddunit(1,2) = 'm/s'
190         IF ( ln_output_clim ) THEN
191            fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity'
192            fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity'
193            fbdata%caddunit(2,1) = 'm/s'
194            fbdata%caddunit(2,2) = 'm/s'
195         ENDIF         
196         fbdata%cgrid(1)      = 'U' 
197         fbdata%cgrid(2)      = 'V'
198         DO ja = 1, iadd
199            fbdata%caddname(iadd_exp+ja) = padd%cdname(ja)
200            fbdata%caddlong(iadd_exp+ja,1) = padd%cdlong(ja,1)
201            fbdata%caddunit(iadd_exp+ja,1) = padd%cdunit(ja,1)
202         END DO
203
204      CASE('PLCHLTOT')
205
206         clfiletype = 'plchltotfb'
207         cllongname = 'log10(chlorophyll concentration)'
208         clunits    = 'log10(mg/m3)'
209         clgrid     = 'T'
210
211      CASE('PCHLTOT')
212
213         clfiletype = 'pchltotfb'
214         cllongname = 'chlorophyll concentration'
215         clunits    = 'mg/m3'
216         clgrid     = 'T'
217
218      CASE('PNO3')
219
220         clfiletype = 'pno3fb'
221         cllongname = 'nitrate'
222         clunits    = 'mmol/m3'
223         clgrid     = 'T'
224
225      CASE('PSI4')
226
227         clfiletype = 'psi4fb'
228         cllongname = 'silicate'
229         clunits    = 'mmol/m3'
230         clgrid     = 'T'
231
232      CASE('PPO4')
233
234         clfiletype = 'ppo4fb'
235         cllongname = 'phosphate'
236         clunits    = 'mmol/m3'
237         clgrid     = 'T'
238
239      CASE('PDIC')
240
241         clfiletype = 'pdicfb'
242         cllongname = 'dissolved inorganic carbon'
243         clunits    = 'mmol/m3'
244         clgrid     = 'T'
245
246      CASE('PALK')
247
248         clfiletype = 'palkfb'
249         cllongname = 'alkalinity'
250         clunits    = 'meq/m3'
251         clgrid     = 'T'
252
253      CASE('PPH')
254
255         clfiletype = 'pphfb'
256         cllongname = 'pH'
257         clunits    = '-'
258         clgrid     = 'T'
259
260      CASE('PO2')
261
262         clfiletype = 'po2fb'
263         cllongname = 'dissolved oxygen'
264         clunits    = 'mmol/m3'
265         clgrid     = 'T'
266
267      END SELECT
268     
269      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. &
270         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
271         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
272            &                 iadd_expt + iadd, iext, .TRUE. )
273         fbdata%cname(1)      = profdata%cvars(1)
274         fbdata%coblong(1)    = cllongname
275         fbdata%cobunit(1)    = clunits
276         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
277         fbdata%caddunit(1,1) = clunits
278         IF ( ln_output_clim ) THEN
279            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname)
280            fbdata%caddunit(2,1) = clunits
281         ENDIF         
282         fbdata%cgrid(:)      = clgrid
283         DO je = 1, iext
284            fbdata%cextname(je) = pext%cdname(je)
285            fbdata%cextlong(je) = pext%cdlong(je,1)
286            fbdata%cextunit(je) = pext%cdunit(je,1)
287         END DO
288         DO ja = 1, iadd
289            fbdata%caddname(iadd_expt+ja) = padd%cdname(ja)
290            fbdata%caddlong(iadd_expt+ja,1) = padd%cdlong(ja,1)
291            fbdata%caddunit(iadd_expt+ja,1) = padd%cdunit(ja,1)
292         END DO
293      ENDIF
294
295      fbdata%caddname(1)   = 'Hx'
296      IF ( ln_output_clim ) fbdata%caddname(2)   = 'CLM'
297     
298      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
299
300      IF(lwp) THEN
301         WRITE(numout,*)
302         WRITE(numout,*)'obs_wri_prof :'
303         WRITE(numout,*)'~~~~~~~~~~~~~'
304         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
305      ENDIF
306
307      ! Transform obs_prof data structure into obfb data structure
308      fbdata%cdjuldref = '19500101000000'
309      DO jo = 1, profdata%nprof
310         fbdata%plam(jo)      = profdata%rlam(jo)
311         fbdata%pphi(jo)      = profdata%rphi(jo)
312         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
313         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
314         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
315         IF ( profdata%nqc(jo) > 255 ) THEN
316            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2)
317            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
318            fbdata%ioqcf(2,jo) = profdata%nqc(jo)
319         ELSE
320            fbdata%ioqc(jo)    = profdata%nqc(jo)
321            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
322         ENDIF
323         fbdata%ipqc(jo)      = profdata%ipqc(jo)
324         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
325         fbdata%itqc(jo)      = profdata%itqc(jo)
326         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
327         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
328         fbdata%kindex(jo)    = profdata%npfil(jo)
329         DO jvar = 1, profdata%nvar
330            IF (ln_grid_global) THEN
331               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
332               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
333            ELSE
334               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
335               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
336            ENDIF
337         END DO
338         CALL greg2jul( 0, &
339            &           profdata%nmin(jo), &
340            &           profdata%nhou(jo), &
341            &           profdata%nday(jo), &
342            &           profdata%nmon(jo), &
343            &           profdata%nyea(jo), &
344            &           fbdata%ptim(jo),   &
345            &           krefdate = 19500101 )
346         ! Reform the profiles arrays for output
347         DO jvar = 1, profdata%nvar
348            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
349               ik = profdata%var(jvar)%nvlidx(jk)
350               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
351               IF ( ln_output_clim ) THEN           
352                  fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vclm(jk)     
353               ENDIF             
354               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
355               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
356               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
357               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
358               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
359                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
360                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
361                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111')
362               ELSE
363                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
364                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
365               ENDIF
366               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
367               DO ja = 1, iadd
368                  fbdata%padd(ik,jo,iadd_exp+ja,jvar) = &
369                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
370               END DO
371               DO je = 1, iext
372                  fbdata%pext(ik,jo,1+je) = &
373                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
374               END DO
375               IF ( ( jvar == 1 ) .AND. &
376                  & ( TRIM(profdata%cvars(1)) == 'POTM' ) ) THEN
377                  fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
378               ENDIF
379            END DO
380         END DO
381      END DO
382
383      IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
384         ! Convert insitu temperature to potential temperature using the model
385         ! salinity if no potential temperature
386         DO jo = 1, fbdata%nobs
387            IF ( fbdata%pphi(jo) < 9999.0 ) THEN
388               DO jk = 1, fbdata%nlev
389                  IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
390                     & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
391                     & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
392                     & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
393                     zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
394                        &              REAL(fbdata%pphi(jo),wp) )
395                     fbdata%pob(jk,jo,1) = potemp( &
396                        &                     REAL(fbdata%padd(jk,jo,1,2), wp),  &
397                        &                     REAL(fbdata%pext(jk,jo,1), wp), &
398                        &                     zpres, 0.0_wp )
399                  ENDIF
400               END DO
401            ENDIF
402         END DO
403      ENDIF
404
405      ! Write the obfbdata structure
406      CALL write_obfbdata( clfname, fbdata )
407
408      ! Output some basic statistics
409      CALL obs_wri_stats( fbdata )
410
411      CALL dealloc_obfbdata( fbdata )
412
413   END SUBROUTINE obs_wri_prof
414
415   SUBROUTINE obs_wri_surf( surfdata, padd, pext )
416      !!-----------------------------------------------------------------------
417      !!
418      !!                     *** ROUTINE obs_wri_surf  ***
419      !!
420      !! ** Purpose : Write surface observation files
421      !!
422      !! ** Method  : NetCDF
423      !!
424      !! ** Action  :
425      !!
426      !!      ! 07-03  (K. Mogensen) Original
427      !!      ! 09-01  (K. Mogensen) New feedback format.
428      !!      ! 15-02  (M. Martin) Combined surface writing routine.
429      !!-----------------------------------------------------------------------
430
431      !! * Modules used
432      IMPLICIT NONE
433
434      !! * Arguments
435      TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data
436      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
437      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info
438
439      !! * Local declarations
440      TYPE(obfbdata) :: fbdata
441      CHARACTER(LEN=40) :: clfname         ! netCDF filename
442      CHARACTER(LEN=10) :: clfiletype
443      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
444      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
445      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
446      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
447      INTEGER :: jo
448      INTEGER :: ja
449      INTEGER :: je
450      INTEGER :: iadd
451      INTEGER :: iadd_exp
452      INTEGER :: iext
453      INTEGER :: indx_std
454      INTEGER :: iadd_std
455      INTEGER :: iadd_clm     
456
457
458      ! Set up number of additional variables to be ouput:
459      ! Hx, CLIM, ...
460      iadd_exp = 1   ! Hx
461      IF ( ln_output_clim ) iadd_exp = iadd_exp + 1
462 
463      IF ( PRESENT( padd ) ) THEN
464         iadd = padd%inum
465      ELSE
466         iadd = 0
467      ENDIF
468
469      IF ( PRESENT( pext ) ) THEN
470         iext = pext%inum
471      ELSE
472         iext = 0
473      ENDIF
474
475      iadd_std = 0
476      indx_std = -1
477      IF ( surfdata%nextra > 0 ) THEN
478         DO je = 1, surfdata%nextra
479           IF ( TRIM( surfdata%cext(je) ) == 'STD' ) THEN
480             iadd_std = 1
481             indx_std = je
482           ENDIF
483         END DO
484      ENDIF
485     
486      iadd_clm = 0
487     
488      CALL init_obfbdata( fbdata )
489
490      SELECT CASE ( TRIM(surfdata%cvars(1)) )
491      CASE('SLA')
492         
493         ! SLA needs special treatment because of MDT, so is all done here
494         ! Other variables are done more generically
495         ! No climatology for SLA, MDT is our best estimate of that and is already output.
496
497         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
498            &                 2 + iadd + iadd_std, 1 + iext, .TRUE. )
499
500         clfiletype = 'slafb'
501         fbdata%cname(1)      = surfdata%cvars(1)
502         fbdata%coblong(1)    = 'Sea level anomaly'
503         fbdata%cobunit(1)    = 'Metres'
504         fbdata%cextname(1)   = 'MDT'
505         fbdata%cextlong(1)   = 'Mean dynamic topography'
506         fbdata%cextunit(1)   = 'Metres'
507         DO je = 1, iext
508            fbdata%cextname(je) = pext%cdname(je)
509            fbdata%cextlong(je) = pext%cdlong(je,1)
510            fbdata%cextunit(je) = pext%cdunit(je,1)
511         END DO
512         fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
513         fbdata%caddunit(1,1) = 'Metres' 
514         fbdata%caddname(2)   = 'SSH'
515         fbdata%caddlong(2,1) = 'Model Sea surface height'
516         fbdata%caddunit(2,1) = 'Metres'
517         fbdata%cgrid(1)      = 'T'
518         DO ja = 1, iadd
519            fbdata%caddname(2+iadd_std+ja) = padd%cdname(ja)
520            fbdata%caddlong(2+iadd_std+ja,1) = padd%cdlong(ja,1)
521            fbdata%caddunit(2+iadd_std+ja,1) = padd%cdunit(ja,1)
522         END DO
523
524      CASE('SST')
525
526         clfiletype = 'sstfb'
527         cllongname = 'Sea surface temperature'
528         clunits    = 'Degree centigrade'
529         clgrid     = 'T'
530         IF ( ln_output_clim ) iadd_clm = 1
531         
532      CASE('ICECONC')
533
534         clfiletype = 'sicfb'
535         cllongname = 'Sea ice'
536         clunits    = 'Fraction'
537         clgrid     = 'T'
538
539      CASE('SSS')
540
541         clfiletype = 'sssfb'
542         cllongname = 'Sea surface salinity'
543         clunits    = 'psu'
544         clgrid     = 'T'
545         IF ( ln_output_clim ) iadd_clm = 1
546         
547      CASE('SLCHLTOT','LOGCHL','LogChl','logchl')
548
549         clfiletype = 'slchltotfb'
550         cllongname = 'Surface total log10(chlorophyll)'
551         clunits    = 'log10(mg/m3)'
552         clgrid     = 'T'
553
554      CASE('SLCHLDIA')
555
556         clfiletype = 'slchldiafb'
557         cllongname = 'Surface diatom log10(chlorophyll)'
558         clunits    = 'log10(mg/m3)'
559         clgrid     = 'T'
560
561      CASE('SLCHLNON')
562
563         clfiletype = 'slchlnonfb'
564         cllongname = 'Surface non-diatom log10(chlorophyll)'
565         clunits    = 'log10(mg/m3)'
566         clgrid     = 'T'
567
568      CASE('SLCHLDIN')
569
570         clfiletype = 'slchldinfb'
571         cllongname = 'Surface dinoflagellate log10(chlorophyll)'
572         clunits    = 'log10(mg/m3)'
573         clgrid     = 'T'
574
575      CASE('SLCHLMIC')
576
577         clfiletype = 'slchlmicfb'
578         cllongname = 'Surface microphytoplankton log10(chlorophyll)'
579         clunits    = 'log10(mg/m3)'
580         clgrid     = 'T'
581
582      CASE('SLCHLNAN')
583
584         clfiletype = 'slchlnanfb'
585         cllongname = 'Surface nanophytoplankton log10(chlorophyll)'
586         clunits    = 'log10(mg/m3)'
587         clgrid     = 'T'
588
589      CASE('SLCHLPIC')
590
591         clfiletype = 'slchlpicfb'
592         cllongname = 'Surface picophytoplankton log10(chlorophyll)'
593         clunits    = 'log10(mg/m3)'
594         clgrid     = 'T'
595
596      CASE('SCHLTOT')
597
598         clfiletype = 'schltotfb'
599         cllongname = 'Surface total chlorophyll'
600         clunits    = 'mg/m3'
601         clgrid     = 'T'
602
603      CASE('SLPHYTOT')
604
605         clfiletype = 'slphytotfb'
606         cllongname = 'Surface total log10(phytoplankton carbon)'
607         clunits    = 'log10(mmolC/m3)'
608         clgrid     = 'T'
609
610      CASE('SLPHYDIA')
611
612         clfiletype = 'slphydiafb'
613         cllongname = 'Surface diatom log10(phytoplankton carbon)'
614         clunits    = 'log10(mmolC/m3)'
615         clgrid     = 'T'
616
617      CASE('SLPHYNON')
618
619         clfiletype = 'slphynonfb'
620         cllongname = 'Surface non-diatom log10(phytoplankton carbon)'
621         clunits    = 'log10(mmolC/m3)'
622         clgrid     = 'T'
623
624      CASE('SSPM')
625
626         clfiletype = 'sspmfb'
627         cllongname = 'Surface suspended particulate matter'
628         clunits    = 'g/m3'
629         clgrid     = 'T'
630
631      CASE('SFCO2','FCO2','fCO2','fco2')
632
633         clfiletype = 'sfco2fb'
634         cllongname = 'Surface fugacity of carbon dioxide'
635         clunits    = 'uatm'
636         clgrid     = 'T'
637
638      CASE('SPCO2')
639
640         clfiletype = 'spco2fb'
641         cllongname = 'Surface partial pressure of carbon dioxide'
642         clunits    = 'uatm'
643         clgrid     = 'T'
644
645      CASE DEFAULT
646
647         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
648
649      END SELECT
650
651      ! SLA needs special treatment because of MDT, so is done above
652      ! Remaining variables treated more generically
653
654      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
655     
656         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
657            &                 1 + iadd + iadd_std + iadd_clm, iext, .TRUE. )
658
659         fbdata%cname(1)      = surfdata%cvars(1)
660         fbdata%coblong(1)    = cllongname
661         fbdata%cobunit(1)    = clunits
662         DO je = 1, iext
663            fbdata%cextname(je) = pext%cdname(je)
664            fbdata%cextlong(je) = pext%cdlong(je,1)
665            fbdata%cextunit(je) = pext%cdunit(je,1)
666         END DO
667         IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
668            fbdata%caddlong(1,1) = 'Model interpolated ICE'
669         ELSE
670            fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1))
671         ENDIF
672         fbdata%caddunit(1,1) = clunits
673         fbdata%cgrid(1)      = clgrid
674         DO ja = 1, iadd
675            fbdata%caddname(1+iadd_std+iadd_clm+ja) = padd%cdname(ja)
676            fbdata%caddlong(1+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1)
677            fbdata%caddunit(1+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1)
678         END DO
679
680      ENDIF
681     
682      fbdata%caddname(1)   = 'Hx'
683      IF ( indx_std /= -1 ) THEN
684         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1
685         fbdata%caddname(1+iadd_std)   = surfdata%cext(indx_std)
686         fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation'
687         fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1)
688      ENDIF
689     
690      IF ( ln_output_clim .AND. ( iadd_clm > 0 ) ) THEN
691         fbdata%caddname(1+iadd_std+iadd_clm)   = 'CLM'
692         fbdata%caddlong(1+iadd_std+iadd_clm,1) = 'Climatology'
693         fbdata%caddunit(1+iadd_std+iadd_clm,1) = fbdata%cobunit(1)
694      ENDIF
695     
696     
697      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
698
699      IF(lwp) THEN
700         WRITE(numout,*)
701         WRITE(numout,*)'obs_wri_surf :'
702         WRITE(numout,*)'~~~~~~~~~~~~~'
703         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
704      ENDIF
705
706      ! Transform surf data structure into obfbdata structure
707      fbdata%cdjuldref = '19500101000000'
708      DO jo = 1, surfdata%nsurf
709         fbdata%plam(jo)      = surfdata%rlam(jo)
710         fbdata%pphi(jo)      = surfdata%rphi(jo)
711         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
712         fbdata%ivqc(jo,:)    = 0
713         fbdata%ivqcf(:,jo,:) = 0
714         IF ( surfdata%nqc(jo) > 255 ) THEN
715            fbdata%ioqc(jo)    = 4
716            fbdata%ioqcf(1,jo) = 0
717            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
718         ELSE
719            fbdata%ioqc(jo)    = surfdata%nqc(jo)
720            fbdata%ioqcf(:,jo) = 0
721         ENDIF
722         fbdata%ipqc(jo)      = 0
723         fbdata%ipqcf(:,jo)   = 0
724         fbdata%itqc(jo)      = 0
725         fbdata%itqcf(:,jo)   = 0
726         fbdata%cdwmo(jo)     = surfdata%cwmo(jo)
727         fbdata%kindex(jo)    = surfdata%nsfil(jo)
728         IF (ln_grid_global) THEN
729            fbdata%iobsi(jo,1) = surfdata%mi(jo)
730            fbdata%iobsj(jo,1) = surfdata%mj(jo)
731         ELSE
732            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
733            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
734         ENDIF
735         CALL greg2jul( 0, &
736            &           surfdata%nmin(jo), &
737            &           surfdata%nhou(jo), &
738            &           surfdata%nday(jo), &
739            &           surfdata%nmon(jo), &
740            &           surfdata%nyea(jo), &
741            &           fbdata%ptim(jo),   &
742            &           krefdate = 19500101 )
743         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
744         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
745            fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
746         ENDIF   
747         IF ( ln_output_clim .AND. ( iadd_clm > 0 ) ) THEN
748            fbdata%padd(1,jo,2,1) = surfdata%rclm(jo,1)
749         ENDIF
750                     
751         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
752         fbdata%pdep(1,jo)     = 0.0
753         fbdata%idqc(1,jo)     = 0
754         fbdata%idqcf(:,1,jo)  = 0
755         IF ( surfdata%nqc(jo) > 255 ) THEN
756            fbdata%ivqc(jo,1)       = 4
757            fbdata%ivlqc(1,jo,1)    = 4
758            fbdata%ivlqcf(1,1,jo,1) = 0
759            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
760         ELSE
761            fbdata%ivqc(jo,1)       = surfdata%nqc(jo)
762            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo)
763            fbdata%ivlqcf(:,1,jo,1) = 0
764         ENDIF
765         fbdata%iobsk(1,jo,1)  = 0
766         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
767         IF ( indx_std /= -1 ) THEN
768            fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std)
769         ENDIF
770         
771         DO ja = 1, iadd
772            fbdata%padd(1,jo,2+iadd_std+ja,1) = &
773               & surfdata%rext(jo,padd%ipoint(ja))
774         END DO
775         DO je = 1, iext
776            fbdata%pext(1,jo,1+je) = &
777               & surfdata%rext(jo,pext%ipoint(je))
778         END DO
779      END DO
780
781      ! Write the obfbdata structure
782      CALL write_obfbdata( clfname, fbdata )
783
784      ! Output some basic statistics
785      CALL obs_wri_stats( fbdata )
786
787      CALL dealloc_obfbdata( fbdata )
788
789   END SUBROUTINE obs_wri_surf
790
791   SUBROUTINE obs_wri_stats( fbdata )
792      !!-----------------------------------------------------------------------
793      !!
794      !!                     *** ROUTINE obs_wri_stats  ***
795      !!
796      !! ** Purpose : Output some basic statistics of the data being written out
797      !!
798      !! ** Method  :
799      !!
800      !! ** Action  :
801      !!
802      !!      ! 2014-08  (D. Lea) Initial version
803      !!-----------------------------------------------------------------------
804
805      !! * Arguments
806      TYPE(obfbdata) :: fbdata
807
808      !! * Local declarations
809      INTEGER :: jvar
810      INTEGER :: jo
811      INTEGER :: jk
812      INTEGER :: inumgoodobs
813      INTEGER :: inumgoodobsmpp
814      REAL(wp) :: zsumx
815      REAL(wp) :: zsumx2
816      REAL(wp) :: zomb
817     
818
819      IF (lwp) THEN
820         WRITE(numout,*) ''
821         WRITE(numout,*) 'obs_wri_stats :'
822         WRITE(numout,*) '~~~~~~~~~~~~~~~'
823      ENDIF
824
825      DO jvar = 1, fbdata%nvar
826         zsumx=0.0_wp
827         zsumx2=0.0_wp
828         inumgoodobs=0
829         DO jo = 1, fbdata%nobs
830            DO jk = 1, fbdata%nlev
831               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
832                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
833                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
834
835                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
836                  zsumx=zsumx+zomb
837                  zsumx2=zsumx2+zomb**2
838                  inumgoodobs=inumgoodobs+1
839               ENDIF
840            ENDDO
841         ENDDO
842
843         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
844         CALL mpp_sum(zsumx)
845         CALL mpp_sum(zsumx2)
846
847         IF (lwp) THEN
848            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
849            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
850            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
851            WRITE(numout,*) ''
852         ENDIF
853
854      ENDDO
855
856   END SUBROUTINE obs_wri_stats
857
858END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.