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

source: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 11455

Last change on this file since 11455 was 11455, checked in by mattmartin, 5 years ago

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