source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 35.2 KB
Line 
1MODULE obs_write
2   !!======================================================================
3   !!                       ***  MODULE obs_write   ***
4   !! Observation diagnosticss: Write observation related diagnostics
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_wri_p3d   : Write profile observation diagnostics in NetCDF format
9   !!   obs_wri_sla   : Write SLA observation related diagnostics
10   !!   obs_wri_sst   : Write SST observation related diagnostics
11   !!   obs_wri_seaice: Write seaice observation related diagnostics
12   !!   obs_wri_vel   : Write velocity observation diagnostics in NetCDF format
13   !!   obs_wri_stats : Print basic statistics on the data being written out
14   !!----------------------------------------------------------------------
15
16   !! * Modules used
17   USE par_kind, ONLY : &   ! Precision variables
18      & wp
19   USE in_out_manager       ! I/O manager
20   USE dom_oce              ! Ocean space and time domain variables
21   USE obs_types            ! Observation type integer to character translation
22   USE julian, ONLY : &         ! Julian date routines
23      & greg2jul
24   USE obs_utils, ONLY : &  ! Observation operator utility functions
25      & chkerr
26   USE obs_profiles_def     ! Type definitions for profiles
27   USE obs_surf_def         ! Type defintions for surface observations
28   USE obs_fbm              ! Observation feedback I/O
29   USE obs_grid             ! Grid tools
30   USE obs_conv             ! Conversion between units
31   USE obs_const
32   USE obs_sla_types
33   USE obs_rot_vel          ! Rotation of velocities
34   USE obs_mpp              ! MPP support routines for observation diagnostics
35   USE lib_mpp        ! MPP routines
36
37   IMPLICIT NONE
38
39   !! * Routine accessibility
40   PRIVATE
41   PUBLIC obs_wri_p3d, &    ! Write profile observation related diagnostics
42      &   obs_wri_sla, &    ! Write SLA observation related diagnostics
43      &   obs_wri_sst, &    ! Write SST observation related diagnostics
44      &   obs_wri_sss, &    ! Write SSS observation related diagnostics
45      &   obs_wri_seaice, & ! Write seaice observation related diagnostics
46      &   obs_wri_vel, &    ! Write velocity observation related diagnostics
47      &   obswriinfo
48   
49   TYPE obswriinfo
50      INTEGER :: inum
51      INTEGER, POINTER, DIMENSION(:) :: ipoint
52      CHARACTER(len=ilenname), POINTER, DIMENSION(:) :: cdname
53      CHARACTER(len=ilenlong), POINTER, DIMENSION(:,:) :: cdlong
54      CHARACTER(len=ilenunit), POINTER, DIMENSION(:,:) :: cdunit
55   END TYPE obswriinfo
56
57   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
59   !! $Id$
60   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62
63CONTAINS
64
65   SUBROUTINE obs_wri_p3d( cprefix, profdata, padd, pext )
66      !!-----------------------------------------------------------------------
67      !!
68      !!                     *** ROUTINE obs_wri_p3d  ***
69      !!
70      !! ** Purpose : Write temperature and salinity (profile) observation
71      !!              related diagnostics
72      !!
73      !! ** Method  : NetCDF
74      !!
75      !! ** Action  :
76      !!
77      !! History :
78      !!      ! 06-04  (A. Vidard) Original
79      !!      ! 06-04  (A. Vidard) Reformatted
80      !!      ! 06-10  (A. Weaver) Cleanup
81      !!      ! 07-01  (K. Mogensen) Use profile data types
82      !!      ! 07-03  (K. Mogensen) General handling of profiles
83      !!      ! 09-01  (K. Mogensen) New feedback format
84      !!-----------------------------------------------------------------------
85
86      !! * Modules used
87
88      !! * Arguments
89      CHARACTER(LEN=*), INTENT(IN) :: cprefix        ! Prefix for output files
90      TYPE(obs_prof), INTENT(INOUT) :: profdata      ! Full set of profile data
91      TYPE(obswriinfo), OPTIONAL :: padd             ! Additional info for each variable
92      TYPE(obswriinfo), OPTIONAL :: pext             ! Extra info
93     
94      !! * Local declarations
95      TYPE(obfbdata) :: fbdata
96      CHARACTER(LEN=40) :: cfname
97      INTEGER :: ilevel
98      INTEGER :: jvar
99      INTEGER :: jo
100      INTEGER :: jk
101      INTEGER :: ik
102      INTEGER :: ja
103      INTEGER :: je
104      REAL(wp) :: zpres
105      INTEGER :: nadd
106      INTEGER :: next
107
108      IF ( PRESENT( padd ) ) THEN
109         nadd = padd%inum
110      ELSE
111         nadd = 0
112      ENDIF
113
114      IF ( PRESENT( pext ) ) THEN
115         next = pext%inum
116      ELSE
117         next = 0
118      ENDIF
119     
120      CALL init_obfbdata( fbdata )
121
122      ! Find maximum level
123      ilevel = 0
124      DO jvar = 1, 2
125         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
126      END DO
127      CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
128         &                 1 + nadd, 1 + next, .TRUE. )
129
130      fbdata%cname(1)      = 'POTM'
131      fbdata%cname(2)      = 'PSAL'
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      DO je = 1, next
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      fbdata%caddname(1)   = 'Hx'
145      fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
146      fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
147      fbdata%caddunit(1,1) = 'Degrees centigrade'
148      fbdata%caddunit(1,2) = 'PSU'
149      fbdata%cgrid(:)      = 'T'
150      DO ja = 1, nadd
151         fbdata%caddname(1+ja) = padd%cdname(ja)
152         DO jvar = 1, 2
153            fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
154            fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)
155         END DO
156      END DO
157         
158      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
159
160      IF(lwp) THEN
161         WRITE(numout,*)
162         WRITE(numout,*)'obs_wri_p3d :'
163         WRITE(numout,*)'~~~~~~~~~~~~~'
164         WRITE(numout,*)'Writing profile feedback file : ',TRIM(cfname)
165      ENDIF
166
167      ! Transform obs_prof data structure into obfbdata structure
168      fbdata%cdjuldref = '19500101000000'
169      DO jo = 1, profdata%nprof
170         fbdata%plam(jo)      = profdata%rlam(jo)
171         fbdata%pphi(jo)      = profdata%rphi(jo)
172         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
173         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
174         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
175         IF ( profdata%nqc(jo) > 10 ) THEN
176            fbdata%ioqc(jo)    = 4
177            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
178            fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10
179         ELSE
180            fbdata%ioqc(jo)    = profdata%nqc(jo)
181            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
182         ENDIF
183         fbdata%ipqc(jo)      = profdata%ipqc(jo)
184         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
185         fbdata%itqc(jo)      = profdata%itqc(jo)
186         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
187         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
188         fbdata%kindex(jo)    = profdata%npfil(jo)
189         DO jvar = 1, profdata%nvar
190            IF (ln_grid_global) THEN
191               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
192               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
193            ELSE
194               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
195               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
196            ENDIF
197         END DO
198         CALL greg2jul( 0, &
199            &           profdata%nmin(jo), &
200            &           profdata%nhou(jo), &
201            &           profdata%nday(jo), &
202            &           profdata%nmon(jo), &
203            &           profdata%nyea(jo), &
204            &           fbdata%ptim(jo),   &
205            &           krefdate = 19500101 )
206         ! Reform the profiles arrays for output
207         DO jvar = 1, 2
208            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
209               ik = profdata%var(jvar)%nvlidx(jk)
210               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
211               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
212               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
213               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
214               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
215               IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN
216                  fbdata%ivlqc(ik,jo,jvar) = 4
217                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
218                  fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10
219               ELSE
220                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
221                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
222               ENDIF
223               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
224               DO ja = 1, nadd
225                  fbdata%padd(ik,jo,1+ja,jvar) = &
226                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
227               END DO
228               DO je = 1, next
229                  fbdata%pext(ik,jo,1+je) = &
230                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
231               END DO
232               IF ( jvar == 1 ) THEN
233                  fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
234               ENDIF
235            END DO
236         END DO
237      END DO
238
239      ! Convert insitu temperature to potential temperature using the model
240      ! salinity if no potential temperature
241      DO jo = 1, fbdata%nobs
242         IF ( fbdata%pphi(jo) < 9999.0 ) THEN
243            DO jk = 1, fbdata%nlev
244               IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
245                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
246                  & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
247                  & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
248                  zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
249                     &              REAL(fbdata%pphi(jo),wp) )
250                  fbdata%pob(jk,jo,1) = potemp( &
251                     &                     REAL(fbdata%padd(jk,jo,1,2), wp),  &
252                     &                     REAL(fbdata%pext(jk,jo,1), wp), &
253                     &                     zpres, 0.0_wp )
254               ENDIF
255            END DO
256         ENDIF
257      END DO
258     
259      ! Write the obfbdata structure
260      CALL write_obfbdata( cfname, fbdata )
261
262      ! Output some basic statistics
263      CALL obs_wri_stats( fbdata )
264
265      CALL dealloc_obfbdata( fbdata )
266     
267   END SUBROUTINE obs_wri_p3d
268
269   SUBROUTINE obs_wri_sla( cprefix, sladata, padd, pext )
270      !!-----------------------------------------------------------------------
271      !!
272      !!                     *** ROUTINE obs_wri_sla  ***
273      !!
274      !! ** Purpose : Write SLA observation diagnostics
275      !!              related
276      !!
277      !! ** Method  : NetCDF
278      !!
279      !! ** Action  :
280      !!
281      !!      ! 07-03  (K. Mogensen) Original
282      !!      ! 09-01  (K. Mogensen) New feedback format.
283      !!-----------------------------------------------------------------------
284
285      !! * Modules used
286      IMPLICIT NONE
287
288      !! * Arguments
289      CHARACTER(LEN=*), INTENT(IN) :: cprefix          ! Prefix for output files
290      TYPE(obs_surf), INTENT(INOUT) :: sladata         ! Full set of SLAa
291      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
292      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info
293
294      !! * Local declarations
295      TYPE(obfbdata) :: fbdata
296      CHARACTER(LEN=40) :: cfname         ! netCDF filename
297      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_sla'
298      INTEGER :: jo
299      INTEGER :: ja
300      INTEGER :: je
301      INTEGER :: nadd
302      INTEGER :: next
303
304      IF ( PRESENT( padd ) ) THEN
305         nadd = padd%inum
306      ELSE
307         nadd = 0
308      ENDIF
309
310      IF ( PRESENT( pext ) ) THEN
311         next = pext%inum
312      ELSE
313         next = 0
314      ENDIF
315
316      CALL init_obfbdata( fbdata )
317
318      CALL alloc_obfbdata( fbdata, 1, sladata%nsurf, 1, &
319         &                 2 + nadd, 1 + next, .TRUE. )
320
321      fbdata%cname(1)      = 'SLA'
322      fbdata%coblong(1)    = 'Sea level anomaly'
323      fbdata%cobunit(1)    = 'Metres'
324      fbdata%cextname(1)   = 'MDT'
325      fbdata%cextlong(1)   = 'Mean dynamic topography'
326      fbdata%cextunit(1)   = 'Metres'
327      DO je = 1, next
328         fbdata%cextname(1+je) = pext%cdname(je)
329         fbdata%cextlong(1+je) = pext%cdlong(je,1)
330         fbdata%cextunit(1+je) = pext%cdunit(je,1)
331      END DO
332      fbdata%caddname(1)   = 'Hx'
333      fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
334      fbdata%caddunit(1,1) = 'Metres' 
335      fbdata%caddname(2)   = 'SSH'
336      fbdata%caddlong(2,1) = 'Model Sea surface height'
337      fbdata%caddunit(2,1) = 'Metres'
338      fbdata%cgrid(1)      = 'T'
339      DO ja = 1, nadd
340         fbdata%caddname(2+ja) = padd%cdname(ja)
341         fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
342         fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
343      END DO
344
345      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
346
347      IF(lwp) THEN
348         WRITE(numout,*)
349         WRITE(numout,*)'obs_wri_sla :'
350         WRITE(numout,*)'~~~~~~~~~~~~~'
351         WRITE(numout,*)'Writing SLA feedback file : ',TRIM(cfname)
352      ENDIF
353
354      ! Transform obs_prof data structure into obfbdata structure
355      fbdata%cdjuldref = '19500101000000'
356      DO jo = 1, sladata%nsurf
357         fbdata%plam(jo)      = sladata%rlam(jo)
358         fbdata%pphi(jo)      = sladata%rphi(jo)
359         WRITE(fbdata%cdtyp(jo),'(I4)') sladata%ntyp(jo)
360         fbdata%ivqc(jo,:)    = 0
361         fbdata%ivqcf(:,jo,:) = 0
362         IF ( sladata%nqc(jo) > 10 ) THEN
363            fbdata%ioqc(jo)    = 4
364            fbdata%ioqcf(1,jo) = 0
365            fbdata%ioqcf(2,jo) = sladata%nqc(jo) - 10
366         ELSE
367            fbdata%ioqc(jo)    = sladata%nqc(jo)
368            fbdata%ioqcf(:,jo) = 0
369         ENDIF
370         fbdata%ipqc(jo)      = 0
371         fbdata%ipqcf(:,jo)   = 0
372         fbdata%itqc(jo)      = 0
373         fbdata%itqcf(:,jo)   = 0
374         fbdata%cdwmo(jo)     = sladata%cwmo(jo)
375         fbdata%kindex(jo)    = sladata%nsfil(jo)
376         IF (ln_grid_global) THEN
377            fbdata%iobsi(jo,1) = sladata%mi(jo)
378            fbdata%iobsj(jo,1) = sladata%mj(jo)
379         ELSE
380            fbdata%iobsi(jo,1) = mig(sladata%mi(jo))
381            fbdata%iobsj(jo,1) = mjg(sladata%mj(jo))
382         ENDIF
383         CALL greg2jul( 0, &
384            &           sladata%nmin(jo), &
385            &           sladata%nhou(jo), &
386            &           sladata%nday(jo), &
387            &           sladata%nmon(jo), &
388            &           sladata%nyea(jo), &
389            &           fbdata%ptim(jo),   &
390            &           krefdate = 19500101 )
391         fbdata%padd(1,jo,1,1) = sladata%rmod(jo,1)
392         fbdata%padd(1,jo,2,1) = sladata%rext(jo,1)
393         fbdata%pob(1,jo,1)    = sladata%robs(jo,1) 
394         fbdata%pdep(1,jo)     = 0.0
395         fbdata%idqc(1,jo)     = 0
396         fbdata%idqcf(:,1,jo)  = 0
397         IF ( sladata%nqc(jo) > 10 ) THEN
398            fbdata%ivqc(jo,1)       = 4
399            fbdata%ivlqc(1,jo,1)    = 4
400            fbdata%ivlqcf(1,1,jo,1) = 0
401            fbdata%ivlqcf(2,1,jo,1) = sladata%nqc(jo) - 10
402         ELSE
403            fbdata%ivqc(jo,1)       = sladata%nqc(jo)
404            fbdata%ivlqc(1,jo,1)    = sladata%nqc(jo)
405            fbdata%ivlqcf(:,1,jo,1) = 0
406         ENDIF
407         fbdata%iobsk(1,jo,1)  = 0
408         fbdata%pext(1,jo,1) = sladata%rext(jo,2)
409         DO ja = 1, nadd
410            fbdata%padd(1,jo,2+ja,1) = &
411               & sladata%rext(jo,padd%ipoint(ja))
412         END DO
413         DO je = 1, next
414            fbdata%pext(1,jo,1+je) = &
415               & sladata%rext(jo,pext%ipoint(je))
416         END DO
417      END DO
418
419      ! Write the obfbdata structure
420      CALL write_obfbdata( cfname, fbdata )
421
422      ! Output some basic statistics
423      CALL obs_wri_stats( fbdata )
424
425      CALL dealloc_obfbdata( fbdata )
426
427   END SUBROUTINE obs_wri_sla
428
429   SUBROUTINE obs_wri_sst( cprefix, sstdata, padd, pext )
430      !!-----------------------------------------------------------------------
431      !!
432      !!                     *** ROUTINE obs_wri_sst  ***
433      !!
434      !! ** Purpose : Write SST observation diagnostics
435      !!              related
436      !!
437      !! ** Method  : NetCDF
438      !!
439      !! ** Action  :
440      !!
441      !!      ! 07-07  (S. Ricci) Original
442      !!      ! 09-01  (K. Mogensen) New feedback format.
443      !!-----------------------------------------------------------------------
444
445      !! * Modules used
446      IMPLICIT NONE
447
448      !! * Arguments
449      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files
450      TYPE(obs_surf), INTENT(INOUT) :: sstdata      ! Full set of SST
451      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable
452      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info
453
454      !! * Local declarations
455      TYPE(obfbdata) :: fbdata
456      CHARACTER(LEN=40) ::  cfname             ! netCDF filename
457      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_sst'
458      INTEGER :: jo
459      INTEGER :: ja
460      INTEGER :: je
461      INTEGER :: nadd
462      INTEGER :: next
463
464      IF ( PRESENT( padd ) ) THEN
465         nadd = padd%inum
466      ELSE
467         nadd = 0
468      ENDIF
469
470      IF ( PRESENT( pext ) ) THEN
471         next = pext%inum
472      ELSE
473         next = 0
474      ENDIF
475
476      CALL init_obfbdata( fbdata )
477
478      CALL alloc_obfbdata( fbdata, 1, sstdata%nsurf, 1, &
479         &                 1 + nadd, next, .TRUE. )
480
481      fbdata%cname(1)      = 'SST'
482      fbdata%coblong(1)    = 'Sea surface temperature'
483      fbdata%cobunit(1)    = 'Degree centigrade'
484      DO je = 1, next
485         fbdata%cextname(je) = pext%cdname(je)
486         fbdata%cextlong(je) = pext%cdlong(je,1)
487         fbdata%cextunit(je) = pext%cdunit(je,1)
488      END DO
489      fbdata%caddname(1)   = 'Hx'
490      fbdata%caddlong(1,1) = 'Model interpolated SST'
491      fbdata%caddunit(1,1) = 'Degree centigrade'
492      fbdata%cgrid(1)      = 'T'
493      DO ja = 1, nadd
494         fbdata%caddname(1+ja) = padd%cdname(ja)
495         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
496         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
497      END DO
498
499      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
500
501      IF(lwp) THEN
502         WRITE(numout,*)
503         WRITE(numout,*)'obs_wri_sst :'
504         WRITE(numout,*)'~~~~~~~~~~~~~'
505         WRITE(numout,*)'Writing SST feedback file : ',TRIM(cfname)
506      ENDIF
507
508      ! Transform obs_prof data structure into obfbdata structure
509      fbdata%cdjuldref = '19500101000000'
510      DO jo = 1, sstdata%nsurf
511         fbdata%plam(jo)      = sstdata%rlam(jo)
512         fbdata%pphi(jo)      = sstdata%rphi(jo)
513         WRITE(fbdata%cdtyp(jo),'(I4)') sstdata%ntyp(jo)
514         fbdata%ivqc(jo,:)    = 0
515         fbdata%ivqcf(:,jo,:) = 0
516         IF ( sstdata%nqc(jo) > 10 ) THEN
517            fbdata%ioqc(jo)    = 4
518            fbdata%ioqcf(1,jo) = 0
519            fbdata%ioqcf(2,jo) = sstdata%nqc(jo) - 10
520         ELSE
521            fbdata%ioqc(jo)    = MAX(sstdata%nqc(jo),1)
522            fbdata%ioqcf(:,jo) = 0
523         ENDIF
524         fbdata%ipqc(jo)      = 0
525         fbdata%ipqcf(:,jo)   = 0
526         fbdata%itqc(jo)      = 0
527         fbdata%itqcf(:,jo)   = 0
528         fbdata%cdwmo(jo)     = ''
529         fbdata%kindex(jo)    = sstdata%nsfil(jo)
530         IF (ln_grid_global) THEN
531            fbdata%iobsi(jo,1) = sstdata%mi(jo)
532            fbdata%iobsj(jo,1) = sstdata%mj(jo)
533         ELSE
534            fbdata%iobsi(jo,1) = mig(sstdata%mi(jo))
535            fbdata%iobsj(jo,1) = mjg(sstdata%mj(jo))
536         ENDIF
537         CALL greg2jul( 0, &
538            &           sstdata%nmin(jo), &
539            &           sstdata%nhou(jo), &
540            &           sstdata%nday(jo), &
541            &           sstdata%nmon(jo), &
542            &           sstdata%nyea(jo), &
543            &           fbdata%ptim(jo),   &
544            &           krefdate = 19500101 )
545         fbdata%padd(1,jo,1,1) = sstdata%rmod(jo,1)
546         fbdata%pob(1,jo,1)    = sstdata%robs(jo,1)
547         fbdata%pdep(1,jo)     = 0.0
548         fbdata%idqc(1,jo)     = 0
549         fbdata%idqcf(:,1,jo)  = 0
550         IF ( sstdata%nqc(jo) > 10 ) THEN
551            fbdata%ivqc(jo,1)       = 4
552            fbdata%ivlqc(1,jo,1)    = 4
553            fbdata%ivlqcf(1,1,jo,1) = 0
554            fbdata%ivlqcf(2,1,jo,1) = sstdata%nqc(jo) - 10
555         ELSE
556            fbdata%ivqc(jo,1)       = MAX(sstdata%nqc(jo),1)
557            fbdata%ivlqc(1,jo,1)    = MAX(sstdata%nqc(jo),1)
558            fbdata%ivlqcf(:,1,jo,1) = 0
559         ENDIF
560         fbdata%iobsk(1,jo,1)  = 0
561         DO ja = 1, nadd
562            fbdata%padd(1,jo,1+ja,1) = &
563               & sstdata%rext(jo,padd%ipoint(ja))
564         END DO
565         DO je = 1, next
566            fbdata%pext(1,jo,je) = &
567               & sstdata%rext(jo,pext%ipoint(je))
568         END DO
569
570      END DO
571
572      ! Write the obfbdata structure
573
574      CALL write_obfbdata( cfname, fbdata )
575
576      ! Output some basic statistics
577      CALL obs_wri_stats( fbdata )
578
579      CALL dealloc_obfbdata( fbdata )
580
581   END SUBROUTINE obs_wri_sst
582
583   SUBROUTINE obs_wri_sss
584   END SUBROUTINE obs_wri_sss
585
586   SUBROUTINE obs_wri_seaice( cprefix, seaicedata, padd, pext )
587      !!-----------------------------------------------------------------------
588      !!
589      !!                     *** ROUTINE obs_wri_seaice  ***
590      !!
591      !! ** Purpose : Write sea ice observation diagnostics
592      !!              related
593      !!
594      !! ** Method  : NetCDF
595      !!
596      !! ** Action  :
597      !!
598      !!      ! 07-07  (S. Ricci) Original
599      !!      ! 09-01  (K. Mogensen) New feedback format.
600      !!-----------------------------------------------------------------------
601
602      !! * Modules used
603      IMPLICIT NONE
604
605      !! * Arguments
606      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files
607      TYPE(obs_surf), INTENT(INOUT) :: seaicedata   ! Full set of sea ice
608      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable
609      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info
610
611      !! * Local declarations
612      TYPE(obfbdata) :: fbdata
613      CHARACTER(LEN=40) :: cfname             ! netCDF filename
614      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_seaice'
615      INTEGER :: jo
616      INTEGER :: ja
617      INTEGER :: je
618      INTEGER :: nadd
619      INTEGER :: next
620
621      IF ( PRESENT( padd ) ) THEN
622         nadd = padd%inum
623      ELSE
624         nadd = 0
625      ENDIF
626
627      IF ( PRESENT( pext ) ) THEN
628         next = pext%inum
629      ELSE
630         next = 0
631      ENDIF
632
633      CALL init_obfbdata( fbdata )
634
635      CALL alloc_obfbdata( fbdata, 1, seaicedata%nsurf, 1, 1, 0, .TRUE. )
636
637      fbdata%cname(1)      = 'SEAICE'
638      fbdata%coblong(1)    = 'Sea ice'
639      fbdata%cobunit(1)    = 'Fraction'
640      DO je = 1, next
641         fbdata%cextname(je) = pext%cdname(je)
642         fbdata%cextlong(je) = pext%cdlong(je,1)
643         fbdata%cextunit(je) = pext%cdunit(je,1)
644      END DO
645      fbdata%caddname(1)   = 'Hx'
646      fbdata%caddlong(1,1) = 'Model interpolated ICE'
647      fbdata%caddunit(1,1) = 'Fraction'
648      fbdata%cgrid(1)      = 'T'
649      DO ja = 1, nadd
650         fbdata%caddname(1+ja) = padd%cdname(ja)
651         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
652         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
653      END DO
654
655      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
656
657      IF(lwp) THEN
658         WRITE(numout,*)
659         WRITE(numout,*)'obs_wri_seaice :'
660         WRITE(numout,*)'~~~~~~~~~~~~~~~~'
661         WRITE(numout,*)'Writing SEAICE feedback file : ',TRIM(cfname)
662      ENDIF
663
664      ! Transform obs_prof data structure into obfbdata structure
665      fbdata%cdjuldref = '19500101000000'
666      DO jo = 1, seaicedata%nsurf
667         fbdata%plam(jo)      = seaicedata%rlam(jo)
668         fbdata%pphi(jo)      = seaicedata%rphi(jo)
669         WRITE(fbdata%cdtyp(jo),'(I4)') seaicedata%ntyp(jo)
670         fbdata%ivqc(jo,:)    = 0
671         fbdata%ivqcf(:,jo,:) = 0
672         IF ( seaicedata%nqc(jo) > 10 ) THEN
673            fbdata%ioqc(jo)    = 4
674            fbdata%ioqcf(1,jo) = 0
675            fbdata%ioqcf(2,jo) = seaicedata%nqc(jo) - 10
676         ELSE
677            fbdata%ioqc(jo)    = MAX(seaicedata%nqc(jo),1)
678            fbdata%ioqcf(:,jo) = 0
679         ENDIF
680         fbdata%ipqc(jo)      = 0
681         fbdata%ipqcf(:,jo)   = 0
682         fbdata%itqc(jo)      = 0
683         fbdata%itqcf(:,jo)   = 0
684         fbdata%cdwmo(jo)     = ''
685         fbdata%kindex(jo)    = seaicedata%nsfil(jo)
686         IF (ln_grid_global) THEN
687            fbdata%iobsi(jo,1) = seaicedata%mi(jo)
688            fbdata%iobsj(jo,1) = seaicedata%mj(jo)
689         ELSE
690            fbdata%iobsi(jo,1) = mig(seaicedata%mi(jo))
691            fbdata%iobsj(jo,1) = mjg(seaicedata%mj(jo))
692         ENDIF
693         CALL greg2jul( 0, &
694            &           seaicedata%nmin(jo), &
695            &           seaicedata%nhou(jo), &
696            &           seaicedata%nday(jo), &
697            &           seaicedata%nmon(jo), &
698            &           seaicedata%nyea(jo), &
699            &           fbdata%ptim(jo),   &
700            &           krefdate = 19500101 )
701         fbdata%padd(1,jo,1,1) = seaicedata%rmod(jo,1)
702         fbdata%pob(1,jo,1)    = seaicedata%robs(jo,1)
703         fbdata%pdep(1,jo)     = 0.0
704         fbdata%idqc(1,jo)     = 0
705         fbdata%idqcf(:,1,jo)  = 0
706         IF ( seaicedata%nqc(jo) > 10 ) THEN
707            fbdata%ivlqc(1,jo,1) = 4
708            fbdata%ivlqcf(1,1,jo,1) = 0
709            fbdata%ivlqcf(2,1,jo,1) = seaicedata%nqc(jo) - 10
710         ELSE
711            fbdata%ivlqc(1,jo,1) = MAX(seaicedata%nqc(jo),1)
712            fbdata%ivlqcf(:,1,jo,1) = 0
713         ENDIF
714         fbdata%iobsk(1,jo,1)  = 0
715         DO ja = 1, nadd
716            fbdata%padd(1,jo,1+ja,1) = &
717               & seaicedata%rext(jo,padd%ipoint(ja))
718         END DO
719         DO je = 1, next
720            fbdata%pext(1,jo,je) = &
721               & seaicedata%rext(jo,pext%ipoint(je))
722         END DO
723
724      END DO
725
726      ! Write the obfbdata structure
727      CALL write_obfbdata( cfname, fbdata )
728
729      ! Output some basic statistics
730      CALL obs_wri_stats( fbdata )
731
732      CALL dealloc_obfbdata( fbdata )
733
734   END SUBROUTINE obs_wri_seaice
735
736   SUBROUTINE obs_wri_vel( cprefix, profdata, k2dint, padd, pext )
737      !!-----------------------------------------------------------------------
738      !!
739      !!                     *** ROUTINE obs_wri_vel  ***
740      !!
741      !! ** Purpose : Write current (profile) observation
742      !!              related diagnostics
743      !!
744      !! ** Method  : NetCDF
745      !!
746      !! ** Action  :
747      !!
748      !! History :
749      !!      ! 09-01  (K. Mogensen) New feedback format routine
750      !!-----------------------------------------------------------------------
751
752      !! * Modules used
753
754      !! * Arguments
755      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files
756      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Full set of profile data
757      INTEGER, INTENT(IN) :: k2dint                 ! Horizontal interpolation method
758      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable
759      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info
760
761      !! * Local declarations
762      TYPE(obfbdata) :: fbdata
763      CHARACTER(LEN=40) :: cfname
764      INTEGER :: ilevel
765      INTEGER :: jvar
766      INTEGER :: jk
767      INTEGER :: ik
768      INTEGER :: jo
769      INTEGER :: ja
770      INTEGER :: je
771      INTEGER :: nadd
772      INTEGER :: next
773      REAL(wp) :: zpres
774      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
775         & zu, &
776         & zv
777
778      IF ( PRESENT( padd ) ) THEN
779         nadd = padd%inum
780      ELSE
781         nadd = 0
782      ENDIF
783
784      IF ( PRESENT( pext ) ) THEN
785         next = pext%inum
786      ELSE
787         next = 0
788      ENDIF
789
790      CALL init_obfbdata( fbdata )
791
792      ! Find maximum level
793      ilevel = 0
794      DO jvar = 1, 2
795         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
796      END DO
797      CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. )
798
799      fbdata%cname(1)      = 'UVEL'
800      fbdata%cname(2)      = 'VVEL'
801      fbdata%coblong(1)    = 'Zonal velocity'
802      fbdata%coblong(2)    = 'Meridional velocity'
803      fbdata%cobunit(1)    = 'm/s'
804      fbdata%cobunit(2)    = 'm/s'
805      DO je = 1, next
806         fbdata%cextname(je) = pext%cdname(je)
807         fbdata%cextlong(je) = pext%cdlong(je,1)
808         fbdata%cextunit(je) = pext%cdunit(je,1)
809      END DO
810      fbdata%caddname(1)   = 'Hx'
811      fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
812      fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
813      fbdata%caddunit(1,1) = 'm/s'
814      fbdata%caddunit(1,2) = 'm/s'
815      fbdata%caddname(2)   = 'HxG'
816      fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)'
817      fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)'
818      fbdata%caddunit(2,1) = 'm/s'
819      fbdata%caddunit(2,2) = 'm/s' 
820      fbdata%cgrid(1)      = 'U' 
821      fbdata%cgrid(2)      = 'V'
822      DO ja = 1, nadd
823         fbdata%caddname(2+ja) = padd%cdname(ja)
824         fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
825         fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
826      END DO
827
828      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
829
830      IF(lwp) THEN
831         WRITE(numout,*)
832         WRITE(numout,*)'obs_wri_vel :'
833         WRITE(numout,*)'~~~~~~~~~~~~~'
834         WRITE(numout,*)'Writing velocuty feedback file : ',TRIM(cfname)
835      ENDIF
836
837      ALLOCATE( &
838         & zu(profdata%nvprot(1)), &
839         & zv(profdata%nvprot(2))  &
840         & )
841      CALL obs_rotvel( profdata, k2dint, zu, zv )
842
843      ! Transform obs_prof data structure into obfbdata structure
844      fbdata%cdjuldref = '19500101000000'
845      DO jo = 1, profdata%nprof
846         fbdata%plam(jo)      = profdata%rlam(jo)
847         fbdata%pphi(jo)      = profdata%rphi(jo)
848         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
849         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
850         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
851         IF ( profdata%nqc(jo) > 10 ) THEN
852            fbdata%ioqc(jo)    = 4
853            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
854            fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10
855         ELSE
856            fbdata%ioqc(jo)    = profdata%nqc(jo)
857            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
858         ENDIF
859         fbdata%ipqc(jo)      = profdata%ipqc(jo)
860         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
861         fbdata%itqc(jo)      = profdata%itqc(jo)
862         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
863         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
864         fbdata%kindex(jo)    = profdata%npfil(jo)
865         DO jvar = 1, profdata%nvar
866            IF (ln_grid_global) THEN
867               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
868               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
869            ELSE
870               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
871               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
872            ENDIF
873         END DO
874         CALL greg2jul( 0, &
875            &           profdata%nmin(jo), &
876            &           profdata%nhou(jo), &
877            &           profdata%nday(jo), &
878            &           profdata%nmon(jo), &
879            &           profdata%nyea(jo), &
880            &           fbdata%ptim(jo),   &
881            &           krefdate = 19500101 )
882         ! Reform the profiles arrays for output
883         DO jvar = 1, 2
884            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
885               ik = profdata%var(jvar)%nvlidx(jk)
886               IF ( jvar == 1 ) THEN
887                  fbdata%padd(ik,jo,1,jvar) = zu(jk)
888               ELSE
889                  fbdata%padd(ik,jo,1,jvar) = zv(jk)
890               ENDIF
891               fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk)
892               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
893               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
894               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
895               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
896               IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN
897                  fbdata%ivlqc(ik,jo,jvar) = 4
898                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
899                  fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10
900               ELSE
901                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
902                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
903               ENDIF
904               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
905               DO ja = 1, nadd
906                  fbdata%padd(ik,jo,2+ja,jvar) = &
907                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
908               END DO
909               DO je = 1, next
910                  fbdata%pext(ik,jo,je) = &
911                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
912               END DO
913            END DO
914         END DO
915      END DO
916
917      ! Write the obfbdata structure
918      CALL write_obfbdata( cfname, fbdata )
919     
920      ! Output some basic statistics
921      CALL obs_wri_stats( fbdata )
922
923      CALL dealloc_obfbdata( fbdata )
924     
925      DEALLOCATE( &
926         & zu, &
927         & zv  &
928         & )
929
930   END SUBROUTINE obs_wri_vel
931
932   SUBROUTINE obs_wri_stats( fbdata )
933      !!-----------------------------------------------------------------------
934      !!
935      !!                     *** ROUTINE obs_wri_stats  ***
936      !!
937      !! ** Purpose : Output some basic statistics of the data being written out
938      !!
939      !! ** Method  :
940      !!
941      !! ** Action  :
942      !!
943      !!      ! 2014-08  (D. Lea) Initial version
944      !!-----------------------------------------------------------------------
945
946      !! * Arguments
947      TYPE(obfbdata) :: fbdata
948
949      !! * Local declarations
950      INTEGER :: jvar
951      INTEGER :: jo
952      INTEGER :: jk
953
954!      INTEGER :: nlev
955!      INTEGER :: nlevmpp
956!      INTEGER :: nobsmpp
957      INTEGER :: numgoodobs
958      INTEGER :: numgoodobsmpp
959      REAL(wp) :: zsumx
960      REAL(wp) :: zsumx2
961      REAL(wp) :: zomb
962
963      IF (lwp) THEN
964         WRITE(numout,*) ''
965         WRITE(numout,*) 'obs_wri_stats :'
966         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
967      ENDIF
968
969      DO jvar = 1, fbdata%nvar
970         zsumx=0.0_wp
971         zsumx2=0.0_wp
972         numgoodobs=0
973         DO jo = 1, fbdata%nobs
974            DO jk = 1, fbdata%nlev
975               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
976                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
977                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
978       
979             zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
980                  zsumx=zsumx+zomb
981                  zsumx2=zsumx2+zomb**2
982                  numgoodobs=numgoodobs+1
983          ENDIF
984            ENDDO
985         ENDDO
986
987         CALL obs_mpp_sum_integer( numgoodobs, numgoodobsmpp )
988         CALL mpp_sum(zsumx)
989         CALL mpp_sum(zsumx2)
990
991         IF (lwp) THEN
992       WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',numgoodobsmpp 
993       WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/numgoodobsmpp
994            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/numgoodobsmpp )
995       WRITE(numout,*) ''
996         ENDIF
997 
998      ENDDO
999
1000   END SUBROUTINE obs_wri_stats
1001
1002END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.