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

source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

Last change on this file was 15764, checked in by jcastill, 2 years ago

Update with the latest changes in branches/UKMO/dev_r5518_obs_oper_update@15400
File size: 31.5 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      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
87      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
88      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
89      INTEGER :: ilevel
90      INTEGER :: jvar
91      INTEGER :: jo
92      INTEGER :: jk
93      INTEGER :: ik
94      INTEGER :: ja
95      INTEGER :: je
96      INTEGER :: iadd
97      INTEGER :: iadd_clm ! 1 if climatology present
98      INTEGER :: iext
99      REAL(wp) :: zpres
100
101      iadd_clm = 0 
102      IF ( profdata%lclim ) iadd_clm = 1
103
104      IF ( PRESENT( padd ) ) THEN
105         iadd = padd%inum
106      ELSE
107         iadd = 0
108      ENDIF
109
110      IF ( PRESENT( pext ) ) THEN
111         iext = pext%inum
112      ELSE
113         iext = 0
114      ENDIF
115
116      CALL init_obfbdata( fbdata )
117
118      ! Find maximum level
119      ilevel = 0
120      DO jvar = 1, profdata%nvar
121         ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
122      END DO
123
124      SELECT CASE ( TRIM(profdata%cvars(1)) )
125      CASE('POTM')
126
127         clfiletype='profb'
128         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
129            &                 1 + iadd, 1 + iext, .TRUE. )
130         fbdata%cname(1)      = profdata%cvars(1)
131         fbdata%cname(2)      = profdata%cvars(2)
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         fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
140         fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
141         fbdata%caddunit(1,1) = 'Degrees centigrade'
142         fbdata%caddunit(1,2) = 'PSU'
143         IF ( profdata%lclim ) THEN
144            fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature' 
145            fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity' 
146            fbdata%caddunit(2,1) = 'Degrees centigrade' 
147            fbdata%caddunit(2,2) = 'PSU' 
148         ENDIF
149         fbdata%cgrid(:)      = 'T'
150         DO je = 1, iext
151            fbdata%cextname(1+je) = pext%cdname(je)
152            fbdata%cextlong(1+je) = pext%cdlong(je,1)
153            fbdata%cextunit(1+je) = pext%cdunit(je,1)
154         END DO
155         DO ja = 1, iadd
156            fbdata%caddname(1+ja) = padd%cdname(ja)
157            DO jvar = 1, 2
158               fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
159               fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)
160            END DO
161         END DO
162
163      CASE('UVEL')
164
165         clfiletype='velfb'
166         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 1, 0, .TRUE. )
167         fbdata%cname(1)      = profdata%cvars(1)
168         fbdata%cname(2)      = profdata%cvars(2)
169         fbdata%coblong(1)    = 'Zonal velocity'
170         fbdata%coblong(2)    = 'Meridional velocity'
171         fbdata%cobunit(1)    = 'm/s'
172         fbdata%cobunit(2)    = 'm/s'
173         DO je = 1, iext
174            fbdata%cextname(je) = pext%cdname(je)
175            fbdata%cextlong(je) = pext%cdlong(je,1)
176            fbdata%cextunit(je) = pext%cdunit(je,1)
177         END DO
178         fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
179         fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
180         fbdata%caddunit(1,1) = 'm/s'
181         fbdata%caddunit(1,2) = 'm/s'
182         IF ( profdata%lclim ) THEN
183            fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity' 
184            fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity' 
185            fbdata%caddunit(2,1) = 'm/s' 
186            fbdata%caddunit(2,2) = 'm/s' 
187         ENDIF
188         fbdata%cgrid(1)      = 'U' 
189         fbdata%cgrid(2)      = 'V'
190         DO ja = 1, iadd
191            fbdata%caddname(1+ja) = padd%cdname(ja)
192            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
193            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
194         END DO
195
196      CASE('PLCHLTOT')
197
198         clfiletype = 'plchltotfb'
199         cllongname = 'log10(chlorophyll concentration)'
200         clunits    = 'log10(mg/m3)'
201         clgrid     = 'T'
202
203      CASE('PCHLTOT')
204
205         clfiletype = 'pchltotfb'
206         cllongname = 'chlorophyll concentration'
207         clunits    = 'mg/m3'
208         clgrid     = 'T'
209
210      CASE('PNO3')
211
212         clfiletype = 'pno3fb'
213         cllongname = 'nitrate'
214         clunits    = 'mmol/m3'
215         clgrid     = 'T'
216
217      CASE('PSI4')
218
219         clfiletype = 'psi4fb'
220         cllongname = 'silicate'
221         clunits    = 'mmol/m3'
222         clgrid     = 'T'
223
224      CASE('PPO4')
225
226         clfiletype = 'ppo4fb'
227         cllongname = 'phosphate'
228         clunits    = 'mmol/m3'
229         clgrid     = 'T'
230
231      CASE('PDIC')
232
233         clfiletype = 'pdicfb'
234         cllongname = 'dissolved inorganic carbon'
235         clunits    = 'mmol/m3'
236         clgrid     = 'T'
237
238      CASE('PALK')
239
240         clfiletype = 'palkfb'
241         cllongname = 'alkalinity'
242         clunits    = 'meq/m3'
243         clgrid     = 'T'
244
245      CASE('PPH')
246
247         clfiletype = 'pphfb'
248         cllongname = 'pH'
249         clunits    = '-'
250         clgrid     = 'T'
251
252      CASE('PO2')
253
254         clfiletype = 'po2fb'
255         cllongname = 'dissolved oxygen'
256         clunits    = 'mmol/m3'
257         clgrid     = 'T'
258
259      END SELECT
260     
261      IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. &
262         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN
263         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, &
264            &                 1 + iadd, iext, .TRUE. )
265         fbdata%cname(1)      = profdata%cvars(1)
266         fbdata%coblong(1)    = cllongname
267         fbdata%cobunit(1)    = clunits
268         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname)
269         fbdata%caddunit(1,1) = clunits
270         IF ( profdata%lclim ) THEN
271            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname) 
272            fbdata%caddunit(2,1) = clunits 
273         ENDIF
274         fbdata%cgrid(:)      = clgrid
275         DO je = 1, iext
276            fbdata%cextname(je) = pext%cdname(je)
277            fbdata%cextlong(je) = pext%cdlong(je,1)
278            fbdata%cextunit(je) = pext%cdunit(je,1)
279         END DO
280         DO ja = 1, iadd
281            fbdata%caddname(1+ja) = padd%cdname(ja)
282            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
283            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
284         END DO
285      ENDIF
286
287      fbdata%caddname(1)   = 'Hx'
288
289      IF ( profdata%lclim ) fbdata%caddname(1+iadd_clm)   = 'CLM'
290
291      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
292
293      IF(lwp) THEN
294         WRITE(numout,*)
295         WRITE(numout,*)'obs_wri_prof :'
296         WRITE(numout,*)'~~~~~~~~~~~~~'
297         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
298      ENDIF
299
300      ! Transform obs_prof data structure into obfb data structure
301      fbdata%cdjuldref = '19500101000000'
302      DO jo = 1, profdata%nprof
303         fbdata%plam(jo)      = profdata%rlam(jo)
304         fbdata%pphi(jo)      = profdata%rphi(jo)
305         WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
306         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:)
307         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
308         IF ( profdata%nqc(jo) > 255 ) THEN
309            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2)
310            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
311            fbdata%ioqcf(2,jo) = profdata%nqc(jo)
312         ELSE
313            fbdata%ioqc(jo)    = profdata%nqc(jo)
314            fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
315         ENDIF
316         fbdata%ipqc(jo)      = profdata%ipqc(jo)
317         fbdata%ipqcf(:,jo)   = profdata%ipqcf(:,jo)
318         fbdata%itqc(jo)      = profdata%itqc(jo)
319         fbdata%itqcf(:,jo)   = profdata%itqcf(:,jo)
320         fbdata%cdwmo(jo)     = profdata%cwmo(jo)
321         fbdata%kindex(jo)    = profdata%npfil(jo)
322         DO jvar = 1, profdata%nvar
323            IF (ln_grid_global) THEN
324               fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
325               fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
326            ELSE
327               fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
328               fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
329            ENDIF
330         END DO
331         CALL greg2jul( 0, &
332            &           profdata%nmin(jo), &
333            &           profdata%nhou(jo), &
334            &           profdata%nday(jo), &
335            &           profdata%nmon(jo), &
336            &           profdata%nyea(jo), &
337            &           fbdata%ptim(jo),   &
338            &           krefdate = 19500101 )
339         ! Reform the profiles arrays for output
340         DO jvar = 1, profdata%nvar
341            DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
342               ik = profdata%var(jvar)%nvlidx(jk)
343               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk)
344               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk)
345               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk)
346               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk)
347               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
348                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
349                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
350                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111')
351               ELSE
352                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
353                  fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
354               ENDIF
355               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk)
356               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 
357               IF ( profdata%lclim ) THEN           
358                  fbdata%padd(ik,jo,1+iadd_clm,jvar) = profdata%var(jvar)%vclm(jk)     
359               ENDIF               
360               DO ja = 1, iadd 
361                  fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = &
362                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja))
363               END DO
364               DO je = 1, iext
365                  fbdata%pext(ik,jo,1+je) = &
366                     & profdata%var(jvar)%vext(jk,pext%ipoint(je))
367               END DO
368               IF ( ( jvar == 1 ) .AND. &
369                  & ( TRIM(profdata%cvars(1)) == 'POTM' ) ) THEN
370                  fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
371               ENDIF
372            END DO
373         END DO
374      END DO
375
376      IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
377         ! Convert insitu temperature to potential temperature using the model
378         ! salinity if no potential temperature
379         DO jo = 1, fbdata%nobs
380            IF ( fbdata%pphi(jo) < 9999.0 ) THEN
381               DO jk = 1, fbdata%nlev
382                  IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
383                     & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
384                     & ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
385                     & ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
386                     zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
387                        &              REAL(fbdata%pphi(jo),wp) )
388                     fbdata%pob(jk,jo,1) = potemp( &
389                        &                     REAL(fbdata%padd(jk,jo,1,2), wp),  &
390                        &                     REAL(fbdata%pext(jk,jo,1), wp), &
391                        &                     zpres, 0.0_wp )
392                  ENDIF
393               END DO
394            ENDIF
395         END DO
396      ENDIF
397
398      ! Write the obfbdata structure
399      CALL write_obfbdata( clfname, fbdata )
400
401      ! Output some basic statistics
402      CALL obs_wri_stats( fbdata )
403
404      CALL dealloc_obfbdata( fbdata )
405
406   END SUBROUTINE obs_wri_prof
407
408   SUBROUTINE obs_wri_surf( surfdata, padd, pext )
409      !!-----------------------------------------------------------------------
410      !!
411      !!                     *** ROUTINE obs_wri_surf  ***
412      !!
413      !! ** Purpose : Write surface observation files
414      !!
415      !! ** Method  : NetCDF
416      !!
417      !! ** Action  :
418      !!
419      !!      ! 07-03  (K. Mogensen) Original
420      !!      ! 09-01  (K. Mogensen) New feedback format.
421      !!      ! 15-02  (M. Martin) Combined surface writing routine.
422      !!-----------------------------------------------------------------------
423
424      !! * Modules used
425      IMPLICIT NONE
426
427      !! * Arguments
428      TYPE(obs_surf), INTENT(INOUT) :: surfdata         ! Full set of surface data
429      TYPE(obswriinfo), OPTIONAL :: padd               ! Additional info for each variable
430      TYPE(obswriinfo), OPTIONAL :: pext               ! Extra info
431
432      !! * Local declarations
433      TYPE(obfbdata) :: fbdata
434      CHARACTER(LEN=40) :: clfname         ! netCDF filename
435      CHARACTER(LEN=10) :: clfiletype
436      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
437      CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable
438      CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable
439      CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable
440      INTEGER :: jo
441      INTEGER :: ja
442      INTEGER :: je
443      INTEGER :: jv
444      INTEGER :: iadd
445      INTEGER :: iext
446      INTEGER :: indx_std
447      INTEGER :: iadd_std
448      INTEGER :: iadd_clm     
449      INTEGER :: iadd_mdt
450
451      IF ( PRESENT( pext ) ) THEN
452         iext = pext%inum 
453      ELSE
454         iext = 0 
455      ENDIF 
456
457      ! Set up number of additional variables to be ouput:
458      ! Hx, CLM, STD, MDT...
459      IF ( PRESENT( padd ) ) THEN
460         iadd = padd%inum
461      ELSE
462         iadd = 0
463      ENDIF
464
465      iadd_std = 0
466      indx_std = -1
467      IF ( surfdata%nextra > 0 ) THEN
468         DO je = 1, surfdata%nextra
469           IF ( TRIM( surfdata%cext(je) ) == 'STD' ) THEN
470             iadd_std = 1
471             indx_std = je
472           ENDIF
473         END DO
474      ENDIF
475     
476      iadd_clm = 0 
477      IF ( surfdata%lclim ) iadd_clm = 1 
478
479      iadd_mdt = 0 
480      IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1 
481
482      CALL init_obfbdata( fbdata )
483
484      SELECT CASE ( TRIM(surfdata%cvars(1)) )
485      CASE('SLA')
486
487         ! SLA needs special treatment because of MDT, so is all done here
488         ! Other variables are done more generically
489         ! No climatology for SLA, MDT is our best estimate of that and is already output.
490
491         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
492            &                 1 + iadd_mdt + iadd_std + iadd, & 
493            &                 1 + iext, .TRUE. )
494
495         clfiletype = 'slafb'
496         fbdata%cname(1)      = surfdata%cvars(1)
497         fbdata%coblong(1)    = 'Sea level anomaly'
498         fbdata%cobunit(1)    = 'Metres'
499         fbdata%cextname(1)   = 'MDT'
500         fbdata%cextlong(1)   = 'Mean dynamic topography'
501         fbdata%cextunit(1)   = 'Metres'
502         DO je = 1, iext
503            fbdata%cextname(je) = pext%cdname(je)
504            fbdata%cextlong(je) = pext%cdlong(je,1)
505            fbdata%cextunit(je) = pext%cdunit(je,1)
506         END DO
507         fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
508         fbdata%caddunit(1,1) = 'Metres' 
509         fbdata%caddname(2)   = 'SSH'
510         fbdata%caddlong(2,1) = 'Model Sea surface height'
511         fbdata%caddunit(2,1) = 'Metres'
512         fbdata%cgrid(1)      = 'T'
513         DO ja = 1, iadd
514            fbdata%caddname(2+iadd_std+ja) = padd%cdname(ja)
515            fbdata%caddlong(2+iadd_std+ja,1) = padd%cdlong(ja,1)
516            fbdata%caddunit(2+iadd_std+ja,1) = padd%cdunit(ja,1)
517         END DO
518
519      CASE('SST')
520
521         clfiletype = 'sstfb'
522         cllongname = 'Sea surface temperature'
523         clunits    = 'Degree centigrade'
524         clgrid     = 'T'
525
526      CASE('ICECONC')
527
528         clfiletype = 'sicfb'
529         cllongname = 'Sea ice'
530         clunits    = 'Fraction'
531         clgrid     = 'T'
532
533      CASE('SIT') 
534
535         clfiletype    = 'sitfb' 
536         cllongname    = 'Sea ice thickness' 
537         clunits       = 'm' 
538         clgrid        = 'T' 
539
540      CASE('FBD') 
541
542         clfiletype = 'fbdfb' 
543         ! Change label from FBD ("freeboard") to SIT ("Sea Ice Thickness")
544         surfdata%cvars(1) = 'SIT' 
545         cllongname = 'Sea ice thickness' 
546         clunits    = 'm' 
547         clgrid     = 'T' 
548
549      CASE('SSS')
550
551         clfiletype = 'sssfb'
552         cllongname = 'Sea surface salinity'
553         clunits    = 'psu'
554         clgrid     = 'T'
555
556      CASE('UVEL') 
557
558         clfiletype    = 'ssvfb' 
559         cllongname    = 'Eastward sea surface velocity' 
560         clunits       = 'm s-1' 
561         clgrid        = 'U' 
562         cllongname    = 'Northward sea surface velocity' 
563         clunits       = 'm s-1' 
564         clgrid        = 'V' 
565
566      CASE('SLCHLTOT') 
567
568         clfiletype    = 'slchltotfb' 
569         cllongname    = 'Surface total log10(chlorophyll)' 
570         clunits       = 'log10(mg/m3)' 
571         clgrid        = 'T' 
572
573      CASE('SLCHLDIA') 
574
575         clfiletype    = 'slchldiafb' 
576         cllongname    = 'Surface diatom log10(chlorophyll)' 
577         clunits       = 'log10(mg/m3)' 
578         clgrid        = 'T' 
579
580      CASE('SLCHLNON') 
581
582         clfiletype    = 'slchlnonfb' 
583         cllongname    = 'Surface non-diatom log10(chlorophyll)' 
584         clunits       = 'log10(mg/m3)' 
585         clgrid        = 'T' 
586
587      CASE('SLCHLDIN') 
588
589         clfiletype    = 'slchldinfb' 
590         cllongname    = 'Surface dinoflagellate log10(chlorophyll)' 
591         clunits       = 'log10(mg/m3)' 
592         clgrid        = 'T' 
593
594      CASE('SLCHLMIC') 
595
596         clfiletype    = 'slchlmicfb' 
597         cllongname    = 'Surface microphytoplankton log10(chlorophyll)' 
598         clunits       = 'log10(mg/m3)' 
599         clgrid        = 'T' 
600
601      CASE('SLCHLNAN') 
602
603         clfiletype    = 'slchlnanfb' 
604         cllongname    = 'Surface nanophytoplankton log10(chlorophyll)' 
605         clunits       = 'log10(mg/m3)' 
606         clgrid        = 'T' 
607
608      CASE('SLCHLPIC') 
609
610         clfiletype    = 'slchlpicfb' 
611         cllongname    = 'Surface picophytoplankton log10(chlorophyll)' 
612         clunits       = 'log10(mg/m3)' 
613         clgrid        = 'T' 
614
615      CASE('SCHLTOT') 
616
617         clfiletype = 'schltotfb' 
618         cllongname = 'Surface total chlorophyll' 
619         clunits    = 'mg/m3' 
620         clgrid     = 'T' 
621
622      CASE('SLPHYTOT') 
623
624         clfiletype    = 'slphytotfb' 
625         cllongname    = 'Surface total log10(phytoplankton carbon)' 
626         clunits       = 'log10(mmolC/m3)' 
627         clgrid        = 'T' 
628
629      CASE('SLPHYDIA') 
630
631         clfiletype    = 'slphydiafb' 
632         cllongname    = 'Surface diatom log10(phytoplankton carbon)' 
633         clunits       = 'log10(mmolC/m3)' 
634         clgrid        = 'T' 
635
636      CASE('SLPHYNON') 
637
638         clfiletype    = 'slphynonfb' 
639         cllongname    = 'Surface non-diatom log10(phytoplankton carbon)' 
640         clunits       = 'log10(mmolC/m3)' 
641         clgrid        = 'T' 
642
643      CASE('SSPM') 
644
645         clfiletype    = 'sspmfb' 
646         cllongname    = 'Surface suspended particulate matter' 
647         clunits       = 'g/m3' 
648         clgrid        = 'T' 
649
650      CASE('SKD490') 
651
652         clfiletype    = 'skd490fb' 
653         cllongname    = 'Surface attenuation coefficient of downwelling radiation at 490 nm' 
654         clunits       = 'm-1' 
655         clgrid        = 'T' 
656
657      CASE('SFCO2') 
658
659         clfiletype    = 'sfco2fb' 
660         cllongname    = 'Surface fugacity of carbon dioxide' 
661         clunits       = 'uatm' 
662         clgrid        = 'T' 
663
664      CASE('SPCO2') 
665
666         clfiletype    = 'spco2fb' 
667         cllongname    = 'Surface partial pressure of carbon dioxide' 
668         clunits       = 'uatm' 
669         clgrid        = 'T'
670
671      CASE DEFAULT
672
673         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
674
675      END SELECT
676
677      ! SLA needs special treatment because of MDT, so is done above
678      ! Remaining variables treated more generically
679
680      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN
681
682         CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, & 
683            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. ) 
684
685         DO jv = 1, surfdata%nvar 
686            fbdata%cname(jv)      = surfdata%cvars(jv) 
687            fbdata%coblong(jv)    = cllongname 
688            fbdata%cobunit(jv)    = clunits
689         END DO
690         DO je = 1, iext 
691            fbdata%cextname(je) = pext%cdname(je) 
692            fbdata%cextlong(je) = pext%cdlong(je,1) 
693            fbdata%cextunit(je) = pext%cdunit(je,1) 
694         END DO
695         DO jv = 1, surfdata%nvar         
696            IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN
697               fbdata%caddlong(1,jv) = 'Model interpolated ICE' 
698            ELSE
699               fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv)) 
700            ENDIF
701            fbdata%caddunit(1,jv) = clunits
702            fbdata%cgrid(jv)      = clgrid
703         END DO             
704         DO ja = 1, iadd 
705            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja) 
706            DO jv = 1, surfdata%nvar                     
707               fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv) 
708               fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv) 
709            END DO
710         END DO
711      ENDIF
712     
713      fbdata%caddname(1)   = 'Hx'
714      IF ( indx_std /= -1 ) THEN
715         fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std) 
716         DO jv = 1, surfdata%nvar                         
717            fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation' 
718            fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1) 
719         END DO
720      ENDIF
721
722      IF ( surfdata%lclim ) THEN
723         fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm)   = 'CLM' 
724         DO jv = 1, surfdata%nvar                                 
725            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology' 
726            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1) 
727         END DO
728      ENDIF
729     
730      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
731
732      IF(lwp) THEN
733         WRITE(numout,*)
734         WRITE(numout,*)'obs_wri_surf :'
735         WRITE(numout,*)'~~~~~~~~~~~~~'
736         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
737      ENDIF
738
739      ! Transform surf data structure into obfbdata structure
740      fbdata%cdjuldref = '19500101000000'
741      DO jo = 1, surfdata%nsurf
742         fbdata%plam(jo)      = surfdata%rlam(jo)
743         fbdata%pphi(jo)      = surfdata%rphi(jo)
744         WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
745         fbdata%ivqc(jo,:)    = 0
746         fbdata%ivqcf(:,jo,:) = 0
747         IF ( surfdata%nqc(jo) > 255 ) THEN
748            fbdata%ioqc(jo)    = 4
749            fbdata%ioqcf(1,jo) = 0
750            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111')
751         ELSE
752            fbdata%ioqc(jo)    = surfdata%nqc(jo)
753            fbdata%ioqcf(:,jo) = 0
754         ENDIF
755         fbdata%ipqc(jo)      = 0
756         fbdata%ipqcf(:,jo)   = 0
757         fbdata%itqc(jo)      = 0
758         fbdata%itqcf(:,jo)   = 0
759         fbdata%cdwmo(jo)     = surfdata%cwmo(jo)
760         fbdata%kindex(jo)    = surfdata%nsfil(jo)
761         IF (ln_grid_global) THEN
762            DO jv = 1, surfdata%nvar 
763               fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv) 
764               fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv) 
765            END DO
766         ELSE
767            DO jv = 1, surfdata%nvar         
768               fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv)) 
769               fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv)) 
770            END DO
771         ENDIF
772         CALL greg2jul( 0, &
773            &           surfdata%nmin(jo), &
774            &           surfdata%nhou(jo), &
775            &           surfdata%nday(jo), &
776            &           surfdata%nmon(jo), &
777            &           surfdata%nyea(jo), &
778            &           fbdata%ptim(jo),   &
779            &           krefdate = 19500101 )
780         fbdata%pdep(1,jo)     = 0.0
781         fbdata%idqc(1,jo)     = 0
782         fbdata%idqcf(:,1,jo)  = 0
783         DO jv = 1, surfdata%nvar 
784            fbdata%pob(1,jo,jv)    = surfdata%robs(jo,jv) 
785
786            IF ( surfdata%nqc(jo) > 255 ) THEN
787               fbdata%ivqc(jo,jv)       = 4 
788               fbdata%ivlqc(1,jo,jv)    = 4 
789               fbdata%ivlqcf(1,1,jo,jv) = 0 
790               fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000000011111111') 
791            ELSE
792               fbdata%ivqc(jo,jv)       = surfdata%nqc(jo) 
793               fbdata%ivlqc(1,jo,jv)    = surfdata%nqc(jo) 
794               fbdata%ivlqcf(:,1,jo,jv) = 0 
795            ENDIF
796            fbdata%iobsk(1,jo,jv)  = 0 
797
798            ! Additional variables.
799            ! Hx is always the first additional variable
800            fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv) 
801            ! MDT is output as an additional variable if SLA obs type
802            IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN
803               fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1) 
804            ENDIF     
805            ! STD is output as an additional variable if available
806            IF ( indx_std /= -1 ) THEN
807               fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std) 
808            ENDIF 
809            ! CLM is output as an additional variable if available
810            IF ( surfdata%lclim ) THEN
811               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv) 
812            ENDIF 
813            ! Then other additional variables are output
814            DO ja = 1, iadd 
815               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = & 
816                  & surfdata%rext(jo,padd%ipoint(ja)) 
817            END DO
818         END DO         
819         ! Extra variables
820         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)         
821         DO je = 1, iext 
822            fbdata%pext(1,jo,1+je) = & 
823               & surfdata%rext(jo,pext%ipoint(je)) 
824         END DO
825      END DO
826
827      ! Write the obfbdata structure
828      CALL write_obfbdata( clfname, fbdata )
829
830      ! Output some basic statistics
831      CALL obs_wri_stats( fbdata )
832
833      CALL dealloc_obfbdata( fbdata )
834
835   END SUBROUTINE obs_wri_surf
836
837   SUBROUTINE obs_wri_stats( fbdata )
838      !!-----------------------------------------------------------------------
839      !!
840      !!                     *** ROUTINE obs_wri_stats  ***
841      !!
842      !! ** Purpose : Output some basic statistics of the data being written out
843      !!
844      !! ** Method  :
845      !!
846      !! ** Action  :
847      !!
848      !!      ! 2014-08  (D. Lea) Initial version
849      !!-----------------------------------------------------------------------
850
851      !! * Arguments
852      TYPE(obfbdata) :: fbdata
853
854      !! * Local declarations
855      INTEGER :: jvar
856      INTEGER :: jo
857      INTEGER :: jk
858      INTEGER :: inumgoodobs
859      INTEGER :: inumgoodobsmpp
860      REAL(wp) :: zsumx
861      REAL(wp) :: zsumx2
862      REAL(wp) :: zomb
863     
864
865      IF (lwp) THEN
866         WRITE(numout,*) ''
867         WRITE(numout,*) 'obs_wri_stats :'
868         WRITE(numout,*) '~~~~~~~~~~~~~~~'
869      ENDIF
870
871      DO jvar = 1, fbdata%nvar
872         zsumx=0.0_wp
873         zsumx2=0.0_wp
874         inumgoodobs=0
875         DO jo = 1, fbdata%nobs
876            DO jk = 1, fbdata%nlev
877               IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
878                  & ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
879                  & ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
880
881                  zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
882                  zsumx=zsumx+zomb
883                  zsumx2=zsumx2+zomb**2
884                  inumgoodobs=inumgoodobs+1
885               ENDIF
886            ENDDO
887         ENDDO
888
889         CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
890         CALL mpp_sum(zsumx)
891         CALL mpp_sum(zsumx2)
892
893         IF (lwp) THEN
894            WRITE(numout,*) 'Type: ',fbdata%cname(jvar),'  Total number of good observations: ',inumgoodobsmpp 
895            WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
896            WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
897            WRITE(numout,*) ''
898         ENDIF
899
900      ENDDO
901
902   END SUBROUTINE obs_wri_stats
903
904END MODULE obs_write
Note: See TracBrowser for help on using the repository browser.