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

source: branches/UKMO/dev_r5518_obs_oper_update_kd490/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 11253

Last change on this file since 11253 was 11253, checked in by dford, 5 years ago

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