source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obs_write.F90 @ 13135

Last change on this file since 13135 was 13135, checked in by orioltp, 3 months ago

dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: merge with trunk@13134, see #2364

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