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_r5518_obs_oper_update_obserr/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_obserr/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 8856

Last change on this file since 8856 was 8856, checked in by mattmartin, 7 years ago

Add first version of code to include obs error standard deviations in output files and read in from input files.

File size: 26.3 KB
Line 
1MODULE obs_write
2   !!======================================================================
3   !!                       ***  MODULE obs_write   ***
4   !! Observation diagnosticss: Write observation related diagnostics
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_wri_prof   : Write profile observations in feedback format
9   !!   obs_wri_surf   : Write surface observations in feedback format
10   !!   obs_wri_stats  : Print basic statistics on the data being written out
11   !!----------------------------------------------------------------------
12
13   !! * Modules used
14   USE par_kind, ONLY : &   ! Precision variables
15      & wp
16   USE in_out_manager       ! I/O manager
17   USE dom_oce              ! Ocean space and time domain variables
18   USE obs_types            ! Observation type integer to character translation
19   USE julian, ONLY : &         ! Julian date routines
20      & greg2jul
21   USE obs_utils, ONLY : &  ! Observation operator utility functions
22      & chkerr
23   USE obs_profiles_def     ! Type definitions for profiles
24   USE obs_surf_def         ! Type defintions for surface observations
25   USE obs_fbm              ! Observation feedback I/O
26   USE obs_grid             ! Grid tools
27   USE obs_conv             ! Conversion between units
28   USE obs_const
29   USE obs_mpp              ! MPP support routines for observation diagnostics
30   USE lib_mpp        ! MPP routines
31
32   IMPLICIT NONE
33
34   !! * Routine accessibility
35   PRIVATE
36   PUBLIC obs_wri_prof, &    ! Write profile observation files
37      &   obs_wri_surf, &    ! Write surface observation files
38      &   obswriinfo
39   
40   TYPE obswriinfo
41      INTEGER :: inum
42      INTEGER, POINTER, DIMENSION(:) :: ipoint
43      CHARACTER(len=ilenname), POINTER, DIMENSION(:) :: cdname
44      CHARACTER(len=ilenlong), POINTER, DIMENSION(:,:) :: cdlong
45      CHARACTER(len=ilenunit), POINTER, DIMENSION(:,:) :: cdunit
46   END TYPE obswriinfo
47
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE obs_wri_prof( profdata, padd, pext )
57      !!-----------------------------------------------------------------------
58      !!
59      !!                     *** ROUTINE obs_wri_prof  ***
60      !!
61      !! ** Purpose : Write profile feedback files
62      !!
63      !! ** Method  : NetCDF
64      !!
65      !! ** Action  :
66      !!
67      !! History :
68      !!      ! 06-04  (A. Vidard) Original
69      !!      ! 06-04  (A. Vidard) Reformatted
70      !!      ! 06-10  (A. Weaver) Cleanup
71      !!      ! 07-01  (K. Mogensen) Use profile data types
72      !!      ! 07-03  (K. Mogensen) General handling of profiles
73      !!      ! 09-01  (K. Mogensen) New feedback format
74      !!      ! 15-02  (M. Martin) Combined routine for writing profiles
75      !!-----------------------------------------------------------------------
76
77      !! * Arguments
78      TYPE(obs_prof), INTENT(INOUT) :: profdata      ! Full set of profile data
79      TYPE(obswriinfo), OPTIONAL :: padd             ! Additional info for each variable
80      TYPE(obswriinfo), OPTIONAL :: pext             ! Extra info
81
82      !! * Local declarations
83      TYPE(obfbdata) :: fbdata
84      CHARACTER(LEN=40) :: clfname
85      CHARACTER(LEN=10) :: clfiletype
86      INTEGER :: ilevel
87      INTEGER :: jvar
88      INTEGER :: jo
89      INTEGER :: jk
90      INTEGER :: ik
91      INTEGER :: ja
92      INTEGER :: je
93      INTEGER :: iadd
94      INTEGER :: iext
95      REAL(wp) :: zpres
96
97      IF ( PRESENT( padd ) ) THEN
98         iadd = padd%inum
99      ELSE
100         iadd = 0
101      ENDIF
102
103      IF ( PRESENT( pext ) ) THEN
104         iext = pext%inum
105      ELSE
106         iext = 0
107      ENDIF
108
109      CALL init_obfbdata( fbdata )
110
111      ! Find maximum level
112      ilevel = 0
113      DO jvar = 1, 2
114         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
115      END DO
116
117      SELECT CASE ( TRIM(profdata%cvars(1)) )
118      CASE('POTM')
119
120         clfiletype='profb'
121         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
122            &                 1 + iadd, 1 + iext, .TRUE. )
123         fbdata%cname(1)      = profdata%cvars(1)
124         fbdata%cname(2)      = profdata%cvars(2)
125         fbdata%coblong(1)    = 'Potential temperature'
126         fbdata%coblong(2)    = 'Practical salinity'
127         fbdata%cobunit(1)    = 'Degrees centigrade'
128         fbdata%cobunit(2)    = 'PSU'
129         fbdata%cextname(1)   = 'TEMP'
130         fbdata%cextlong(1)   = 'Insitu temperature'
131         fbdata%cextunit(1)   = 'Degrees centigrade'
132         fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
133         fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
134         fbdata%caddunit(1,1) = 'Degrees centigrade'
135         fbdata%caddunit(1,2) = 'PSU'
136         fbdata%cgrid(:)      = 'T'
137         DO je = 1, iext
138            fbdata%cextname(1+je) = pext%cdname(je)
139            fbdata%cextlong(1+je) = pext%cdlong(je,1)
140            fbdata%cextunit(1+je) = pext%cdunit(je,1)
141         END DO
142         DO ja = 1, iadd
143            fbdata%caddname(1+ja) = padd%cdname(ja)
144            DO jvar = 1, 2
145               fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
146               fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)
147            END DO
148         END DO
149
150      CASE('UVEL')
151
152         clfiletype='velfb'
153         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 1, 0, .TRUE. )
154         fbdata%cname(1)      = profdata%cvars(1)
155         fbdata%cname(2)      = profdata%cvars(2)
156         fbdata%coblong(1)    = 'Zonal velocity'
157         fbdata%coblong(2)    = 'Meridional velocity'
158         fbdata%cobunit(1)    = 'm/s'
159         fbdata%cobunit(2)    = 'm/s'
160         DO je = 1, iext
161            fbdata%cextname(je) = pext%cdname(je)
162            fbdata%cextlong(je) = pext%cdlong(je,1)
163            fbdata%cextunit(je) = pext%cdunit(je,1)
164         END DO
165         fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
166         fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
167         fbdata%caddunit(1,1) = 'm/s'
168         fbdata%caddunit(1,2) = 'm/s'
169         fbdata%cgrid(1)      = 'U' 
170         fbdata%cgrid(2)      = 'V'
171         DO ja = 1, iadd
172            fbdata%caddname(1+ja) = padd%cdname(ja)
173            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
174            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
175         END DO
176
177      END SELECT
178
179      fbdata%caddname(1)   = 'Hx'
180
181      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
182
183      IF(lwp) THEN
184         WRITE(numout,*)
185         WRITE(numout,*)'obs_wri_prof :'
186         WRITE(numout,*)'~~~~~~~~~~~~~'
187         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
188      ENDIF
189
190      ! Transform obs_prof data structure into obfb data structure
191      fbdata%cdjuldref = '19500101000000'
192      DO jo = 1, profdata%nprof
193         fbdata%plam(jo)      = profdata%rlam(jo)
194         fbdata%pphi(jo)      = profdata%rphi(jo)
195         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
196         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
197         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
198         IF ( profdata%nqc(jo) > 255 ) THEN
199            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2)
200            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
201            fbdata%ioqcf(2,jo) = profdata%nqc(jo)
202         ELSE
203            fbdata%ioqc(jo)    = profdata%nqc(jo)
204            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
205         ENDIF
206         fbdata%ipqc(jo)      = profdata%ipqc(jo)
207         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
208         fbdata%itqc(jo)      = profdata%itqc(jo)
209         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
210         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
211         fbdata%kindex(jo)    = profdata%npfil(jo)
212         DO jvar = 1, profdata%nvar
213            IF (ln_grid_global) THEN
214               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
215               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
216            ELSE
217               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
218               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
219            ENDIF
220         END DO
221         CALL greg2jul( 0, &
222            &           profdata%nmin(jo), &
223            &           profdata%nhou(jo), &
224            &           profdata%nday(jo), &
225            &           profdata%nmon(jo), &
226            &           profdata%nyea(jo), &
227            &           fbdata%ptim(jo),   &
228            &           krefdate = 19500101 )
229         ! Reform the profiles arrays for output
230         DO jvar = 1, 2
231            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
232               ik = profdata%var(jvar)%nvlidx(jk)
233               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
234               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
235               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
236               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
237               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
238               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
239                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
240                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
241                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111')
242               ELSE
243                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
244                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
245               ENDIF
246               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
247               DO ja = 1, iadd
248                  fbdata%padd(ik,jo,1+ja,jvar) = &
249                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
250               END DO
251               DO je = 1, iext
252                  fbdata%pext(ik,jo,1+je) = &
253                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
254               END DO
255               IF ( ( jvar == 1 ) .AND. &
256                  & ( TRIM(profdata%cvars(1)) == 'POTM' ) ) THEN
257                  fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
258               ENDIF
259            END DO
260         END DO
261      END DO
262
263      IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
264         ! Convert insitu temperature to potential temperature using the model
265         ! salinity if no potential temperature
266         DO jo = 1, fbdata%nobs
267            IF ( fbdata%pphi(jo) < 9999.0 ) THEN
268               DO jk = 1, fbdata%nlev
269                  IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
270                     & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
271                     & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
272                     & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
273                     zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
274                        &              REAL(fbdata%pphi(jo),wp) )
275                     fbdata%pob(jk,jo,1) = potemp( &
276                        &                     REAL(fbdata%padd(jk,jo,1,2), wp),  &
277                        &                     REAL(fbdata%pext(jk,jo,1), wp), &
278                        &                     zpres, 0.0_wp )
279                  ENDIF
280               END DO
281            ENDIF
282         END DO
283      ENDIF
284
285      ! Write the obfbdata structure
286      CALL write_obfbdata( clfname, fbdata )
287
288      ! Output some basic statistics
289      CALL obs_wri_stats( fbdata )
290
291      CALL dealloc_obfbdata( fbdata )
292
293   END SUBROUTINE obs_wri_prof
294
295   SUBROUTINE obs_wri_surf( surfdata, padd, pext )
296      !!-----------------------------------------------------------------------
297      !!
298      !!                     *** ROUTINE obs_wri_surf  ***
299      !!
300      !! ** Purpose : Write surface observation files
301      !!
302      !! ** Method  : NetCDF
303      !!
304      !! ** Action  :
305      !!
306      !!      ! 07-03  (K. Mogensen) Original
307      !!      ! 09-01  (K. Mogensen) New feedback format.
308      !!      ! 15-02  (M. Martin) Combined surface writing routine.
309      !!-----------------------------------------------------------------------
310
311      !! * Modules used
312      IMPLICIT NONE
313
314      !! * Arguments
315      TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data
316      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
317      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info
318
319      !! * Local declarations
320      TYPE(obfbdata) :: fbdata
321      CHARACTER(LEN=40) :: clfname         ! netCDF filename
322      CHARACTER(LEN=10) :: clfiletype
323      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
324      INTEGER :: jo
325      INTEGER :: ja
326      INTEGER :: je
327      INTEGER :: iadd
328      INTEGER :: iext
329      INTEGER :: indx_std
330      INTEGER :: iadd_std
331
332      IF ( PRESENT( padd ) ) THEN
333         iadd = padd%inum
334      ELSE
335         iadd = 0
336      ENDIF
337
338      IF ( PRESENT( pext ) ) THEN
339         iext = pext%inum
340      ELSE
341         iext = 0
342      ENDIF
343
344      IF ( surfdata%nextra > 0 ) THEN
345
346         indx_std = -1
347         iadd_std = 0
348         DO je = 1, surfdata%nextra
349           IF ( TRIM( surfdata%cext(je) ) == 'STD' ) THEN
350             iadd_std = 1
351             indx_std = je
352           ENDIF
353         END DO
354
355      CALL init_obfbdata( fbdata )
356
357      SELECT CASE ( TRIM(surfdata%cvars(1)) )
358      CASE('SLA')
359
360         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
361            &                 2 + iadd + iadd_std, 1 + iext, .TRUE. )
362
363         clfiletype = 'slafb'
364         fbdata%cname(1)      = surfdata%cvars(1)
365         fbdata%coblong(1)    = 'Sea level anomaly'
366         fbdata%cobunit(1)    = 'Metres'
367         fbdata%cextname(1)   = 'MDT'
368         fbdata%cextlong(1)   = 'Mean dynamic topography'
369         fbdata%cextunit(1)   = 'Metres'
370         DO je = 1, iext
371            fbdata%cextname(je) = pext%cdname(je)
372            fbdata%cextlong(je) = pext%cdlong(je,1)
373            fbdata%cextunit(je) = pext%cdunit(je,1)
374         END DO
375         fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
376         fbdata%caddunit(1,1) = 'Metres' 
377         fbdata%caddname(2)   = 'SSH'
378         fbdata%caddlong(2,1) = 'Model Sea surface height'
379         fbdata%caddunit(2,1) = 'Metres'
380         fbdata%cgrid(1)      = 'T'
381         DO ja = 1, iadd
382            fbdata%caddname(2+iadd_std+ja) = padd%cdname(ja)
383            fbdata%caddlong(2+iadd_std+ja,1) = padd%cdlong(ja,1)
384            fbdata%caddunit(2+iadd_std+ja,1) = padd%cdunit(ja,1)
385         END DO
386
387      CASE('SST')
388
389         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
390            &                 1 + iadd + iadd_std, iext, .TRUE. )
391
392         clfiletype = 'sstfb'
393         fbdata%cname(1)      = surfdata%cvars(1)
394         fbdata%coblong(1)    = 'Sea surface temperature'
395         fbdata%cobunit(1)    = 'Degree centigrade'
396         DO je = 1, iext
397            fbdata%cextname(je) = pext%cdname(je)
398            fbdata%cextlong(je) = pext%cdlong(je,1)
399            fbdata%cextunit(je) = pext%cdunit(je,1)
400         END DO
401         fbdata%caddlong(1,1) = 'Model interpolated SST'
402         fbdata%caddunit(1,1) = 'Degree centigrade'
403         fbdata%cgrid(1)      = 'T'
404         DO ja = 1, iadd
405            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
406            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
407            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
408         END DO
409
410      CASE('ICECONC')
411
412         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
413            &                 1 + iadd + iadd_std, iext, .TRUE. )
414
415         clfiletype = 'sicfb'
416         fbdata%cname(1)      = surfdata%cvars(1)
417         fbdata%coblong(1)    = 'Sea ice'
418         fbdata%cobunit(1)    = 'Fraction'
419         DO je = 1, iext
420            fbdata%cextname(je) = pext%cdname(je)
421            fbdata%cextlong(je) = pext%cdlong(je,1)
422            fbdata%cextunit(je) = pext%cdunit(je,1)
423         END DO
424         fbdata%caddlong(1,1) = 'Model interpolated ICE'
425         fbdata%caddunit(1,1) = 'Fraction'
426         fbdata%cgrid(1)      = 'T'
427         DO ja = 1, iadd
428            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
429            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
430            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
431         END DO
432
433      CASE('SSS')
434
435         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
436            &                 1 + iadd + iadd_std, iext, .TRUE. )
437
438         clfiletype = 'sssfb'
439         fbdata%cname(1)      = surfdata%cvars(1)
440         fbdata%coblong(1)    = 'Sea surface salinity'
441         fbdata%cobunit(1)    = 'psu'
442         DO je = 1, iext
443            fbdata%cextname(je) = pext%cdname(je)
444            fbdata%cextlong(je) = pext%cdlong(je,1)
445            fbdata%cextunit(je) = pext%cdunit(je,1)
446         END DO
447         fbdata%caddlong(1,1) = 'Model interpolated SSS'
448         fbdata%caddunit(1,1) = 'psu'
449         fbdata%cgrid(1)      = 'T'
450         DO ja = 1, iadd
451            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
452            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
453            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
454         END DO
455
456      CASE('LOGCHL','LogChl','logchl')
457
458         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
459            &                 1 + iadd + iadd_std, iext, .TRUE. )
460
461         clfiletype = 'logchlfb'
462         fbdata%cname(1)      = surfdata%cvars(1)
463         fbdata%coblong(1)    = 'logchl concentration'
464         fbdata%cobunit(1)    = 'mg/m3'
465         DO je = 1, iext
466            fbdata%cextname(je) = pext%cdname(je)
467            fbdata%cextlong(je) = pext%cdlong(je,1)
468            fbdata%cextunit(je) = pext%cdunit(je,1)
469         END DO
470         fbdata%caddlong(1,1) = 'Model interpolated LOGCHL'
471         fbdata%caddunit(1,1) = 'mg/m3'
472         fbdata%cgrid(1)      = 'T'
473         DO ja = 1, iadd
474            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
475            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
476            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
477         END DO
478
479      CASE('SPM','Spm','spm')
480
481         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
482            &                 1 + iadd + iadd_std, iext, .TRUE. )
483
484         clfiletype = 'spmfb'
485         fbdata%cname(1)      = surfdata%cvars(1)
486         fbdata%coblong(1)    = 'spm'
487         fbdata%cobunit(1)    = 'g/m3'
488         DO je = 1, iext
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%caddlong(1,1) = 'Model interpolated spm'
494         fbdata%caddunit(1,1) = 'g/m3'
495         fbdata%cgrid(1)      = 'T'
496         DO ja = 1, iadd
497            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
498            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
499            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
500         END DO
501
502      CASE('FCO2','fCO2','fco2')
503
504         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
505            &                 1 + iadd + iadd_std, iext, .TRUE. )
506
507         clfiletype = 'fco2fb'
508         fbdata%cname(1)      = surfdata%cvars(1)
509         fbdata%coblong(1)    = 'fco2'
510         fbdata%cobunit(1)    = 'uatm'
511         DO je = 1, iext
512            fbdata%cextname(je) = pext%cdname(je)
513            fbdata%cextlong(je) = pext%cdlong(je,1)
514            fbdata%cextunit(je) = pext%cdunit(je,1)
515         END DO
516         fbdata%caddlong(1,1) = 'Model interpolated fco2'
517         fbdata%caddunit(1,1) = 'uatm'
518         fbdata%cgrid(1)      = 'T'
519         DO ja = 1, iadd
520            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
521            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
522            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
523         END DO
524
525      CASE('PCO2','pCO2','pco2')
526
527         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
528            &                 1 + iadd + iadd_std, iext, .TRUE. )
529
530         clfiletype = 'pco2fb'
531         fbdata%cname(1)      = surfdata%cvars(1)
532         fbdata%coblong(1)    = 'pco2'
533         fbdata%cobunit(1)    = 'uatm'
534         DO je = 1, iext
535            fbdata%cextname(je) = pext%cdname(je)
536            fbdata%cextlong(je) = pext%cdlong(je,1)
537            fbdata%cextunit(je) = pext%cdunit(je,1)
538         END DO
539         fbdata%caddlong(1,1) = 'Model interpolated pco2'
540         fbdata%caddunit(1,1) = 'uatm'
541         fbdata%cgrid(1)      = 'T'
542         DO ja = 1, iadd
543            fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja)
544            fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1)
545            fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1)
546         END DO
547
548      CASE DEFAULT
549
550         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
551
552      END SELECT
553
554      fbdata%caddname(1)   = 'Hx'
555      IF ( indx_std /= -1 ) THEN
556         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1
557         fbdata%caddname(1+iadd_std)   = surfdata%cext(indx_std)
558         fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation'
559         fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1)
560      ENDIF
561
562      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
563
564      IF(lwp) THEN
565         WRITE(numout,*)
566         WRITE(numout,*)'obs_wri_surf :'
567         WRITE(numout,*)'~~~~~~~~~~~~~'
568         WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
569      ENDIF
570
571      ! Transform surf data structure into obfbdata structure
572      fbdata%cdjuldref = '19500101000000'
573      DO jo = 1, surfdata%nsurf
574         fbdata%plam(jo)      = surfdata%rlam(jo)
575         fbdata%pphi(jo)      = surfdata%rphi(jo)
576         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
577         fbdata%ivqc(jo,:)    = 0
578         fbdata%ivqcf(:,jo,:) = 0
579         IF ( surfdata%nqc(jo) > 255 ) THEN
580            fbdata%ioqc(jo)    = 4
581            fbdata%ioqcf(1,jo) = 0
582            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
583         ELSE
584            fbdata%ioqc(jo)    = surfdata%nqc(jo)
585            fbdata%ioqcf(:,jo) = 0
586         ENDIF
587         fbdata%ipqc(jo)      = 0
588         fbdata%ipqcf(:,jo)   = 0
589         fbdata%itqc(jo)      = 0
590         fbdata%itqcf(:,jo)   = 0
591         fbdata%cdwmo(jo)     = surfdata%cwmo(jo)
592         fbdata%kindex(jo)    = surfdata%nsfil(jo)
593         IF (ln_grid_global) THEN
594            fbdata%iobsi(jo,1) = surfdata%mi(jo)
595            fbdata%iobsj(jo,1) = surfdata%mj(jo)
596         ELSE
597            fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
598            fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
599         ENDIF
600         CALL greg2jul( 0, &
601            &           surfdata%nmin(jo), &
602            &           surfdata%nhou(jo), &
603            &           surfdata%nday(jo), &
604            &           surfdata%nmon(jo), &
605            &           surfdata%nyea(jo), &
606            &           fbdata%ptim(jo),   &
607            &           krefdate = 19500101 )
608         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
609         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
610         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1) 
611         fbdata%pdep(1,jo)     = 0.0
612         fbdata%idqc(1,jo)     = 0
613         fbdata%idqcf(:,1,jo)  = 0
614         IF ( surfdata%nqc(jo) > 255 ) THEN
615            fbdata%ivqc(jo,1)       = 4
616            fbdata%ivlqc(1,jo,1)    = 4
617            fbdata%ivlqcf(1,1,jo,1) = 0
618            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
619         ELSE
620            fbdata%ivqc(jo,1)       = surfdata%nqc(jo)
621            fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo)
622            fbdata%ivlqcf(:,1,jo,1) = 0
623         ENDIF
624         fbdata%iobsk(1,jo,1)  = 0
625         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
626         IF ( indx_std /= -1 ) THEN
627            fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std)
628         ENDIF
629
630         DO ja = 1, iadd
631            fbdata%padd(1,jo,2+iadd_std+ja,1) = &
632               & surfdata%rext(jo,padd%ipoint(ja))
633         END DO
634         DO je = 1, iext
635            fbdata%pext(1,jo,1+je) = &
636               & surfdata%rext(jo,pext%ipoint(je))
637         END DO
638      END DO
639
640      ! Write the obfbdata structure
641      CALL write_obfbdata( clfname, fbdata )
642
643      ! Output some basic statistics
644      CALL obs_wri_stats( fbdata )
645
646      CALL dealloc_obfbdata( fbdata )
647
648   END SUBROUTINE obs_wri_surf
649
650   SUBROUTINE obs_wri_stats( fbdata )
651      !!-----------------------------------------------------------------------
652      !!
653      !!                     *** ROUTINE obs_wri_stats  ***
654      !!
655      !! ** Purpose : Output some basic statistics of the data being written out
656      !!
657      !! ** Method  :
658      !!
659      !! ** Action  :
660      !!
661      !!      ! 2014-08  (D. Lea) Initial version
662      !!-----------------------------------------------------------------------
663
664      !! * Arguments
665      TYPE(obfbdata) :: fbdata
666
667      !! * Local declarations
668      INTEGER :: jvar
669      INTEGER :: jo
670      INTEGER :: jk
671      INTEGER :: inumgoodobs
672      INTEGER :: inumgoodobsmpp
673      REAL(wp) :: zsumx
674      REAL(wp) :: zsumx2
675      REAL(wp) :: zomb
676     
677
678      IF (lwp) THEN
679         WRITE(numout,*) ''
680         WRITE(numout,*) 'obs_wri_stats :'
681         WRITE(numout,*) '~~~~~~~~~~~~~~~'
682      ENDIF
683
684      DO jvar = 1, fbdata%nvar
685         zsumx=0.0_wp
686         zsumx2=0.0_wp
687         inumgoodobs=0
688         DO jo = 1, fbdata%nobs
689            DO jk = 1, fbdata%nlev
690               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
691                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
692                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
693
694                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
695                  zsumx=zsumx+zomb
696                  zsumx2=zsumx2+zomb**2
697                  inumgoodobs=inumgoodobs+1
698               ENDIF
699            ENDDO
700         ENDDO
701
702         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
703         CALL mpp_sum(zsumx)
704         CALL mpp_sum(zsumx2)
705
706         IF (lwp) THEN
707            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
708            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
709            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
710            WRITE(numout,*) ''
711         ENDIF
712
713      ENDDO
714
715   END SUBROUTINE obs_wri_stats
716
717END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.