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 NEMO/trunk/src/OCE/OBS – NEMO

source: NEMO/trunk/src/OCE/OBS/obs_write.F90

Last change on this file was 14275, checked in by smasson, 3 years ago

trunk: suppress nproc ( = mpprank = narea-1)

  • Property svn:keywords set to Id
File size: 23.4 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   !! * Substitutions
49#  include "do_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
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      CHARACTER(LEN=12) :: clfmt            ! writing format
92      INTEGER :: idg                        ! number of digits
93      INTEGER :: ilevel
94      INTEGER :: jvar
95      INTEGER :: jo
96      INTEGER :: jk
97      INTEGER :: ik
98      INTEGER :: ja
99      INTEGER :: je
100      INTEGER :: iadd
101      INTEGER :: iext
102      REAL(wp) :: zpres
103
104      IF ( PRESENT( padd ) ) THEN
105         iadd = padd%inum
106      ELSE
107         iadd = 0
108      ENDIF
109
110      IF ( PRESENT( pext ) ) THEN
111         iext = pext%inum
112      ELSE
113         iext = 0
114      ENDIF
115
116      CALL init_obfbdata( fbdata )
117
118      ! Find maximum level
119      ilevel = 0
120      DO jvar = 1, profdata%nvar
121         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
122      END DO
123
124      SELECT CASE ( TRIM(profdata%cvars(1)) )
125      CASE('POTM')
126
127         clfiletype='profb'
128         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
129            &                 1 + iadd, 1 + iext, .TRUE. )
130         fbdata%cname(1)      = profdata%cvars(1)
131         fbdata%cname(2)      = profdata%cvars(2)
132         fbdata%coblong(1)    = 'Potential temperature'
133         fbdata%coblong(2)    = 'Practical salinity'
134         fbdata%cobunit(1)    = 'Degrees centigrade'
135         fbdata%cobunit(2)    = 'PSU'
136         fbdata%cextname(1)   = 'TEMP'
137         fbdata%cextlong(1)   = 'Insitu temperature'
138         fbdata%cextunit(1)   = 'Degrees centigrade'
139         fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
140         fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
141         fbdata%caddunit(1,1) = 'Degrees centigrade'
142         fbdata%caddunit(1,2) = 'PSU'
143         fbdata%cgrid(:)      = 'T'
144         DO je = 1, iext
145            fbdata%cextname(1+je) = pext%cdname(je)
146            fbdata%cextlong(1+je) = pext%cdlong(je,1)
147            fbdata%cextunit(1+je) = pext%cdunit(je,1)
148         END DO
149         DO ja = 1, iadd
150            fbdata%caddname(1+ja) = padd%cdname(ja)
151            DO jvar = 1, 2
152               fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
153               fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)
154            END DO
155         END DO
156
157      CASE('UVEL')
158
159         clfiletype='velfb'
160         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 1, 0, .TRUE. )
161         fbdata%cname(1)      = profdata%cvars(1)
162         fbdata%cname(2)      = profdata%cvars(2)
163         fbdata%coblong(1)    = 'Zonal velocity'
164         fbdata%coblong(2)    = 'Meridional velocity'
165         fbdata%cobunit(1)    = 'm/s'
166         fbdata%cobunit(2)    = 'm/s'
167         DO je = 1, iext
168            fbdata%cextname(je) = pext%cdname(je)
169            fbdata%cextlong(je) = pext%cdlong(je,1)
170            fbdata%cextunit(je) = pext%cdunit(je,1)
171         END DO
172         fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
173         fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
174         fbdata%caddunit(1,1) = 'm/s'
175         fbdata%caddunit(1,2) = 'm/s'
176         fbdata%cgrid(1)      = 'U' 
177         fbdata%cgrid(2)      = 'V'
178         DO ja = 1, iadd
179            fbdata%caddname(1+ja) = padd%cdname(ja)
180            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
181            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
182         END DO
183
184      END SELECT
185     
186      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. &
187         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
188         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
189            &                 1 + iadd, iext, .TRUE. )
190         fbdata%cname(1)      = profdata%cvars(1)
191         fbdata%coblong(1)    = cllongname
192         fbdata%cobunit(1)    = clunits
193         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
194         fbdata%caddunit(1,1) = clunits
195         fbdata%cgrid(:)      = clgrid
196         DO je = 1, iext
197            fbdata%cextname(je) = pext%cdname(je)
198            fbdata%cextlong(je) = pext%cdlong(je,1)
199            fbdata%cextunit(je) = pext%cdunit(je,1)
200         END DO
201         DO ja = 1, iadd
202            fbdata%caddname(1+ja) = padd%cdname(ja)
203            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
204            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
205         END DO
206      ENDIF
207
208      fbdata%caddname(1)   = 'Hx'
209
210      idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )            ! how many digits to we need to write? min=4, max=9
211      WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg   ! '(a,a,ix.x,a)'
212      WRITE(clfname,clfmt) TRIM(clfiletype), '_fdbk_', narea-1, '.nc'
213
214      IF(lwp) THEN
215         WRITE(numout,*)
216         WRITE(numout,*)'obs_wri_prof :'
217         WRITE(numout,*)'~~~~~~~~~~~~~'
218         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
219      ENDIF
220
221      ! Transform obs_prof data structure into obfb data structure
222      fbdata%cdjuldref = '19500101000000'
223      DO jo = 1, profdata%nprof
224         fbdata%plam(jo)      = profdata%rlam(jo)
225         fbdata%pphi(jo)      = profdata%rphi(jo)
226         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
227         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
228         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
229         IF ( profdata%nqc(jo) > 255 ) THEN
230            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2)
231            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
232            fbdata%ioqcf(2,jo) = profdata%nqc(jo)
233         ELSE
234            fbdata%ioqc(jo)    = profdata%nqc(jo)
235            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
236         ENDIF
237         fbdata%ipqc(jo)      = profdata%ipqc(jo)
238         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
239         fbdata%itqc(jo)      = profdata%itqc(jo)
240         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
241         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
242         fbdata%kindex(jo)    = profdata%npfil(jo)
243         DO jvar = 1, profdata%nvar
244            IF (ln_grid_global) THEN
245               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
246               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
247            ELSE
248               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
249               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
250            ENDIF
251         END DO
252         CALL greg2jul( 0, &
253            &           profdata%nmin(jo), &
254            &           profdata%nhou(jo), &
255            &           profdata%nday(jo), &
256            &           profdata%nmon(jo), &
257            &           profdata%nyea(jo), &
258            &           fbdata%ptim(jo),   &
259            &           krefdate = 19500101 )
260         ! Reform the profiles arrays for output
261         DO jvar = 1, profdata%nvar
262            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
263               ik = profdata%var(jvar)%nvlidx(jk)
264               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
265               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
266               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
267               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
268               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
269               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
270                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
271                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
272!$AGRIF_DO_NOT_TREAT
273                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111')
274!$AGRIF_END_DO_NOT_TREAT
275               ELSE
276                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
277                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
278               ENDIF
279               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
280               DO ja = 1, iadd
281                  fbdata%padd(ik,jo,1+ja,jvar) = &
282                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
283               END DO
284               DO je = 1, iext
285                  fbdata%pext(ik,jo,1+je) = &
286                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
287               END DO
288               IF ( ( jvar == 1 ) .AND. &
289                  & ( TRIM(profdata%cvars(1)) == 'POTM' ) ) THEN
290                  fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
291               ENDIF
292            END DO
293         END DO
294      END DO
295
296      IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
297         ! Convert insitu temperature to potential temperature using the model
298         ! salinity if no potential temperature
299         DO jo = 1, fbdata%nobs
300            IF ( fbdata%pphi(jo) < 9999.0 ) THEN
301               DO jk = 1, fbdata%nlev
302                  IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
303                     & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
304                     & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
305                     & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
306                     zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
307                        &              REAL(fbdata%pphi(jo),wp) )
308                     fbdata%pob(jk,jo,1) = potemp( &
309                        &                     REAL(fbdata%padd(jk,jo,1,2), wp),  &
310                        &                     REAL(fbdata%pext(jk,jo,1), wp), &
311                        &                     zpres, 0.0_wp )
312                  ENDIF
313               END DO
314            ENDIF
315         END DO
316      ENDIF
317
318      ! Write the obfbdata structure
319      CALL write_obfbdata( clfname, fbdata )
320
321      ! Output some basic statistics
322      CALL obs_wri_stats( fbdata )
323
324      CALL dealloc_obfbdata( fbdata )
325
326   END SUBROUTINE obs_wri_prof
327
328   SUBROUTINE obs_wri_surf( surfdata, padd, pext )
329      !!-----------------------------------------------------------------------
330      !!
331      !!                     *** ROUTINE obs_wri_surf  ***
332      !!
333      !! ** Purpose : Write surface observation files
334      !!
335      !! ** Method  : NetCDF
336      !!
337      !! ** Action  :
338      !!
339      !!      ! 07-03  (K. Mogensen) Original
340      !!      ! 09-01  (K. Mogensen) New feedback format.
341      !!      ! 15-02  (M. Martin) Combined surface writing routine.
342      !!-----------------------------------------------------------------------
343
344      !! * Modules used
345      IMPLICIT NONE
346
347      !! * Arguments
348      TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data
349      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
350      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info
351
352      !! * Local declarations
353      TYPE(obfbdata) :: fbdata
354      CHARACTER(LEN=40) :: clfname         ! netCDF filename
355      CHARACTER(LEN=10) :: clfiletype
356      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
357      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
358      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
359      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
360      CHARACTER(LEN=12) :: clfmt           ! writing format
361      INTEGER :: idg                       ! number of digits
362      INTEGER :: jo
363      INTEGER :: ja
364      INTEGER :: je
365      INTEGER :: iadd
366      INTEGER :: iext
367
368      IF ( PRESENT( padd ) ) THEN
369         iadd = padd%inum
370      ELSE
371         iadd = 0
372      ENDIF
373
374      IF ( PRESENT( pext ) ) THEN
375         iext = pext%inum
376      ELSE
377         iext = 0
378      ENDIF
379
380      CALL init_obfbdata( fbdata )
381
382      SELECT CASE ( TRIM(surfdata%cvars(1)) )
383      CASE('SLA')
384         
385         ! SLA needs special treatment because of MDT, so is all done here
386         ! Other variables are done more generically
387         ! No climatology for SLA, MDT is our best estimate of that and is already output.
388
389         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
390            &                 2 + iadd, 1 + iext, .TRUE. )
391
392         clfiletype = 'slafb'
393         fbdata%cname(1)      = surfdata%cvars(1)
394         fbdata%coblong(1)    = 'Sea level anomaly'
395         fbdata%cobunit(1)    = 'Metres'
396         fbdata%cextname(1)   = 'MDT'
397         fbdata%cextlong(1)   = 'Mean dynamic topography'
398         fbdata%cextunit(1)   = 'Metres'
399         DO je = 1, iext
400            fbdata%cextname(je) = pext%cdname(je)
401            fbdata%cextlong(je) = pext%cdlong(je,1)
402            fbdata%cextunit(je) = pext%cdunit(je,1)
403         END DO
404         fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
405         fbdata%caddunit(1,1) = 'Metres' 
406         fbdata%caddname(2)   = 'SSH'
407         fbdata%caddlong(2,1) = 'Model Sea surface height'
408         fbdata%caddunit(2,1) = 'Metres'
409         fbdata%cgrid(1)      = 'T'
410         DO ja = 1, iadd
411            fbdata%caddname(2+ja) = padd%cdname(ja)
412            fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
413            fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
414         END DO
415
416      CASE('SST')
417
418         clfiletype = 'sstfb'
419         cllongname = 'Sea surface temperature'
420         clunits    = 'Degree centigrade'
421         clgrid     = 'T'
422         
423      CASE('ICECONC')
424
425         clfiletype = 'sicfb'
426         cllongname = 'Sea ice concentration'
427         clunits    = 'Fraction'
428         clgrid     = 'T'
429
430      CASE('SSS')
431
432         clfiletype = 'sssfb'
433         cllongname = 'Sea surface salinity'
434         clunits    = 'psu'
435         clgrid     = 'T'
436
437      CASE DEFAULT
438
439         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
440
441      END SELECT
442
443      ! SLA needs special treatment because of MDT, so is done above
444      ! Remaining variables treated more generically
445
446      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
447     
448         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
449            &                 1 + iadd, iext, .TRUE. )
450
451         fbdata%cname(1)      = surfdata%cvars(1)
452         fbdata%coblong(1)    = cllongname
453         fbdata%cobunit(1)    = clunits
454         DO je = 1, iext
455            fbdata%cextname(je) = pext%cdname(je)
456            fbdata%cextlong(je) = pext%cdlong(je,1)
457            fbdata%cextunit(je) = pext%cdunit(je,1)
458         END DO       
459         IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
460            fbdata%caddlong(1,1) = 'Model interpolated ICE'
461         ELSE
462            fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1))
463         ENDIF
464         fbdata%caddunit(1,1) = clunits
465         fbdata%cgrid(1)      = clgrid
466         DO ja = 1, iadd
467            fbdata%caddname(1+ja) = padd%cdname(ja)
468            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
469            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
470         END DO
471      ENDIF
472
473      fbdata%caddname(1)   = 'Hx'
474
475      idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )            ! how many digits to we need to write? min=4, max=9
476      WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg   ! '(a,a,ix.x,a)'
477      WRITE(clfname,clfmt) TRIM(clfiletype), '_fdbk_', narea-1, '.nc'
478
479      IF(lwp) THEN
480         WRITE(numout,*)
481         WRITE(numout,*)'obs_wri_surf :'
482         WRITE(numout,*)'~~~~~~~~~~~~~'
483         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
484      ENDIF
485
486      ! Transform surf data structure into obfbdata structure
487      fbdata%cdjuldref = '19500101000000'
488      DO jo = 1, surfdata%nsurf
489         fbdata%plam(jo)      = surfdata%rlam(jo)
490         fbdata%pphi(jo)      = surfdata%rphi(jo)
491         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
492         fbdata%ivqc(jo,:)    = 0
493         fbdata%ivqcf(:,jo,:) = 0
494         IF ( surfdata%nqc(jo) > 255 ) THEN
495            fbdata%ioqc(jo)    = 4
496            fbdata%ioqcf(1,jo) = 0
497!$AGRIF_DO_NOT_TREAT
498            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111')
499!$AGRIF_END_DO_NOT_TREAT
500         ELSE
501            fbdata%ioqc(jo)    = surfdata%nqc(jo)
502            fbdata%ioqcf(:,jo) = 0
503         ENDIF
504         fbdata%ipqc(jo)      = 0
505         fbdata%ipqcf(:,jo)   = 0
506         fbdata%itqc(jo)      = 0
507         fbdata%itqcf(:,jo)   = 0
508         fbdata%cdwmo(jo)     = surfdata%cwmo(jo)
509         fbdata%kindex(jo)    = surfdata%nsfil(jo)
510         IF (ln_grid_global) THEN
511            fbdata%iobsi(jo,1) = surfdata%mi(jo)
512            fbdata%iobsj(jo,1) = surfdata%mj(jo)
513         ELSE
514            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
515            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
516         ENDIF
517         CALL greg2jul( 0, &
518            &           surfdata%nmin(jo), &
519            &           surfdata%nhou(jo), &
520            &           surfdata%nday(jo), &
521            &           surfdata%nmon(jo), &
522            &           surfdata%nyea(jo), &
523            &           fbdata%ptim(jo),   &
524            &           krefdate = 19500101 )
525         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
526         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
527         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
528         fbdata%pdep(1,jo)     = 0.0
529         fbdata%idqc(1,jo)     = 0
530         fbdata%idqcf(:,1,jo)  = 0
531         IF ( surfdata%nqc(jo) > 255 ) THEN
532            fbdata%ivqc(jo,1)       = 4
533            fbdata%ivlqc(1,jo,1)    = 4
534            fbdata%ivlqcf(1,1,jo,1) = 0
535!$AGRIF_DO_NOT_TREAT
536            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000000011111111')
537!$AGRIF_END_DO_NOT_TREAT
538         ELSE
539            fbdata%ivqc(jo,1)       = surfdata%nqc(jo)
540            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo)
541            fbdata%ivlqcf(:,1,jo,1) = 0
542         ENDIF
543         fbdata%iobsk(1,jo,1)  = 0
544         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
545         DO ja = 1, iadd
546            fbdata%padd(1,jo,2+ja,1) = &
547               & surfdata%rext(jo,padd%ipoint(ja))
548         END DO
549         DO je = 1, iext
550            fbdata%pext(1,jo,1+je) = &
551               & surfdata%rext(jo,pext%ipoint(je))
552         END DO
553      END DO
554
555      ! Write the obfbdata structure
556      CALL write_obfbdata( clfname, fbdata )
557
558      ! Output some basic statistics
559      CALL obs_wri_stats( fbdata )
560
561      CALL dealloc_obfbdata( fbdata )
562
563   END SUBROUTINE obs_wri_surf
564
565   SUBROUTINE obs_wri_stats( fbdata )
566      !!-----------------------------------------------------------------------
567      !!
568      !!                     *** ROUTINE obs_wri_stats  ***
569      !!
570      !! ** Purpose : Output some basic statistics of the data being written out
571      !!
572      !! ** Method  :
573      !!
574      !! ** Action  :
575      !!
576      !!      ! 2014-08  (D. Lea) Initial version
577      !!-----------------------------------------------------------------------
578
579      !! * Arguments
580      TYPE(obfbdata) :: fbdata
581
582      !! * Local declarations
583      INTEGER :: jvar
584      INTEGER :: jo
585      INTEGER :: jk
586      INTEGER :: inumgoodobs
587      INTEGER :: inumgoodobsmpp
588      REAL(wp) :: zsumx
589      REAL(wp) :: zsumx2
590      REAL(wp) :: zomb
591     
592
593      IF (lwp) THEN
594         WRITE(numout,*) ''
595         WRITE(numout,*) 'obs_wri_stats :'
596         WRITE(numout,*) '~~~~~~~~~~~~~~~'
597      ENDIF
598
599      DO jvar = 1, fbdata%nvar
600         zsumx=0.0_wp
601         zsumx2=0.0_wp
602         inumgoodobs=0
603         DO jo = 1, fbdata%nobs
604            DO jk = 1, fbdata%nlev
605               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
606                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
607                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
608
609                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
610                  zsumx=zsumx+zomb
611                  zsumx2=zsumx2+zomb**2
612                  inumgoodobs=inumgoodobs+1
613               ENDIF
614            ENDDO
615         ENDDO
616
617         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
618         CALL mpp_sum('obs_write', zsumx)
619         CALL mpp_sum('obs_write', zsumx2)
620
621         IF (lwp) THEN
622            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
623            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
624            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
625            WRITE(numout,*) ''
626         ENDIF
627
628      ENDDO
629
630   END SUBROUTINE obs_wri_stats
631
632END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.