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 branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 6855

Last change on this file since 6855 was 6855, checked in by dford, 8 years ago

Initial implementation of observation operator for SPM.

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