source: branches/UKMO/dev_r5518_seaiceobs/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 6283

Last change on this file since 6283 was 6283, checked in by drew, 6 years ago

Commit changes Jennie made at VN3.5.

File size: 35.3 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      !!      ! 2011-07 (D. Lea) Change SEAICE to ICECONC
601      !!-----------------------------------------------------------------------
602
603      !! * Modules used
604      IMPLICIT NONE
605
606      !! * Arguments
607      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files
608      TYPE(obs_surf), INTENT(INOUT) :: seaicedata   ! Full set of sea ice
609      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable
610      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info
611
612      !! * Local declarations
613      TYPE(obfbdata) :: fbdata
614      CHARACTER(LEN=40) :: cfname             ! netCDF filename
615      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_seaice'
616      INTEGER :: jo
617      INTEGER :: ja
618      INTEGER :: je
619      INTEGER :: nadd
620      INTEGER :: next
621
622      IF ( PRESENT( padd ) ) THEN
623         nadd = padd%inum
624      ELSE
625         nadd = 0
626      ENDIF
627
628      IF ( PRESENT( pext ) ) THEN
629         next = pext%inum
630      ELSE
631         next = 0
632      ENDIF
633
634      CALL init_obfbdata( fbdata )
635
636      CALL alloc_obfbdata( fbdata, 1, seaicedata%nsurf, 1, &
637         &                 1 + nadd, next, .TRUE. )
638
639      fbdata%cname(1)      = 'ICECONC'
640      fbdata%coblong(1)    = 'Sea ice concentration'
641      fbdata%cobunit(1)    = 'Fraction'
642      DO je = 1, next
643         fbdata%cextname(je) = pext%cdname(je)
644         fbdata%cextlong(je) = pext%cdlong(je,1)
645         fbdata%cextunit(je) = pext%cdunit(je,1)
646      END DO
647      fbdata%caddname(1)   = 'Hx'
648      fbdata%caddlong(1,1) = 'Model interpolated ICE'
649      fbdata%caddunit(1,1) = 'Fraction'
650      fbdata%cgrid(1)      = 'T'
651      DO ja = 1, nadd
652         fbdata%caddname(1+ja) = padd%cdname(ja)
653         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
654         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
655      END DO
656
657      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
658
659      IF(lwp) THEN
660         WRITE(numout,*)
661         WRITE(numout,*)'obs_wri_seaice :'
662         WRITE(numout,*)'~~~~~~~~~~~~~~~~'
663         WRITE(numout,*)'Writing SEAICE feedback file : ',TRIM(cfname)
664      ENDIF
665
666      ! Transform obs_prof data structure into obfbdata structure
667      fbdata%cdjuldref = '19500101000000'
668      DO jo = 1, seaicedata%nsurf
669         fbdata%plam(jo)      = seaicedata%rlam(jo)
670         fbdata%pphi(jo)      = seaicedata%rphi(jo)
671         WRITE(fbdata%cdtyp(jo),'(I4)') seaicedata%ntyp(jo)
672         fbdata%ivqc(jo,:)    = 0
673         fbdata%ivqcf(:,jo,:) = 0
674         IF ( seaicedata%nqc(jo) > 10 ) THEN
675            fbdata%ioqc(jo)    = 4
676            fbdata%ioqcf(1,jo) = 0
677            fbdata%ioqcf(2,jo) = seaicedata%nqc(jo) - 10
678         ELSE
679            fbdata%ioqc(jo)    = MAX(seaicedata%nqc(jo),1)
680            fbdata%ioqcf(:,jo) = 0
681         ENDIF
682         fbdata%ipqc(jo)      = 0
683         fbdata%ipqcf(:,jo)   = 0
684         fbdata%itqc(jo)      = 0
685         fbdata%itqcf(:,jo)   = 0
686         fbdata%cdwmo(jo)     = ''
687         fbdata%kindex(jo)    = seaicedata%nsfil(jo)
688         IF (ln_grid_global) THEN
689            fbdata%iobsi(jo,1) = seaicedata%mi(jo)
690            fbdata%iobsj(jo,1) = seaicedata%mj(jo)
691         ELSE
692            fbdata%iobsi(jo,1) = mig(seaicedata%mi(jo))
693            fbdata%iobsj(jo,1) = mjg(seaicedata%mj(jo))
694         ENDIF
695         CALL greg2jul( 0, &
696            &           seaicedata%nmin(jo), &
697            &           seaicedata%nhou(jo), &
698            &           seaicedata%nday(jo), &
699            &           seaicedata%nmon(jo), &
700            &           seaicedata%nyea(jo), &
701            &           fbdata%ptim(jo),   &
702            &           krefdate = 19500101 )
703         fbdata%padd(1,jo,1,1) = seaicedata%rmod(jo,1)
704         fbdata%pob(1,jo,1)    = seaicedata%robs(jo,1)
705         fbdata%pdep(1,jo)     = 0.0
706         fbdata%idqc(1,jo)     = 0
707         fbdata%idqcf(:,1,jo)  = 0
708         IF ( seaicedata%nqc(jo) > 10 ) THEN
709            fbdata%ivlqc(1,jo,1) = 4
710            fbdata%ivlqcf(1,1,jo,1) = 0
711            fbdata%ivlqcf(2,1,jo,1) = seaicedata%nqc(jo) - 10
712         ELSE
713            fbdata%ivlqc(1,jo,1) = MAX(seaicedata%nqc(jo),1)
714            fbdata%ivlqcf(:,1,jo,1) = 0
715         ENDIF
716         fbdata%iobsk(1,jo,1)  = 0
717         DO ja = 1, nadd
718            fbdata%padd(1,jo,1+ja,1) = &
719               & seaicedata%rext(jo,padd%ipoint(ja))
720         END DO
721         DO je = 1, next
722            fbdata%pext(1,jo,je) = &
723               & seaicedata%rext(jo,pext%ipoint(je))
724         END DO
725
726      END DO
727
728      ! Write the obfbdata structure
729      CALL write_obfbdata( cfname, fbdata )
730
731      ! Output some basic statistics
732      CALL obs_wri_stats( fbdata )
733
734      CALL dealloc_obfbdata( fbdata )
735
736   END SUBROUTINE obs_wri_seaice
737
738   SUBROUTINE obs_wri_vel( cprefix, profdata, k2dint, padd, pext )
739      !!-----------------------------------------------------------------------
740      !!
741      !!                     *** ROUTINE obs_wri_vel  ***
742      !!
743      !! ** Purpose : Write current (profile) observation
744      !!              related diagnostics
745      !!
746      !! ** Method  : NetCDF
747      !!
748      !! ** Action  :
749      !!
750      !! History :
751      !!      ! 09-01  (K. Mogensen) New feedback format routine
752      !!-----------------------------------------------------------------------
753
754      !! * Modules used
755
756      !! * Arguments
757      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files
758      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Full set of profile data
759      INTEGER, INTENT(IN) :: k2dint                 ! Horizontal interpolation method
760      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable
761      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info
762
763      !! * Local declarations
764      TYPE(obfbdata) :: fbdata
765      CHARACTER(LEN=40) :: cfname
766      INTEGER :: ilevel
767      INTEGER :: jvar
768      INTEGER :: jk
769      INTEGER :: ik
770      INTEGER :: jo
771      INTEGER :: ja
772      INTEGER :: je
773      INTEGER :: nadd
774      INTEGER :: next
775      REAL(wp) :: zpres
776      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
777         & zu, &
778         & zv
779
780      IF ( PRESENT( padd ) ) THEN
781         nadd = padd%inum
782      ELSE
783         nadd = 0
784      ENDIF
785
786      IF ( PRESENT( pext ) ) THEN
787         next = pext%inum
788      ELSE
789         next = 0
790      ENDIF
791
792      CALL init_obfbdata( fbdata )
793
794      ! Find maximum level
795      ilevel = 0
796      DO jvar = 1, 2
797         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
798      END DO
799      CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 2, 0, .TRUE. )
800
801      fbdata%cname(1)      = 'UVEL'
802      fbdata%cname(2)      = 'VVEL'
803      fbdata%coblong(1)    = 'Zonal velocity'
804      fbdata%coblong(2)    = 'Meridional velocity'
805      fbdata%cobunit(1)    = 'm/s'
806      fbdata%cobunit(2)    = 'm/s'
807      DO je = 1, next
808         fbdata%cextname(je) = pext%cdname(je)
809         fbdata%cextlong(je) = pext%cdlong(je,1)
810         fbdata%cextunit(je) = pext%cdunit(je,1)
811      END DO
812      fbdata%caddname(1)   = 'Hx'
813      fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
814      fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
815      fbdata%caddunit(1,1) = 'm/s'
816      fbdata%caddunit(1,2) = 'm/s'
817      fbdata%caddname(2)   = 'HxG'
818      fbdata%caddlong(2,1) = 'Model interpolated zonal velocity (model grid)'
819      fbdata%caddlong(2,2) = 'Model interpolated meridional velocity (model grid)'
820      fbdata%caddunit(2,1) = 'm/s'
821      fbdata%caddunit(2,2) = 'm/s' 
822      fbdata%cgrid(1)      = 'U' 
823      fbdata%cgrid(2)      = 'V'
824      DO ja = 1, nadd
825         fbdata%caddname(2+ja) = padd%cdname(ja)
826         fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
827         fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
828      END DO
829
830      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc
831
832      IF(lwp) THEN
833         WRITE(numout,*)
834         WRITE(numout,*)'obs_wri_vel :'
835         WRITE(numout,*)'~~~~~~~~~~~~~'
836         WRITE(numout,*)'Writing velocuty feedback file : ',TRIM(cfname)
837      ENDIF
838
839      ALLOCATE( &
840         & zu(profdata%nvprot(1)), &
841         & zv(profdata%nvprot(2))  &
842         & )
843      CALL obs_rotvel( profdata, k2dint, zu, zv )
844
845      ! Transform obs_prof data structure into obfbdata structure
846      fbdata%cdjuldref = '19500101000000'
847      DO jo = 1, profdata%nprof
848         fbdata%plam(jo)      = profdata%rlam(jo)
849         fbdata%pphi(jo)      = profdata%rphi(jo)
850         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
851         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
852         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
853         IF ( profdata%nqc(jo) > 10 ) THEN
854            fbdata%ioqc(jo)    = 4
855            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
856            fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10
857         ELSE
858            fbdata%ioqc(jo)    = profdata%nqc(jo)
859            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
860         ENDIF
861         fbdata%ipqc(jo)      = profdata%ipqc(jo)
862         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
863         fbdata%itqc(jo)      = profdata%itqc(jo)
864         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
865         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
866         fbdata%kindex(jo)    = profdata%npfil(jo)
867         DO jvar = 1, profdata%nvar
868            IF (ln_grid_global) THEN
869               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
870               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
871            ELSE
872               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
873               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
874            ENDIF
875         END DO
876         CALL greg2jul( 0, &
877            &           profdata%nmin(jo), &
878            &           profdata%nhou(jo), &
879            &           profdata%nday(jo), &
880            &           profdata%nmon(jo), &
881            &           profdata%nyea(jo), &
882            &           fbdata%ptim(jo),   &
883            &           krefdate = 19500101 )
884         ! Reform the profiles arrays for output
885         DO jvar = 1, 2
886            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
887               ik = profdata%var(jvar)%nvlidx(jk)
888               IF ( jvar == 1 ) THEN
889                  fbdata%padd(ik,jo,1,jvar) = zu(jk)
890               ELSE
891                  fbdata%padd(ik,jo,1,jvar) = zv(jk)
892               ENDIF
893               fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vmod(jk)
894               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
895               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
896               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
897               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
898               IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN
899                  fbdata%ivlqc(ik,jo,jvar) = 4
900                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
901                  fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10
902               ELSE
903                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
904                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
905               ENDIF
906               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
907               DO ja = 1, nadd
908                  fbdata%padd(ik,jo,2+ja,jvar) = &
909                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
910               END DO
911               DO je = 1, next
912                  fbdata%pext(ik,jo,je) = &
913                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
914               END DO
915            END DO
916         END DO
917      END DO
918
919      ! Write the obfbdata structure
920      CALL write_obfbdata( cfname, fbdata )
921     
922      ! Output some basic statistics
923      CALL obs_wri_stats( fbdata )
924
925      CALL dealloc_obfbdata( fbdata )
926     
927      DEALLOCATE( &
928         & zu, &
929         & zv  &
930         & )
931
932   END SUBROUTINE obs_wri_vel
933
934   SUBROUTINE obs_wri_stats( fbdata )
935      !!-----------------------------------------------------------------------
936      !!
937      !!                     *** ROUTINE obs_wri_stats  ***
938      !!
939      !! ** Purpose : Output some basic statistics of the data being written out
940      !!
941      !! ** Method  :
942      !!
943      !! ** Action  :
944      !!
945      !!      ! 2014-08  (D. Lea) Initial version
946      !!-----------------------------------------------------------------------
947
948      !! * Arguments
949      TYPE(obfbdata) :: fbdata
950
951      !! * Local declarations
952      INTEGER :: jvar
953      INTEGER :: jo
954      INTEGER :: jk
955
956!      INTEGER :: nlev
957!      INTEGER :: nlevmpp
958!      INTEGER :: nobsmpp
959      INTEGER :: numgoodobs
960      INTEGER :: numgoodobsmpp
961      REAL(wp) :: zsumx
962      REAL(wp) :: zsumx2
963      REAL(wp) :: zomb
964
965      IF (lwp) THEN
966         WRITE(numout,*) ''
967         WRITE(numout,*) 'obs_wri_stats :'
968         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
969      ENDIF
970
971      DO jvar = 1, fbdata%nvar
972         zsumx=0.0_wp
973         zsumx2=0.0_wp
974         numgoodobs=0
975         DO jo = 1, fbdata%nobs
976            DO jk = 1, fbdata%nlev
977               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
978                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
979                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
980       
981             zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
982                  zsumx=zsumx+zomb
983                  zsumx2=zsumx2+zomb**2
984                  numgoodobs=numgoodobs+1
985          ENDIF
986            ENDDO
987         ENDDO
988
989         CALL obs_mpp_sum_integer( numgoodobs, numgoodobsmpp )
990         CALL mpp_sum(zsumx)
991         CALL mpp_sum(zsumx2)
992
993         IF (lwp) THEN
994       WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',numgoodobsmpp 
995       WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/numgoodobsmpp
996            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/numgoodobsmpp )
997       WRITE(numout,*) ''
998         ENDIF
999 
1000      ENDDO
1001
1002   END SUBROUTINE obs_wri_stats
1003
1004END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.