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/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/OBS – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/OBS/obs_write.F90 @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

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