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

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 5682

Last change on this file since 5682 was 5682, checked in by mattmartin, 9 years ago

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.
  • Property svn:keywords set to Id
File size: 21.7 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_rot_vel          ! Rotation of velocities
30   USE obs_mpp              ! MPP support routines for observation diagnostics
31   USE lib_mpp        ! MPP routines
32
33   IMPLICIT NONE
34
35   !! * Routine accessibility
36   PRIVATE
37   PUBLIC obs_wri_prof, &    ! Write profile observation files
38      &   obs_wri_surf, &    ! Write surface observation files
39      &   obswriinfo
40   
41   TYPE obswriinfo
42      INTEGER :: inum
43      INTEGER, POINTER, DIMENSION(:) :: ipoint
44      CHARACTER(len=ilenname), POINTER, DIMENSION(:) :: cdname
45      CHARACTER(len=ilenlong), POINTER, DIMENSION(:,:) :: cdlong
46      CHARACTER(len=ilenunit), POINTER, DIMENSION(:,:) :: cdunit
47   END TYPE obswriinfo
48
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE obs_wri_prof( profdata, k2dint, padd, pext )
58      !!-----------------------------------------------------------------------
59      !!
60      !!                     *** ROUTINE obs_wri_prof  ***
61      !!
62      !! ** Purpose : Write profile feedback files
63      !!
64      !! ** Method  : NetCDF
65      !!
66      !! ** Action  :
67      !!
68      !! History :
69      !!      ! 06-04  (A. Vidard) Original
70      !!      ! 06-04  (A. Vidard) Reformatted
71      !!      ! 06-10  (A. Weaver) Cleanup
72      !!      ! 07-01  (K. Mogensen) Use profile data types
73      !!      ! 07-03  (K. Mogensen) General handling of profiles
74      !!      ! 09-01  (K. Mogensen) New feedback format
75      !!      ! 15-02  (M. Martin) Combined routine for writing profiles
76      !!-----------------------------------------------------------------------
77
78      !! * Arguments
79      TYPE(obs_prof), INTENT(INOUT) :: profdata      ! Full set of profile data
80      INTEGER, INTENT(IN)        :: k2dint           ! Horizontal interpolation method
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=6) :: clfiletype
88      INTEGER :: ilevel
89      INTEGER :: jvar
90      INTEGER :: jo
91      INTEGER :: jk
92      INTEGER :: ik
93      INTEGER :: ja
94      INTEGER :: je
95      INTEGER :: iadd
96      INTEGER :: iext
97      REAL(wp) :: zpres
98      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
99         & zu, &
100         & zv
101
102      IF ( PRESENT( padd ) ) THEN
103         iadd = padd%inum
104      ELSE
105         iadd = 0
106      ENDIF
107
108      IF ( PRESENT( pext ) ) THEN
109         iext = pext%inum
110      ELSE
111         iext = 0
112      ENDIF
113
114      CALL init_obfbdata( fbdata )
115
116      ! Find maximum level
117      ilevel = 0
118      DO jvar = 1, 2
119         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
120      END DO
121
122      SELECT CASE ( TRIM(profdata%cvars(1)) )
123      CASE('POTM')
124
125         clfiletype='profb'
126         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
127            &                 1 + iadd, 1 + iext, .TRUE. )
128         fbdata%cname(1)      = profdata%cvars(1)
129         fbdata%cname(2)      = profdata%cvars(2)
130         fbdata%coblong(1)    = 'Potential temperature'
131         fbdata%coblong(2)    = 'Practical salinity'
132         fbdata%cobunit(1)    = 'Degrees centigrade'
133         fbdata%cobunit(2)    = 'PSU'
134         fbdata%cextname(1)   = 'TEMP'
135         fbdata%cextlong(1)   = 'Insitu temperature'
136         fbdata%cextunit(1)   = 'Degrees centigrade'
137         fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
138         fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
139         fbdata%caddunit(1,1) = 'Degrees centigrade'
140         fbdata%caddunit(1,2) = 'PSU'
141         fbdata%cgrid(:)      = 'T'
142         DO je = 1, iext
143            fbdata%cextname(1+je) = pext%cdname(je)
144            fbdata%cextlong(1+je) = pext%cdlong(je,1)
145            fbdata%cextunit(1+je) = pext%cdunit(je,1)
146         END DO
147         DO ja = 1, iadd
148            fbdata%caddname(1+ja) = padd%cdname(ja)
149            DO jvar = 1, 2
150               fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
151               fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)
152            END DO
153         END DO
154
155      CASE('UVEL')
156
157         clfiletype='velfb'
158         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. )
159         fbdata%cname(1)      = profdata%cvars(1)
160         fbdata%cname(2)      = profdata%cvars(2)
161         fbdata%coblong(1)    = 'Zonal velocity'
162         fbdata%coblong(2)    = 'Meridional velocity'
163         fbdata%cobunit(1)    = 'm/s'
164         fbdata%cobunit(2)    = 'm/s'
165         DO je = 1, iext
166            fbdata%cextname(je) = pext%cdname(je)
167            fbdata%cextlong(je) = pext%cdlong(je,1)
168            fbdata%cextunit(je) = pext%cdunit(je,1)
169         END DO
170         fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
171         fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
172         fbdata%caddunit(1,1) = 'm/s'
173         fbdata%caddunit(1,2) = 'm/s'
174         fbdata%caddname(2)   = 'HxG'
175         fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)'
176         fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)'
177         fbdata%caddunit(2,1) = 'm/s'
178         fbdata%caddunit(2,2) = 'm/s' 
179         fbdata%cgrid(1)      = 'U' 
180         fbdata%cgrid(2)      = 'V'
181         DO ja = 1, iadd
182            fbdata%caddname(2+ja) = padd%cdname(ja)
183            fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
184            fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
185         END DO
186         ALLOCATE( &
187            & zu(profdata%nvprot(1)), &
188            & zv(profdata%nvprot(2))  &
189            & )
190         CALL obs_rotvel( profdata, k2dint, zu, zv )
191
192      END SELECT
193
194      fbdata%caddname(1)   = 'Hx'
195
196      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
197
198      IF(lwp) THEN
199         WRITE(numout,*)
200         WRITE(numout,*)'obs_wri_prof :'
201         WRITE(numout,*)'~~~~~~~~~~~~~'
202         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
203      ENDIF
204
205      ! Transform obs_prof data structure into obfbdata structure
206      fbdata%cdjuldref = '19500101000000'
207      DO jo = 1, profdata%nprof
208         fbdata%plam(jo)      = profdata%rlam(jo)
209         fbdata%pphi(jo)      = profdata%rphi(jo)
210         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
211         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
212         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
213         IF ( profdata%nqc(jo) > 10 ) THEN
214            fbdata%ioqc(jo)    = 4
215            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
216            fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10
217         ELSE
218            fbdata%ioqc(jo)    = profdata%nqc(jo)
219            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
220         ENDIF
221         fbdata%ipqc(jo)      = profdata%ipqc(jo)
222         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
223         fbdata%itqc(jo)      = profdata%itqc(jo)
224         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
225         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
226         fbdata%kindex(jo)    = profdata%npfil(jo)
227         DO jvar = 1, profdata%nvar
228            IF (ln_grid_global) THEN
229               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
230               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
231            ELSE
232               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
233               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
234            ENDIF
235         END DO
236         CALL greg2jul( 0, &
237            &           profdata%nmin(jo), &
238            &           profdata%nhou(jo), &
239            &           profdata%nday(jo), &
240            &           profdata%nmon(jo), &
241            &           profdata%nyea(jo), &
242            &           fbdata%ptim(jo),   &
243            &           krefdate = 19500101 )
244         ! Reform the profiles arrays for output
245         DO jvar = 1, 2
246            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
247               ik = profdata%var(jvar)%nvlidx(jk)
248               IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
249                  fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
250               ELSE IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
251                  IF ( jvar == 1 ) THEN
252                     fbdata%padd(ik,jo,1,jvar) = zu(jk)
253                  ELSE
254                     fbdata%padd(ik,jo,1,jvar) = zv(jk)
255                  ENDIF
256                  fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk)
257               ENDIF
258               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
259               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
260               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
261               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
262               IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN
263                  fbdata%ivlqc(ik,jo,jvar) = 4
264                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
265                  fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10
266               ELSE
267                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
268                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
269               ENDIF
270               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
271               DO ja = 1, iadd
272                  fbdata%padd(ik,jo,1+ja,jvar) = &
273                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
274               END DO
275               DO je = 1, iext
276                  fbdata%pext(ik,jo,1+je) = &
277                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
278               END DO
279               IF ( jvar == 1 ) THEN
280                  fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
281               ENDIF
282            END DO
283         END DO
284      END DO
285
286      IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
287         ! Convert insitu temperature to potential temperature using the model
288         ! salinity if no potential temperature
289         DO jo = 1, fbdata%nobs
290            IF ( fbdata%pphi(jo) < 9999.0 ) THEN
291               DO jk = 1, fbdata%nlev
292                  IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
293                     & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
294                     & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
295                     & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
296                     zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
297                        &              REAL(fbdata%pphi(jo),wp) )
298                     fbdata%pob(jk,jo,1) = potemp( &
299                        &                     REAL(fbdata%padd(jk,jo,1,2), wp),  &
300                        &                     REAL(fbdata%pext(jk,jo,1), wp), &
301                        &                     zpres, 0.0_wp )
302                  ENDIF
303               END DO
304            ENDIF
305         END DO
306      ENDIF
307
308      ! Write the obfbdata structure
309      CALL write_obfbdata( clfname, fbdata )
310
311      ! Output some basic statistics
312      CALL obs_wri_stats( fbdata )
313
314      CALL dealloc_obfbdata( fbdata )
315
316   END SUBROUTINE obs_wri_prof
317
318   SUBROUTINE obs_wri_surf( surfdata, padd, pext )
319      !!-----------------------------------------------------------------------
320      !!
321      !!                     *** ROUTINE obs_wri_surf  ***
322      !!
323      !! ** Purpose : Write surface observation files
324      !!
325      !! ** Method  : NetCDF
326      !!
327      !! ** Action  :
328      !!
329      !!      ! 07-03  (K. Mogensen) Original
330      !!      ! 09-01  (K. Mogensen) New feedback format.
331      !!      ! 15-02  (M. Martin) Combined surface writing routine.
332      !!-----------------------------------------------------------------------
333
334      !! * Modules used
335      IMPLICIT NONE
336
337      !! * Arguments
338      TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data
339      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
340      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info
341
342      !! * Local declarations
343      TYPE(obfbdata) :: fbdata
344      CHARACTER(LEN=40) :: clfname         ! netCDF filename
345      CHARACTER(LEN=6)  :: clfiletype
346      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
347      INTEGER :: jo
348      INTEGER :: ja
349      INTEGER :: je
350      INTEGER :: iadd
351      INTEGER :: iext
352
353      IF ( PRESENT( padd ) ) THEN
354         iadd = padd%inum
355      ELSE
356         iadd = 0
357      ENDIF
358
359      IF ( PRESENT( pext ) ) THEN
360         iext = pext%inum
361      ELSE
362         iext = 0
363      ENDIF
364
365      CALL init_obfbdata( fbdata )
366
367      CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
368         &                 2 + iadd, 1 + iext, .TRUE. )
369
370      SELECT CASE ( TRIM(surfdata%cvars(1)) )
371      CASE('SLA')
372
373         clfiletype = 'slafb'
374         fbdata%cname(1)      = surfdata%cvars(1)
375         fbdata%coblong(1)    = 'Sea level anomaly'
376         fbdata%cobunit(1)    = 'Metres'
377         fbdata%cextname(1)   = 'MDT'
378         fbdata%cextlong(1)   = 'Mean dynamic topography'
379         fbdata%cextunit(1)   = 'Metres'
380         DO je = 1, iext
381            fbdata%cextname(je) = pext%cdname(je)
382            fbdata%cextlong(je) = pext%cdlong(je,1)
383            fbdata%cextunit(je) = pext%cdunit(je,1)
384         END DO
385         fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
386         fbdata%caddunit(1,1) = 'Metres' 
387         fbdata%caddname(2)   = 'SSH'
388         fbdata%caddlong(2,1) = 'Model Sea surface height'
389         fbdata%caddunit(2,1) = 'Metres'
390         fbdata%cgrid(1)      = 'T'
391         DO ja = 1, iadd
392            fbdata%caddname(2+ja) = padd%cdname(ja)
393            fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
394            fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
395         END DO
396
397      CASE('SST')
398
399         clfiletype = 'sstfb'
400         fbdata%cname(1)      = surfdata%cvars(1)
401         fbdata%coblong(1)    = 'Sea surface temperature'
402         fbdata%cobunit(1)    = 'Degree centigrade'
403         DO je = 1, iext
404            fbdata%cextname(je) = pext%cdname(je)
405            fbdata%cextlong(je) = pext%cdlong(je,1)
406            fbdata%cextunit(je) = pext%cdunit(je,1)
407         END DO
408         fbdata%caddlong(1,1) = 'Model interpolated SST'
409         fbdata%caddunit(1,1) = 'Degree centigrade'
410         fbdata%cgrid(1)      = 'T'
411         DO ja = 1, iadd
412            fbdata%caddname(1+ja) = padd%cdname(ja)
413            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
414            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
415         END DO
416
417      CASE('SEAICE')
418
419         clfiletype = 'sicfb'
420         fbdata%cname(1)      = surfdata%cvars(1)
421         fbdata%coblong(1)    = 'Sea ice'
422         fbdata%cobunit(1)    = 'Fraction'
423         DO je = 1, iext
424            fbdata%cextname(je) = pext%cdname(je)
425            fbdata%cextlong(je) = pext%cdlong(je,1)
426            fbdata%cextunit(je) = pext%cdunit(je,1)
427         END DO
428         fbdata%caddlong(1,1) = 'Model interpolated ICE'
429         fbdata%caddunit(1,1) = 'Fraction'
430         fbdata%cgrid(1)      = 'T'
431         DO ja = 1, iadd
432            fbdata%caddname(1+ja) = padd%cdname(ja)
433            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
434            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
435         END DO
436
437      END SELECT
438
439      fbdata%caddname(1)   = 'Hx'
440
441      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
442
443      IF(lwp) THEN
444         WRITE(numout,*)
445         WRITE(numout,*)'obs_wri_surf :'
446         WRITE(numout,*)'~~~~~~~~~~~~~'
447         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
448      ENDIF
449
450      ! Transform obs_prof data structure into obfbdata structure
451      fbdata%cdjuldref = '19500101000000'
452      DO jo = 1, surfdata%nsurf
453         fbdata%plam(jo)      = surfdata%rlam(jo)
454         fbdata%pphi(jo)      = surfdata%rphi(jo)
455         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
456         fbdata%ivqc(jo,:)    = 0
457         fbdata%ivqcf(:,jo,:) = 0
458         IF ( surfdata%nqc(jo) > 10 ) THEN
459            fbdata%ioqc(jo)    = 4
460            fbdata%ioqcf(1,jo) = 0
461            fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10
462         ELSE
463            fbdata%ioqc(jo)    = surfdata%nqc(jo)
464            fbdata%ioqcf(:,jo) = 0
465         ENDIF
466         fbdata%ipqc(jo)      = 0
467         fbdata%ipqcf(:,jo)   = 0
468         fbdata%itqc(jo)      = 0
469         fbdata%itqcf(:,jo)   = 0
470         fbdata%cdwmo(jo)     = surfdata%cwmo(jo)
471         fbdata%kindex(jo)    = surfdata%nsfil(jo)
472         IF (ln_grid_global) THEN
473            fbdata%iobsi(jo,1) = surfdata%mi(jo)
474            fbdata%iobsj(jo,1) = surfdata%mj(jo)
475         ELSE
476            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
477            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
478         ENDIF
479         CALL greg2jul( 0, &
480            &           surfdata%nmin(jo), &
481            &           surfdata%nhou(jo), &
482            &           surfdata%nday(jo), &
483            &           surfdata%nmon(jo), &
484            &           surfdata%nyea(jo), &
485            &           fbdata%ptim(jo),   &
486            &           krefdate = 19500101 )
487         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
488         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
489         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
490         fbdata%pdep(1,jo)     = 0.0
491         fbdata%idqc(1,jo)     = 0
492         fbdata%idqcf(:,1,jo)  = 0
493         IF ( surfdata%nqc(jo) > 10 ) THEN
494            fbdata%ivqc(jo,1)       = 4
495            fbdata%ivlqc(1,jo,1)    = 4
496            fbdata%ivlqcf(1,1,jo,1) = 0
497            fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10
498         ELSE
499            fbdata%ivqc(jo,1)       = surfdata%nqc(jo)
500            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo)
501            fbdata%ivlqcf(:,1,jo,1) = 0
502         ENDIF
503         fbdata%iobsk(1,jo,1)  = 0
504         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
505         DO ja = 1, iadd
506            fbdata%padd(1,jo,2+ja,1) = &
507               & surfdata%rext(jo,padd%ipoint(ja))
508         END DO
509         DO je = 1, iext
510            fbdata%pext(1,jo,1+je) = &
511               & surfdata%rext(jo,pext%ipoint(je))
512         END DO
513      END DO
514
515      ! Write the obfbdata structure
516      CALL write_obfbdata( clfname, fbdata )
517
518      ! Output some basic statistics
519      CALL obs_wri_stats( fbdata )
520
521      CALL dealloc_obfbdata( fbdata )
522
523   END SUBROUTINE obs_wri_surf
524
525   SUBROUTINE obs_wri_stats( fbdata )
526      !!-----------------------------------------------------------------------
527      !!
528      !!                     *** ROUTINE obs_wri_stats  ***
529      !!
530      !! ** Purpose : Output some basic statistics of the data being written out
531      !!
532      !! ** Method  :
533      !!
534      !! ** Action  :
535      !!
536      !!      ! 2014-08  (D. Lea) Initial version
537      !!-----------------------------------------------------------------------
538
539      !! * Arguments
540      TYPE(obfbdata) :: fbdata
541
542      !! * Local declarations
543      INTEGER :: jvar
544      INTEGER :: jo
545      INTEGER :: jk
546      INTEGER :: inumgoodobs
547      INTEGER :: inumgoodobsmpp
548      REAL(wp) :: zsumx
549      REAL(wp) :: zsumx2
550      REAL(wp) :: zomb
551
552      IF (lwp) THEN
553         WRITE(numout,*) ''
554         WRITE(numout,*) 'obs_wri_stats :'
555         WRITE(numout,*) '~~~~~~~~~~~~~~~'
556      ENDIF
557
558      DO jvar = 1, fbdata%nvar
559         zsumx=0.0_wp
560         zsumx2=0.0_wp
561         inumgoodobs=0
562         DO jo = 1, fbdata%nobs
563            DO jk = 1, fbdata%nlev
564               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
565                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
566                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
567
568                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
569                  zsumx=zsumx+zomb
570                  zsumx2=zsumx2+zomb**2
571                  inumgoodobs=inumgoodobs+1
572               ENDIF
573            ENDDO
574         ENDDO
575
576         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
577         CALL mpp_sum(zsumx)
578         CALL mpp_sum(zsumx2)
579
580         IF (lwp) THEN
581            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
582            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
583            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
584            WRITE(numout,*) ''
585         ENDIF
586
587      ENDDO
588
589   END SUBROUTINE obs_wri_stats
590
591END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.