source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_pfts/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8105

Last change on this file since 8105 was 8105, checked in by dford, 4 years ago

Implement initial version of PFT logchl obs operator, developed by Jozef S.

File size: 99.5 KB
Line 
1MODULE diaobs
2   !!======================================================================
3   !!                       ***  MODULE diaobs  ***
4   !! Observation diagnostics: Computation of the misfit between data and
5   !!                          their model equivalent
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   'key_diaobs' : Switch on the observation diagnostic computation
10   !!----------------------------------------------------------------------
11   !!   dia_obs_init : Reading and prepare observations
12   !!   dia_obs      : Compute model equivalent to observations
13   !!   dia_obs_wri  : Write observational diagnostics
14   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS
15   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
16   !!----------------------------------------------------------------------
17   !! * Modules used   
18   USE trc,         ONLY: trn            ! ERSEM state variables
19   USE wrk_nemo                 ! Memory Allocation
20   USE par_kind                 ! Precision variables
21   USE in_out_manager           ! I/O manager
22   USE par_oce
23   USE dom_oce                  ! Ocean space and time domain variables
24   USE obs_const, ONLY: obfillflt ! Fill value
25   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
26   USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
27   USE obs_read_sla             ! Reading and allocation of SLA observations 
28   USE obs_read_sst             ! Reading and allocation of SST observations 
29   USE obs_sstbias              ! Bias correction routine for SST
30   USE obs_readmdt              ! Reading and allocation of MDT for SLA.
31   USE obs_read_seaice          ! Reading and allocation of Sea Ice observations 
32   USE obs_read_vel             ! Reading and allocation of velocity component observations
33   USE obs_read_logchl          ! Reading and allocation of logchl observations
34   USE obs_read_spm             ! Reading and allocation of spm observations
35   USE obs_read_fco2            ! Reading and allocation of fco2 observations
36   USE obs_read_pco2            ! Reading and allocation of pco2 observations
37   USE obs_prep                 ! Preparation of obs. (grid search etc).
38   USE obs_oper                 ! Observation operators
39   USE obs_write                ! Writing of observation related diagnostics
40   USE obs_grid                 ! Grid searching
41   USE obs_read_altbias         ! Bias treatment for altimeter
42   USE obs_profiles_def         ! Profile data definitions
43   USE obs_profiles             ! Profile data storage
44   USE obs_surf_def             ! Surface data definitions
45   USE obs_sla                  ! SLA data storage
46   USE obs_sst                  ! SST data storage
47   USE obs_seaice               ! Sea Ice data storage
48   USE obs_logchl               ! logchl data storage
49   USE obs_spm                  ! spm data storage
50   USE obs_fco2                 ! fco2 data storage
51   USE obs_pco2                 ! pco2 data storage
52   USE obs_types                ! Definitions for observation types
53   USE mpp_map                  ! MPP mapping
54   USE lib_mpp                  ! For ctl_warn/stop
55
56   IMPLICIT NONE
57
58   !! * Routine accessibility
59   PRIVATE
60   PUBLIC dia_obs_init, &  ! Initialize and read observations
61      &   dia_obs,      &  ! Compute model equivalent to observations
62      &   dia_obs_wri,  &  ! Write model equivalent to observations
63      &   dia_obs_dealloc  ! Deallocate dia_obs data
64
65   !! * Shared Module variables
66   LOGICAL, PUBLIC, PARAMETER :: &
67#if defined key_diaobs
68      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
69#else
70      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
71#endif
72
73   !! * Module variables
74   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
75   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
76   LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set
77   LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set
78   LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles
79   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies
80   LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files
81   LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files
82   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
83   LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature
84   LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data
85   LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files
86   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
87   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
88   LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data
89   LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data
90   LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data
91   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data
92   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
93   LOGICAL, PUBLIC :: ln_logchl      !: Logical switch for log10(chlorophyll)
94   LOGICAL, PUBLIC :: ln_logchlfb    !: Logical switch for logchl from feedback files
95
96
97
98
99
100
101
102
103
104
105
106! JOZEF ADDITION
107
108
109   LOGICAL, PUBLIC :: ln_logchlpft      !: Logical switch for log10(chlorophyll) PFTs
110   LOGICAL, PUBLIC :: ln_logchlpftfb    !: Logical switch for logchl PFTs from feedback files
111
112   INTEGER, PUBLIC :: nn_logchlpftscc = 4 ! Number of logchl PFTs
113
114
115
116! END
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131   LOGICAL, PUBLIC :: ln_spm         !: Logical switch for spm
132   LOGICAL, PUBLIC :: ln_spmfb       !: Logical switch for spm from feedback files
133   LOGICAL, PUBLIC :: ln_fco2        !: Logical switch for fco2
134   LOGICAL, PUBLIC :: ln_fco2fb      !: Logical switch for fco2 from feedback files
135   LOGICAL, PUBLIC :: ln_pco2        !: Logical switch for pco2
136   LOGICAL, PUBLIC :: ln_pco2fb      !: Logical switch for pco2 from feedback files
137   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
138   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
139   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
140   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
141   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias 
142   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
143   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
144   LOGICAL, PUBLIC :: ln_sstbias     !: Logical switch for bias corection of SST
145
146   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
147   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
148 
149   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
150   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method
151
152   INTEGER, DIMENSION(imaxavtypes) :: &
153      & endailyavtypes !: ENACT data types which are daily average
154
155   INTEGER, PARAMETER :: MaxNumFiles = 1000
156   
157
158   
159 
160   LOGICAL, DIMENSION(MaxNumFiles) :: &
161      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
162   !                    !: If so use endailyavtypes
163      & ln_profb_enatim !: Change tim for 820 enact data set.
164   
165   INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type
166
167   LOGICAL, DIMENSION(MaxNumFiles) :: &
168      & ln_velfb_av   !: Is the velocity feedback files daily average?
169   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
170      & ld_enact     !: Profile data is ENACT so use endailyavtypes
171   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
172      & ld_velav     !: Velocity data is daily averaged
173   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
174      & ld_sstnight  !: SST observation corresponds to night mean
175
176   !!----------------------------------------------------------------------
177   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
178   !! $Id$
179   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
180   !!----------------------------------------------------------------------
181
182   !! * Substitutions
183#  include "domzgr_substitute.h90"
184CONTAINS
185
186   SUBROUTINE dia_obs_init
187      !!----------------------------------------------------------------------
188      !!                    ***  ROUTINE dia_obs_init  ***
189      !!         
190      !! ** Purpose : Initialize and read observations
191      !!
192      !! ** Method  : Read the namelist and call reading routines
193      !!
194      !! ** Action  : Read the namelist and call reading routines
195      !!
196      !! History :
197      !!        !  06-03  (K. Mogensen) Original code
198      !!        !  06-05  (A. Weaver) Reformatted
199      !!        !  06-10  (A. Weaver) Cleaning and add controls
200      !!        !  07-03  (K. Mogensen) General handling of profiles
201      !!        !  14-08  (J.While) Incorporated SST bias correction
202      !!----------------------------------------------------------------------
203
204      IMPLICIT NONE
205
206     
207      !! * Local declarations
208      CHARACTER(len=128) :: enactfiles(MaxNumFiles)
209      CHARACTER(len=128) :: coriofiles(MaxNumFiles)
210      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
211      CHARACTER(len=128) :: sstfiles(MaxNumFiles)     
212      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
213      CHARACTER(len=128) :: sstbias_files(MaxNumFiles) 
214      CHARACTER(len=128) :: slafilesact(MaxNumFiles)     
215      CHARACTER(len=128) :: slafilespas(MaxNumFiles)     
216      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
217      CHARACTER(len=128) :: seaicefiles(MaxNumFiles)           
218      CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 
219      CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)   
220      CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
221      CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
222      CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
223      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
224      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
225      CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
226      CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
227
228
229
230
231
232! JOZEF ADDITION
233
234
235
236      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: logchlpftfiles
237      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: logchlpftfbfiles
238
239 !      CHARACTER(len=128) :: logchlpftfiles(MaxNumFiles, nn_logchlpftscc)
240 !      CHARACTER(len=128) :: logchlpftfbfiles(MaxNumFiles, nn_logchlpftscc)
241
242! END
243
244
245
246
247
248
249
250      CHARACTER(len=128) :: spmfiles(MaxNumFiles)
251      CHARACTER(len=128) :: spmfbfiles(MaxNumFiles)
252      CHARACTER(len=128) :: fco2files(MaxNumFiles)
253      CHARACTER(len=128) :: fco2fbfiles(MaxNumFiles)
254      CHARACTER(len=128) :: pco2files(MaxNumFiles)
255      CHARACTER(len=128) :: pco2fbfiles(MaxNumFiles)
256      CHARACTER(LEN=128) :: reysstname
257      CHARACTER(LEN=12)  :: reysstfmt
258      CHARACTER(LEN=128) :: bias_file
259      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
260      NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       &
261         &            ln_sla, ln_sladt, ln_slafb,                     &
262         &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       &
263         &            ln_bound_reject,                                &
264         &            enactfiles, coriofiles, profbfiles,             &
265         &            slafilesact, slafilespas, slafbfiles,           &
266         &            sstfiles, sstfbfiles,                           &
267         &            ln_seaice, seaicefiles,                         &
268         &            dobsini, dobsend, n1dint, n2dint,               &
269         &            nmsshc, mdtcorr, mdtcutoff,                     &
270         &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    &
271         &            ln_sstnight,                                    &
272         &            ln_grid_search_lookup,                          &
273         &            grid_search_file, grid_search_res,              &
274         &            ln_grid_global, bias_file, ln_altbias,          &
275         &            endailyavtypes, ln_s_at_t, ln_profb_ena,        &
276         &            ln_vel3d, ln_velavcur, velavcurfiles,           &
277         &            ln_velhrcur, velhrcurfiles,                     &
278         &            ln_velavadcp, velavadcpfiles,                   &
279         &            ln_velhradcp, velhradcpfiles,                   &
280         &            ln_velfb, velfbfiles, ln_velfb_av,              &
281         &            ln_logchl, ln_logchlfb,                         &
282         &            logchlfiles, logchlfbfiles,                     &
283
284
285
286! JOZEF ADDITION
287
288         &            ln_logchlpft, ln_logchlpftfb,                     &
289         &          logchlpftfiles, logchlpftfbfiles,                   &
290
291! END
292
293
294
295
296         &            ln_spm, ln_spmfb,                               &
297         &            spmfiles, spmfbfiles,                           &
298         &            ln_fco2, ln_fco2fb,                             &
299         &            fco2files, fco2fbfiles,                         &
300         &            ln_pco2, ln_pco2fb,                             &
301         &            pco2files, pco2fbfiles,                         &
302         &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
303         &            ln_sstbias, sstbias_files
304
305      INTEGER :: jprofset
306      INTEGER :: jveloset
307      INTEGER :: jvar
308      INTEGER :: jnumenact
309      INTEGER :: jnumcorio
310      INTEGER :: jnumprofb
311      INTEGER :: jnumslaact
312      INTEGER :: jnumslapas
313      INTEGER :: jnumslafb
314      INTEGER :: jnumsst
315      INTEGER :: jnumsstfb
316      INTEGER :: jnumsstbias
317      INTEGER :: jnumseaice
318      INTEGER :: jnumvelavcur
319      INTEGER :: jnumvelhrcur 
320      INTEGER :: jnumvelavadcp
321      INTEGER :: jnumvelhradcp   
322      INTEGER :: jnumvelfb
323      INTEGER :: jnumlogchl
324      INTEGER :: jnumlogchlfb
325     
326
327     
328
329! JOZEF ADDITION
330 !     INTEGER :: nn_logchlpftscc
331
332  !    INTEGER, DIMENSION(:), ALLOCATABLE :: jnumlogchlpft
333  !    INTEGER, DIMENSION(:), ALLOCATABLE :: jnumlogchlpftfb
334
335     INTEGER :: jnumlogchlpft(nn_logchlpftscc)
336     INTEGER :: jnumlogchlpftfb(nn_logchlpftscc)
337
338
339! END
340
341
342
343
344
345      INTEGER :: jnumspm
346      INTEGER :: jnumspmfb
347      INTEGER :: jnumfco2
348      INTEGER :: jnumfco2fb
349      INTEGER :: jnumpco2
350      INTEGER :: jnumpco2fb
351      INTEGER :: ji
352      INTEGER :: jset
353      INTEGER :: ios                 ! Local integer output status for namelist read
354      LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
355
356
357      ! JOZEF ADDITION
358
359      INTEGER :: pft
360
361      !END
362
363
364
365      !-----------------------------------------------------------------------
366      ! Read namelist parameters
367      !-----------------------------------------------------------------------
368
369      ln_logchl   = .FALSE.
370      ln_logchlfb = .FALSE.
371
372
373
374
375
376
377   ! JOZEF ADDITION
378      ln_logchlpft   = .FALSE.
379      ln_logchlpftfb = .FALSE.
380
381  !    #if defined key_fabm
382  !     nn_logchlpftscc = 4
383  !    #elif defined key_medusa
384  !     nn_logchlpftscc = 2
385  !    #elif defined key_hadocc
386  !     nn_logchlpftscc = 1
387  !    #else
388  !     nn_logchlpftscc = 0
389  !    #fi
390   
391   ! END
392
393
394
395
396
397      ln_spm      = .FALSE.
398      ln_spmfb    = .FALSE.
399      ln_fco2     = .FALSE.
400      ln_fco2fb   = .FALSE.
401      ln_pco2     = .FALSE.
402      ln_pco2fb   = .FALSE.
403     
404      !Initalise all values in namelist arrays
405      enactfiles(:) = ''
406      coriofiles(:) = ''
407      profbfiles(:) = ''
408      slafilesact(:) = ''
409      slafilespas(:) = ''
410      slafbfiles(:) = ''
411      sstfiles(:)   = ''
412      sstfbfiles(:) = ''
413      seaicefiles(:) = ''
414      velcurfiles(:) = ''
415      veladcpfiles(:) = ''
416      velavcurfiles(:) = ''
417      velhrcurfiles(:) = ''
418      velavadcpfiles(:) = ''
419      velhradcpfiles(:) = ''
420      velfbfiles(:) = ''
421      velcurfiles(:) = ''
422      veladcpfiles(:) = ''
423      logchlfiles(:) = ''
424      logchlfbfiles(:) = ''
425
426
427
428
429! JOZEF ADDITION
430
431      ALLOCATE(logchlpftfiles(MaxNumFiles, nn_logchlpftscc))
432      ALLOCATE(logchlpftfbfiles(MaxNumFiles, nn_logchlpftscc))
433
434      DO pft = 1,nn_logchlpftscc
435      logchlpftfiles(:,pft) = ''
436      logchlpftfbfiles(:,pft) = ''
437      END DO
438! END
439
440
441
442
443
444
445
446      spmfiles(:) = ''
447      spmfbfiles(:) = ''
448      fco2files(:) = ''
449      fco2fbfiles(:) = ''
450      pco2files(:) = ''
451      pco2fbfiles(:) = ''
452      sstbias_files(:) = ''
453      endailyavtypes(:) = -1
454      endailyavtypes(1) = 820
455      ln_profb_ena(:) = .FALSE.
456      ln_profb_enatim(:) = .TRUE.
457      ln_velfb_av(:) = .FALSE.
458      ln_ignmis = .FALSE.
459      ln_bound_reject = .TRUE.
460
461      ! Read Namelist namobs : control observation diagnostics
462      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
463      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
464901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
465
466      CALL ini_date( dobsini )
467      CALL fin_date( dobsend )
468 
469      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
470      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
471902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
472      IF(lwm) WRITE ( numond, namobs )
473
474      ! Count number of files for each type
475      IF (ln_ena) THEN
476         lmask(:) = .FALSE.
477         WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
478         jnumenact = COUNT(lmask)
479      ENDIF
480      IF (ln_cor) THEN
481         lmask(:) = .FALSE.
482         WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
483         jnumcorio = COUNT(lmask)
484      ENDIF
485      IF (ln_profb) THEN
486         lmask(:) = .FALSE.
487         WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
488         jnumprofb = COUNT(lmask)
489      ENDIF
490      IF (ln_sladt) THEN
491         lmask(:) = .FALSE.
492         WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
493         jnumslaact = COUNT(lmask)
494         lmask(:) = .FALSE.
495         WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
496         jnumslapas = COUNT(lmask)
497      ENDIF
498      IF (ln_slafb) THEN
499         lmask(:) = .FALSE.
500         WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
501         jnumslafb = COUNT(lmask)
502         lmask(:) = .FALSE.
503      ENDIF
504      IF (ln_ghrsst) THEN
505         lmask(:) = .FALSE.
506         WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
507         jnumsst = COUNT(lmask)
508      ENDIF     
509      IF (ln_sstfb) THEN
510         lmask(:) = .FALSE.
511         WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
512         jnumsstfb = COUNT(lmask)
513         lmask(:) = .FALSE.
514      ENDIF
515      IF (ln_sstbias) THEN 
516         lmask(:) = .FALSE. 
517         WHERE (sstbias_files(:) /= '') lmask(:) = .TRUE. 
518         jnumsstbias = COUNT(lmask) 
519         lmask(:) = .FALSE. 
520      ENDIF       
521      IF (ln_seaice) THEN
522         lmask(:) = .FALSE.
523         WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
524         jnumseaice = COUNT(lmask)
525      ENDIF
526      IF (ln_velavcur) THEN
527         lmask(:) = .FALSE.
528         WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
529         jnumvelavcur = COUNT(lmask)
530      ENDIF
531      IF (ln_velhrcur) THEN
532         lmask(:) = .FALSE.
533         WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
534         jnumvelhrcur = COUNT(lmask)
535      ENDIF
536      IF (ln_velavadcp) THEN
537         lmask(:) = .FALSE.
538         WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
539         jnumvelavadcp = COUNT(lmask)
540      ENDIF
541      IF (ln_velhradcp) THEN
542         lmask(:) = .FALSE.
543         WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
544         jnumvelhradcp = COUNT(lmask)
545      ENDIF
546      IF (ln_velfb) THEN
547         lmask(:) = .FALSE.
548         WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
549         jnumvelfb = COUNT(lmask)
550         lmask(:) = .FALSE.
551      ENDIF
552      IF (ln_logchl) THEN
553         lmask(:) = .FALSE.
554         WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE.
555         jnumlogchl = COUNT(lmask)
556      ENDIF
557      IF (ln_logchlfb) THEN
558         lmask(:) = .FALSE.
559         WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE.
560         jnumlogchlfb = COUNT(lmask)
561      ENDIF
562
563
564
565
566
567
568
569
570
571
572! JOZEF ADDITION
573
574     DO pft=1,nn_logchlpftscc
575
576     IF (ln_logchlpft) THEN
577
578         lmask(:) = .FALSE.
579         WHERE (logchlpftfiles(:,pft) /= '') lmask(:) = .TRUE.
580         jnumlogchlpft(pft) = COUNT(lmask)
581
582      ENDIF
583      IF (ln_logchlpftfb) THEN
584         lmask(:) = .FALSE.
585         WHERE (logchlpftfbfiles(:,pft) /= '') lmask(:) = .TRUE.
586         jnumlogchlpftfb(pft) = COUNT(lmask)
587       
588
589      ENDIF
590
591      !WRITE(numout,*)'jnumlogchlpftfb is', jnumlogchlpftfbfiles(:,pft)
592
593    END DO
594! END
595
596
597
598WRITE(numout,*)'jnumlogchlpftfb is', logchlpftfbfiles
599
600
601
602
603
604      IF (ln_spm) THEN
605         lmask(:) = .FALSE.
606         WHERE (spmfiles(:) /= '') lmask(:) = .TRUE.
607         jnumspm = COUNT(lmask)
608      ENDIF
609      IF (ln_spmfb) THEN
610         lmask(:) = .FALSE.
611         WHERE (spmfbfiles(:) /= '') lmask(:) = .TRUE.
612         jnumspmfb = COUNT(lmask)
613      ENDIF
614      IF (ln_fco2) THEN
615         lmask(:) = .FALSE.
616         WHERE (fco2files(:) /= '') lmask(:) = .TRUE.
617         jnumfco2 = COUNT(lmask)
618      ENDIF
619      IF (ln_fco2fb) THEN
620         lmask(:) = .FALSE.
621         WHERE (fco2fbfiles(:) /= '') lmask(:) = .TRUE.
622         jnumfco2fb = COUNT(lmask)
623      ENDIF
624      IF (ln_pco2) THEN
625         lmask(:) = .FALSE.
626         WHERE (pco2files(:) /= '') lmask(:) = .TRUE.
627         jnumpco2 = COUNT(lmask)
628      ENDIF
629      IF (ln_pco2fb) THEN
630         lmask(:) = .FALSE.
631         WHERE (pco2fbfiles(:) /= '') lmask(:) = .TRUE.
632         jnumpco2fb = COUNT(lmask)
633      ENDIF
634     
635      ! Control print
636      IF(lwp) THEN
637         WRITE(numout,*)
638         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
639         WRITE(numout,*) '~~~~~~~~~~~~'
640         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
641         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
642         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
643         WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena
644         WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor
645         WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb
646         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
647         WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt
648         WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb
649         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
650         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
651         WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst   
652         WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst
653         WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb
654         WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias
655         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
656         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
657         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
658         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
659         WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur
660         WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur
661         WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp
662         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
663         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
664         WRITE(numout,*) '             Logical switch for logchl observations          ln_logchl = ', ln_logchl
665         WRITE(numout,*) '             Logical switch for feedback logchl data       ln_logchlfb = ', ln_logchlfb
666
667
668
669
670
671! JOZEF ADDITION
672
673         WRITE(numout,*) '             Logical switch for logchlpft observations          ln_logchlpft = ', ln_logchlpft
674         WRITE(numout,*) '             Logical switch for feedback logchlpft data       ln_logchlpftfb = ', ln_logchlpftfb
675
676! END
677
678
679
680
681
682         WRITE(numout,*) '             Logical switch for spm observations                ln_spm = ', ln_spm
683         WRITE(numout,*) '             Logical switch for feedback spm data             ln_spmfb = ', ln_spmfb
684         WRITE(numout,*) '             Logical switch for fco2 observations              ln_fco2 = ', ln_fco2
685         WRITE(numout,*) '             Logical switch for pco2 observations              ln_pco2 = ', ln_pco2
686         WRITE(numout,*) '             Logical switch for feedback pco2 data           ln_pco2fb = ', ln_pco2fb
687         WRITE(numout,*) '             Logical switch for feedback fco2 data           ln_fco2fb = ', ln_fco2fb
688         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
689         WRITE(numout,*) &
690   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
691         IF (ln_grid_search_lookup) &
692            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file
693         IF (ln_ena) THEN
694            DO ji = 1, jnumenact
695               WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', &
696                  TRIM(enactfiles(ji))
697            END DO
698         ENDIF
699         IF (ln_cor) THEN
700            DO ji = 1, jnumcorio
701               WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', &
702                  TRIM(coriofiles(ji))
703            END DO
704         ENDIF
705         IF (ln_profb) THEN
706            DO ji = 1, jnumprofb
707               IF (ln_profb_ena(ji)) THEN
708                  WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', &
709                     TRIM(profbfiles(ji))
710               ELSE
711                  WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', &
712                     TRIM(profbfiles(ji))
713               ENDIF
714               WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
715            END DO
716         ENDIF
717         IF (ln_sladt) THEN
718            DO ji = 1, jnumslaact
719               WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', &
720                  TRIM(slafilesact(ji))
721            END DO
722            DO ji = 1, jnumslapas
723               WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', &
724                  TRIM(slafilespas(ji))
725            END DO
726         ENDIF
727         IF (ln_slafb) THEN
728            DO ji = 1, jnumslafb
729               WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', &
730                  TRIM(slafbfiles(ji))
731            END DO
732         ENDIF
733         IF (ln_ghrsst) THEN
734            DO ji = 1, jnumsst
735               WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', &
736                  TRIM(sstfiles(ji))
737            END DO
738         ENDIF
739         IF (ln_sstfb) THEN
740            DO ji = 1, jnumsstfb
741               WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', &
742                  TRIM(sstfbfiles(ji))
743            END DO
744         ENDIF
745         IF (ln_seaice) THEN
746            DO ji = 1, jnumseaice
747               WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', &
748                  TRIM(seaicefiles(ji))
749            END DO
750         ENDIF
751         IF (ln_velavcur) THEN
752            DO ji = 1, jnumvelavcur
753               WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', &
754                  TRIM(velavcurfiles(ji))
755            END DO
756         ENDIF
757         IF (ln_velhrcur) THEN
758            DO ji = 1, jnumvelhrcur
759               WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', &
760                  TRIM(velhrcurfiles(ji))
761            END DO
762         ENDIF
763         IF (ln_velavadcp) THEN
764            DO ji = 1, jnumvelavadcp
765               WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', &
766                  TRIM(velavadcpfiles(ji))
767            END DO
768         ENDIF
769         IF (ln_velhradcp) THEN
770            DO ji = 1, jnumvelhradcp
771               WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', &
772                  TRIM(velhradcpfiles(ji))
773            END DO
774         ENDIF
775         IF (ln_velfb) THEN
776            DO ji = 1, jnumvelfb
777               IF (ln_velfb_av(ji)) THEN
778                  WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', &
779                     TRIM(velfbfiles(ji))
780               ELSE
781                  WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', &
782                     TRIM(velfbfiles(ji))
783               ENDIF
784            END DO
785         ENDIF
786
787
788         IF (ln_logchl) THEN
789            DO ji = 1, jnumlogchl
790               WRITE(numout,'(1X,2A)') '             logchl input observation file name        logchlfiles = ', &
791                  TRIM(logchlfiles(ji))
792            END DO
793         ENDIF
794         IF (ln_logchlfb) THEN
795            DO ji = 1, jnumlogchlfb
796               WRITE(numout,'(1X,2A)') '        Feedback logchl input observation file name  logchlfbfiles = ', &
797                  TRIM(logchlfbfiles(ji))
798            END DO
799         ENDIF
800
801
802
803   
804
805! JOZEF ADDITION
806
807          IF (ln_logchlpft) THEN
808         
809            DO pft=1,nn_logchlpftscc
810
811            DO ji = 1, jnumlogchlpft(pft)
812               WRITE(numout,'(1X,2A)') '             logchlpft input observation file name        logchlfiles = ', &
813                  TRIM(logchlpftfiles(ji,pft))
814            END DO
815            END DO
816         ENDIF
817         IF (ln_logchlpftfb) THEN
818
819            DO pft=1,nn_logchlpftscc
820
821            DO ji = 1, jnumlogchlpftfb(pft)
822               WRITE(numout,'(1X,2A)') '        Feedback logchlpft input observation file name  logchlfbfiles = ', &
823                  TRIM(logchlpftfbfiles(ji,pft))
824            END DO
825            END DO
826
827           
828         ENDIF
829
830
831!END
832
833
834
835
836
837
838
839
840         IF (ln_spm) THEN
841            DO ji = 1, jnumspm
842               WRITE(numout,'(1X,2A)') '             spm input observation file name              spmfiles = ', &
843                  TRIM(spmfiles(ji))
844            END DO
845         ENDIF
846         IF (ln_spmfb) THEN
847            DO ji = 1, jnumspmfb
848               WRITE(numout,'(1X,2A)') '             Feedback spm input observation file name   spmfbfiles = ', &
849                  TRIM(spmfbfiles(ji))
850            END DO
851         ENDIF
852         IF (ln_fco2) THEN
853            DO ji = 1, jnumfco2
854               WRITE(numout,'(1X,2A)') '             fco2 input observation file name            fco2files = ', &
855                  TRIM(fco2files(ji))
856            END DO
857         ENDIF
858         IF (ln_fco2fb) THEN
859            DO ji = 1, jnumfco2fb
860               WRITE(numout,'(1X,2A)') '            Feedback fco2 input observation file name  fco2fbfiles = ', &
861                  TRIM(fco2fbfiles(ji))
862            END DO
863         ENDIF
864         IF (ln_pco2) THEN
865            DO ji = 1, jnumpco2
866               WRITE(numout,'(1X,2A)') '             pco2 input observation file name            pco2files = ', &
867                  TRIM(pco2files(ji))
868            END DO
869         ENDIF
870         IF (ln_pco2fb) THEN
871            DO ji = 1, jnumpco2fb
872               WRITE(numout,'(1X,2A)') '            Feedback pco2 input observation file name  pco2fbfiles = ', &
873                  TRIM(pco2fbfiles(ji))
874            END DO
875         ENDIF
876         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini
877         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
878         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
879         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
880         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
881         WRITE(numout,*) '             Rejection of obs near open bdys       ln_bound_reject = ', ln_bound_reject
882         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
883         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
884         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
885         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
886         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
887         WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes
888
889      ENDIF
890     
891      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
892         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
893         RETURN
894      ENDIF
895
896      IF ( ln_grid_global ) THEN
897         CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
898      ENDIF
899
900      CALL obs_typ_init
901     
902      IF ( ln_grid_global ) THEN     
903         CALL mppmap_init
904      ENDIF
905     
906      ! Parameter control
907#if defined key_diaobs
908      IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
909         & ( .NOT. ln_vel3d ).AND.                                         &
910         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
911         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
912         & ( .NOT. ln_spm ).AND.( .NOT. ln_fco2 ).AND.( .NOT. ln_pco2 ) ) THEN
913         IF(lwp) WRITE(numout,cform_war)
914         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
915            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
916            &                    ' ln_logchl, ln_spm, ln_fco2, ln_pco2 are all set to .FALSE.'
917         nwarn = nwarn + 1
918      ENDIF
919#endif
920
921      CALL obs_grid_setup( )
922      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
923         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
924            &                    ' is not available')
925      ENDIF
926      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
927         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
928            &                    ' is not available')
929      ENDIF
930
931      !-----------------------------------------------------------------------
932      ! Depending on switches read the various observation types
933      !-----------------------------------------------------------------------
934      !  - Temperature/salinity profiles
935
936      IF ( ln_t3d .OR. ln_s3d ) THEN
937
938         ! Set the number of variables for profiles to 2 (T and S)
939         nprofvars = 2
940         ! Set the number of extra variables for profiles to 1 (insitu temp).
941         nprofextr = 1
942
943         ! Count how may insitu data sets we have and allocate data.
944         jprofset = 0
945         IF ( ln_ena ) jprofset = jprofset + 1
946         IF ( ln_cor ) jprofset = jprofset + 1
947         IF ( ln_profb ) jprofset = jprofset + jnumprofb
948         nprofsets = jprofset
949         IF ( nprofsets > 0 ) THEN
950            ALLOCATE(ld_enact(nprofsets))
951            ALLOCATE(profdata(nprofsets))
952            ALLOCATE(prodatqc(nprofsets))
953         ENDIF
954
955         jprofset = 0
956         
957         ! ENACT insitu data
958
959         IF ( ln_ena ) THEN
960
961            jprofset = jprofset + 1
962           
963            ld_enact(jprofset) = .TRUE.
964
965            CALL obs_rea_pro_dri( 1, profdata(jprofset),          &
966               &                  jnumenact, enactfiles(1:jnumenact), &
967               &                  nprofvars, nprofextr,        &
968               &                  nitend-nit000+2,             &
969               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
970               &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
971               &                  kdailyavtypes = endailyavtypes )
972
973            DO jvar = 1, 2
974
975               CALL obs_prof_staend( profdata(jprofset), jvar )
976
977            END DO
978
979            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
980               &              ln_t3d, ln_s3d, ln_nea, &
981               &              kdailyavtypes=endailyavtypes )
982           
983         ENDIF
984
985         ! Coriolis insitu data
986
987         IF ( ln_cor ) THEN
988           
989            jprofset = jprofset + 1
990
991            ld_enact(jprofset) = .FALSE.
992
993            CALL obs_rea_pro_dri( 2, profdata(jprofset),          &
994               &                  jnumcorio, coriofiles(1:jnumcorio), &
995               &                  nprofvars, nprofextr,        &
996               &                  nitend-nit000+2,             &
997               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
998               &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
999
1000            DO jvar = 1, 2
1001
1002               CALL obs_prof_staend( profdata(jprofset), jvar )
1003
1004            END DO
1005
1006            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
1007                 &            ln_t3d, ln_s3d, ln_nea )
1008           
1009         ENDIF
1010 
1011         ! Feedback insitu data
1012
1013         IF ( ln_profb ) THEN
1014           
1015            DO jset = 1, jnumprofb
1016               
1017               jprofset = jprofset + 1
1018               ld_enact (jprofset) = ln_profb_ena(jset)
1019
1020               CALL obs_rea_pro_dri( 0, profdata(jprofset),          &
1021                  &                  1, profbfiles(jset:jset), &
1022                  &                  nprofvars, nprofextr,        &
1023                  &                  nitend-nit000+2,             &
1024                  &                  dobsini, dobsend, ln_t3d, ln_s3d, &
1025                  &                  ln_ignmis, ln_s_at_t, &
1026                  &                  ld_enact(jprofset).AND.&
1027                  &                  ln_profb_enatim(jset), &
1028                  &                  .FALSE., kdailyavtypes = endailyavtypes )
1029               
1030               DO jvar = 1, 2
1031                 
1032                  CALL obs_prof_staend( profdata(jprofset), jvar )
1033                 
1034               END DO
1035               
1036               IF ( ld_enact(jprofset) ) THEN
1037                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
1038                     &              ln_t3d, ln_s3d, ln_nea, &
1039                     &              kdailyavtypes = endailyavtypes )
1040               ELSE
1041                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
1042                     &              ln_t3d, ln_s3d, ln_nea )
1043               ENDIF
1044               
1045            END DO
1046
1047         ENDIF
1048
1049      ENDIF
1050
1051      !  - Sea level anomalies
1052      IF ( ln_sla ) THEN
1053        ! Set the number of variables for sla to 1
1054         nslavars = 1
1055
1056         ! Set the number of extra variables for sla to 2
1057         nslaextr = 2
1058         
1059         ! Set the number of sla data sets to 2
1060         nslasets = 0
1061         IF ( ln_sladt ) THEN
1062            nslasets = nslasets + 2
1063         ENDIF
1064         IF ( ln_slafb ) THEN
1065            nslasets = nslasets + jnumslafb
1066         ENDIF
1067         
1068         ALLOCATE(sladata(nslasets))
1069         ALLOCATE(sladatqc(nslasets))
1070         sladata(:)%nsurf=0
1071         sladatqc(:)%nsurf=0
1072
1073         nslasets = 0
1074
1075         ! AVISO SLA data
1076
1077         IF ( ln_sladt ) THEN
1078
1079            ! Active SLA observations
1080           
1081            nslasets = nslasets + 1
1082           
1083            CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
1084               &              slafilesact(1:jnumslaact), &
1085               &              nslavars, nslaextr, nitend-nit000+2, &
1086               &              dobsini, dobsend, ln_ignmis, .FALSE. )
1087            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
1088               &              ln_sla, ln_nea )
1089           
1090            ! Passive SLA observations
1091           
1092            nslasets = nslasets + 1
1093           
1094            CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
1095               &              slafilespas(1:jnumslapas), &
1096               &              nslavars, nslaextr, nitend-nit000+2, &
1097               &              dobsini, dobsend, ln_ignmis, .FALSE. )
1098           
1099            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
1100               &              ln_sla, ln_nea )
1101
1102         ENDIF
1103         
1104         ! Feedback SLA data
1105
1106         IF ( ln_slafb ) THEN
1107
1108            DO jset = 1, jnumslafb
1109           
1110               nslasets = nslasets + 1
1111           
1112               CALL obs_rea_sla( 0, sladata(nslasets), 1, &
1113                  &              slafbfiles(jset:jset), &
1114                  &              nslavars, nslaextr, nitend-nit000+2, &
1115                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
1116               CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
1117                  &              ln_sla, ln_nea )
1118
1119            END DO               
1120
1121         ENDIF
1122         
1123         CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
1124           
1125         ! read in altimeter bias
1126         
1127         IF ( ln_altbias ) THEN     
1128            CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
1129         ENDIF
1130     
1131      ENDIF
1132
1133      !  - Sea surface height
1134      IF ( ln_ssh ) THEN
1135         IF(lwp) WRITE(numout,*) ' SSH currently not available'
1136      ENDIF
1137
1138      !  - Sea surface temperature
1139      IF ( ln_sst ) THEN
1140
1141         ! Set the number of variables for sst to 1
1142         nsstvars = 1
1143
1144         ! Set the number of extra variables for sst to 0
1145         nsstextr = 0
1146
1147         nsstsets = 0
1148
1149         IF (ln_reysst) nsstsets = nsstsets + 1
1150         IF (ln_ghrsst) nsstsets = nsstsets + 1
1151         IF ( ln_sstfb ) THEN
1152            nsstsets = nsstsets + jnumsstfb
1153         ENDIF
1154
1155         ALLOCATE(sstdata(nsstsets))
1156         ALLOCATE(sstdatqc(nsstsets))
1157         ALLOCATE(ld_sstnight(nsstsets))
1158         sstdata(:)%nsurf=0
1159         sstdatqc(:)%nsurf=0   
1160         ld_sstnight(:)=.false.
1161
1162         nsstsets = 0
1163
1164         IF (ln_reysst) THEN
1165
1166            nsstsets = nsstsets + 1
1167
1168            ld_sstnight(nsstsets) = ln_sstnight
1169
1170            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
1171               &                  nsstvars, nsstextr, &
1172               &                  nitend-nit000+2, dobsini, dobsend )
1173            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
1174               &              ln_nea )
1175
1176        ENDIF
1177       
1178        IF (ln_ghrsst) THEN
1179       
1180            nsstsets = nsstsets + 1
1181
1182            ld_sstnight(nsstsets) = ln_sstnight
1183         
1184            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
1185               &              sstfiles(1:jnumsst), &
1186               &              nsstvars, nsstextr, nitend-nit000+2, &
1187               &              dobsini, dobsend, ln_ignmis, .FALSE. )
1188            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
1189               &              ln_nea )
1190
1191        ENDIF
1192               
1193         ! Feedback SST data
1194
1195         IF ( ln_sstfb ) THEN
1196
1197            DO jset = 1, jnumsstfb
1198           
1199               nsstsets = nsstsets + 1
1200
1201               ld_sstnight(nsstsets) = ln_sstnight
1202           
1203               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
1204                  &              sstfbfiles(jset:jset), &
1205                  &              nsstvars, nsstextr, nitend-nit000+2, &
1206                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
1207               CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
1208                  &              ln_sst, ln_nea )
1209
1210            END DO               
1211
1212         ENDIF
1213         
1214         !Read in bias field and correct SST.
1215         IF ( ln_sstbias ) THEN
1216            IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
1217                                             "  but no bias"// & 
1218                                             " files to read in")   
1219            CALL obs_app_sstbias( nsstsets, sstdatqc, n2dint, & 
1220                                  jnumsstbias, & 
1221                                  sstbias_files(1:jnumsstbias) ) 
1222         ENDIF
1223
1224      ENDIF
1225
1226      !  - Sea surface salinity
1227      IF ( ln_sss ) THEN
1228         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1229      ENDIF
1230
1231      !  - Sea Ice Concentration
1232     
1233      IF ( ln_seaice ) THEN
1234
1235         ! Set the number of variables for seaice to 1
1236         nseaicevars = 1
1237
1238         ! Set the number of extra variables for seaice to 0
1239         nseaiceextr = 0
1240         
1241         ! Set the number of data sets to 1
1242         nseaicesets = 1
1243
1244         ALLOCATE(seaicedata(nseaicesets))
1245         ALLOCATE(seaicedatqc(nseaicesets))
1246         seaicedata(:)%nsurf=0
1247         seaicedatqc(:)%nsurf=0
1248
1249         CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
1250            &                 seaicefiles(1:jnumseaice), &
1251            &                 nseaicevars, nseaiceextr, nitend-nit000+2, &
1252            &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1253
1254         CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
1255            &                 ln_seaice, ln_nea )
1256 
1257      ENDIF
1258
1259      IF (ln_vel3d) THEN
1260
1261         ! Set the number of variables for profiles to 2 (U and V)
1262         nvelovars = 2
1263
1264         ! Set the number of extra variables for profiles to 2 to store
1265         ! rotation parameters
1266         nveloextr = 2
1267
1268         jveloset = 0
1269         
1270         IF ( ln_velavcur ) jveloset = jveloset + 1
1271         IF ( ln_velhrcur ) jveloset = jveloset + 1
1272         IF ( ln_velavadcp ) jveloset = jveloset + 1
1273         IF ( ln_velhradcp ) jveloset = jveloset + 1
1274         IF (ln_velfb) jveloset = jveloset + jnumvelfb
1275
1276         nvelosets = jveloset
1277         IF ( nvelosets > 0 ) THEN
1278            ALLOCATE( velodata(nvelosets) )
1279            ALLOCATE( veldatqc(nvelosets) )
1280            ALLOCATE( ld_velav(nvelosets) )
1281         ENDIF
1282         
1283         jveloset = 0
1284         
1285         ! Daily averaged data
1286
1287         IF ( ln_velavcur ) THEN
1288           
1289            jveloset = jveloset + 1
1290           
1291            ld_velav(jveloset) = .TRUE.
1292           
1293            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
1294               &                  velavcurfiles(1:jnumvelavcur), &
1295               &                  nvelovars, nveloextr, &
1296               &                  nitend-nit000+2,              &
1297               &                  dobsini, dobsend, ln_ignmis, &
1298               &                  ld_velav(jveloset), &
1299               &                  .FALSE. )
1300           
1301            DO jvar = 1, 2
1302               CALL obs_prof_staend( velodata(jveloset), jvar )
1303            END DO
1304           
1305            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1306               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1307           
1308         ENDIF
1309
1310         ! High frequency data
1311
1312         IF ( ln_velhrcur ) THEN
1313           
1314            jveloset = jveloset + 1
1315           
1316            ld_velav(jveloset) = .FALSE.
1317               
1318            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
1319               &                  velhrcurfiles(1:jnumvelhrcur), &
1320               &                  nvelovars, nveloextr, &
1321               &                  nitend-nit000+2,              &
1322               &                  dobsini, dobsend, ln_ignmis, &
1323               &                  ld_velav(jveloset), &
1324               &                  .FALSE. )
1325           
1326            DO jvar = 1, 2
1327               CALL obs_prof_staend( velodata(jveloset), jvar )
1328            END DO
1329           
1330            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1331               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1332           
1333         ENDIF
1334
1335         ! Daily averaged data
1336
1337         IF ( ln_velavadcp ) THEN
1338           
1339            jveloset = jveloset + 1
1340           
1341            ld_velav(jveloset) = .TRUE.
1342           
1343            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
1344               &                  velavadcpfiles(1:jnumvelavadcp), &
1345               &                  nvelovars, nveloextr, &
1346               &                  nitend-nit000+2,              &
1347               &                  dobsini, dobsend, ln_ignmis, &
1348               &                  ld_velav(jveloset), &
1349               &                  .FALSE. )
1350           
1351            DO jvar = 1, 2
1352               CALL obs_prof_staend( velodata(jveloset), jvar )
1353            END DO
1354           
1355            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1356               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1357           
1358         ENDIF
1359
1360         ! High frequency data
1361
1362         IF ( ln_velhradcp ) THEN
1363           
1364            jveloset = jveloset + 1
1365           
1366            ld_velav(jveloset) = .FALSE.
1367               
1368            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
1369               &                  velhradcpfiles(1:jnumvelhradcp), &
1370               &                  nvelovars, nveloextr, &
1371               &                  nitend-nit000+2,              &
1372               &                  dobsini, dobsend, ln_ignmis, &
1373               &                  ld_velav(jveloset), &
1374               &                  .FALSE. )
1375           
1376            DO jvar = 1, 2
1377               CALL obs_prof_staend( velodata(jveloset), jvar )
1378            END DO
1379           
1380            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1381               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1382           
1383         ENDIF
1384
1385         IF ( ln_velfb ) THEN
1386
1387            DO jset = 1, jnumvelfb
1388           
1389               jveloset = jveloset + 1
1390
1391               ld_velav(jveloset) = ln_velfb_av(jset)
1392               
1393               CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
1394                  &                  velfbfiles(jset:jset), &
1395                  &                  nvelovars, nveloextr, &
1396                  &                  nitend-nit000+2,              &
1397                  &                  dobsini, dobsend, ln_ignmis, &
1398                  &                  ld_velav(jveloset), &
1399                  &                  .FALSE. )
1400               
1401               DO jvar = 1, 2
1402                  CALL obs_prof_staend( velodata(jveloset), jvar )
1403               END DO
1404               
1405               CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1406                  &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1407
1408
1409            END DO
1410           
1411         ENDIF
1412
1413      ENDIF
1414
1415      !  - log10(chlorophyll)
1416     
1417      IF ( ln_logchl ) THEN
1418
1419         ! Set the number of variables for logchl to 1
1420         nlogchlvars = 1
1421
1422         ! Set the number of extra variables for logchl to 0
1423         nlogchlextr = 0
1424         
1425         IF ( ln_logchlfb ) THEN
1426            nlogchlsets = jnumlogchlfb
1427         ELSE
1428            nlogchlsets = 1
1429         ENDIF
1430
1431         ALLOCATE(logchldata(nlogchlsets))
1432         ALLOCATE(logchldatqc(nlogchlsets))
1433         logchldata(:)%nsurf=0
1434         logchldatqc(:)%nsurf=0
1435
1436         nlogchlsets = 0
1437
1438
1439
1440
1441
1442
1443         IF ( ln_logchlfb ) THEN             ! Feedback file format
1444
1445            DO jset = 1, jnumlogchlfb
1446           
1447               nlogchlsets = nlogchlsets + 1
1448
1449               CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, &
1450                  &                 logchlfbfiles(jset:jset), &
1451                  &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
1452                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1453
1454               CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1455                  &                 ln_logchl, ln_nea )
1456           
1457            ENDDO
1458
1459         ELSE                              ! Original file format
1460
1461            nlogchlsets = nlogchlsets + 1
1462
1463            CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, &
1464               &                 logchlfiles(1:jnumlogchl), &
1465               &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
1466               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1467
1468            CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1469               &                 ln_logchl, ln_nea )
1470
1471         ENDIF
1472 
1473         ENDIF
1474     
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484! JOZEF ADDITION
1485
1486
1487 IF ( ln_logchlpft ) THEN
1488
1489         ! Set the number of variables for logchl to 1
1490         nlogchlpftvars = 1   ! ??????? ASK DAVID
1491
1492         ! Set the number of extra variables for logchl to 0
1493         nlogchlpftextr = 0
1494         
1495         DO pft = 1, nn_logchlpftscc 
1496
1497
1498
1499         IF ( ln_logchlpftfb ) THEN
1500            nlogchlpftsets = jnumlogchlpftfb(pft)
1501         ELSE
1502            nlogchlpftsets = 1
1503         ENDIF
1504
1505
1506
1507         END DO
1508
1509         ALLOCATE(logchlpftdata(nn_logchlpftscc)) 
1510         ALLOCATE(logchlpftdatqc(nn_logchlpftscc))
1511         
1512
1513
1514         DO pft = 1, nn_logchlpftscc
1515
1516         logchlpftdata(pft)%nsurf=0
1517         logchlpftdatqc(pft)%nsurf=0
1518
1519         nlogchlpftsets = 0
1520
1521         WRITE(numout,*)'jnumlogchlpftfb is', jnumlogchlpftfb, ln_logchlpftfb
1522
1523         IF ( ln_logchlpftfb ) THEN             ! Feedback file format
1524
1525            DO jset = 1, jnumlogchlpftfb(pft)
1526           
1527               nlogchlpftsets = nlogchlpftsets + 1
1528
1529               CALL obs_rea_logchl( 0, logchlpftdata(pft), 1, &
1530                  &                 logchlpftfbfiles(jset:jset,pft), &
1531                  &                 nlogchlpftvars, nlogchlpftextr, nitend-nit000+2, &
1532                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1533
1534               CALL obs_pre_logchl( logchlpftdata(pft), logchlpftdatqc(pft), &
1535                  &                 ln_logchlpft, ln_nea )
1536           
1537            ENDDO
1538
1539         ELSE                              ! Original file format
1540
1541            nlogchlpftsets = 1
1542
1543            CALL obs_rea_logchl( 1, logchlpftdata(pft), jnumlogchlpft(pft), &
1544               &                 logchlpftfiles(1:jnumlogchlpft(pft),pft), &
1545               &                 nlogchlpftvars, nlogchlpftextr, nitend-nit000+2, &
1546               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1547
1548            CALL obs_pre_logchl( logchlpftdata(pft), logchlpftdatqc(pft), &
1549               &                 ln_logchlpft, ln_nea )
1550
1551         ENDIF
1552 
1553       ENDDO
1554
1555
1556
1557
1558
1559      ENDIF
1560
1561
1562nlogchlpftsets = 1   ! THIS SPECIAL CASE NEEDS TO BE GENERALIZED LATER
1563
1564WRITE(numout,*)'Producing nlogchlpftsets', nlogchlpftsets
1565
1566! END
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586      !  - spm
1587     
1588      IF ( ln_spm ) THEN
1589
1590         ! Set the number of variables for spm to 1
1591         nspmvars = 1
1592
1593         ! Set the number of extra variables for spm to 0
1594         nspmextr = 0
1595         
1596         IF ( ln_spmfb ) THEN
1597            nspmsets = jnumspmfb
1598         ELSE
1599            nspmsets = 1
1600         ENDIF
1601
1602         ALLOCATE(spmdata(nspmsets))
1603         ALLOCATE(spmdatqc(nspmsets))
1604         spmdata(:)%nsurf=0
1605         spmdatqc(:)%nsurf=0
1606
1607         nspmsets = 0
1608
1609         IF ( ln_spmfb ) THEN             ! Feedback file format
1610
1611            DO jset = 1, jnumspmfb
1612           
1613               nspmsets = nspmsets + 1
1614
1615               CALL obs_rea_spm( 0, spmdata(nspmsets), 1, &
1616                  &                 spmfbfiles(jset:jset), &
1617                  &                 nspmvars, nspmextr, nitend-nit000+2, &
1618                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1619
1620               CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1621                  &                 ln_spm, ln_nea )
1622           
1623            ENDDO
1624
1625         ELSE                              ! Original file format
1626
1627            nspmsets = nspmsets + 1
1628
1629            CALL obs_rea_spm( 1, spmdata(nspmsets), jnumspm, &
1630               &                 spmfiles(1:jnumspm), &
1631               &                 nspmvars, nspmextr, nitend-nit000+2, &
1632               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1633
1634            CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1635               &                 ln_spm, ln_nea )
1636
1637         ENDIF
1638 
1639      ENDIF
1640
1641      !  - fco2
1642     
1643      IF ( ln_fco2 ) THEN
1644
1645         ! Set the number of variables for fco2 to 1
1646         nfco2vars = 1
1647
1648         ! Set the number of extra variables for fco2 to 0
1649         nfco2extr = 0
1650         
1651         IF ( ln_fco2fb ) THEN
1652            nfco2sets = jnumfco2fb
1653         ELSE
1654            nfco2sets = 1
1655         ENDIF
1656
1657         ALLOCATE(fco2data(nfco2sets))
1658         ALLOCATE(fco2datqc(nfco2sets))
1659         fco2data(:)%nsurf=0
1660         fco2datqc(:)%nsurf=0
1661
1662         nfco2sets = 0
1663
1664         IF ( ln_fco2fb ) THEN             ! Feedback file format
1665
1666            DO jset = 1, jnumfco2fb
1667           
1668               nfco2sets = nfco2sets + 1
1669
1670               CALL obs_rea_fco2( 0, fco2data(nfco2sets), 1, &
1671                  &                 fco2fbfiles(jset:jset), &
1672                  &                 nfco2vars, nfco2extr, nitend-nit000+2, &
1673                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1674
1675               CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
1676                  &                 ln_fco2, ln_nea )
1677           
1678            ENDDO
1679
1680         ELSE                              ! Original file format
1681
1682            nfco2sets = nfco2sets + 1
1683
1684            CALL obs_rea_fco2( 1, fco2data(nfco2sets), jnumfco2, &
1685               &                 fco2files(1:jnumfco2), &
1686               &                 nfco2vars, nfco2extr, nitend-nit000+2, &
1687               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1688
1689            CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
1690               &                 ln_fco2, ln_nea )
1691
1692         ENDIF
1693 
1694      ENDIF
1695
1696      !  - pco2
1697     
1698      IF ( ln_pco2 ) THEN
1699
1700         ! Set the number of variables for pco2 to 1
1701         npco2vars = 1
1702
1703         ! Set the number of extra variables for pco2 to 0
1704         npco2extr = 0
1705         
1706         IF ( ln_pco2fb ) THEN
1707            npco2sets = jnumpco2fb
1708         ELSE
1709            npco2sets = 1
1710         ENDIF
1711
1712         ALLOCATE(pco2data(npco2sets))
1713         ALLOCATE(pco2datqc(npco2sets))
1714         pco2data(:)%nsurf=0
1715         pco2datqc(:)%nsurf=0
1716
1717         npco2sets = 0
1718
1719         IF ( ln_pco2fb ) THEN             ! Feedback file format
1720
1721            DO jset = 1, jnumpco2fb
1722           
1723               npco2sets = npco2sets + 1
1724
1725               CALL obs_rea_pco2( 0, pco2data(npco2sets), 1, &
1726                  &                 pco2fbfiles(jset:jset), &
1727                  &                 npco2vars, npco2extr, nitend-nit000+2, &
1728                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1729
1730               CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
1731                  &                 ln_pco2, ln_nea )
1732           
1733            ENDDO
1734
1735         ELSE                              ! Original file format
1736
1737            npco2sets = npco2sets + 1
1738
1739            CALL obs_rea_pco2( 1, pco2data(npco2sets), jnumpco2, &
1740               &                 pco2files(1:jnumpco2), &
1741               &                 npco2vars, npco2extr, nitend-nit000+2, &
1742               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1743
1744            CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
1745               &                 ln_pco2, ln_nea )
1746
1747         ENDIF
1748 
1749      ENDIF
1750     
1751   END SUBROUTINE dia_obs_init
1752
1753   SUBROUTINE dia_obs( kstp )
1754      !!----------------------------------------------------------------------
1755      !!                    ***  ROUTINE dia_obs  ***
1756      !!         
1757      !! ** Purpose : Call the observation operators on each time step
1758      !!
1759      !! ** Method  : Call the observation operators on each time step to
1760      !!              compute the model equivalent of the following date:
1761      !!               - T profiles
1762      !!               - S profiles
1763      !!               - Sea surface height (referenced to a mean)
1764      !!               - Sea surface temperature
1765      !!               - Sea surface salinity
1766      !!               - Velocity component (U,V) profiles
1767      !!               - Sea surface log10(chlorophyll)
1768      !!               - Sea surface spm
1769      !!               - Sea surface fco2
1770      !!               - Sea surface pco2
1771      !!
1772      !! ** Action  :
1773      !!
1774      !! History :
1775      !!        !  06-03  (K. Mogensen) Original code
1776      !!        !  06-05  (K. Mogensen) Reformatted
1777      !!        !  06-10  (A. Weaver) Cleaning
1778      !!        !  07-03  (K. Mogensen) General handling of profiles
1779      !!        !  07-04  (G. Smith) Generalized surface operators
1780      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
1781      !!        !  14-08  (J. While) observation operator for profiles in
1782      !!                             generalised vertical coordinates
1783      !!----------------------------------------------------------------------
1784      !! * Modules used
1785      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
1786         & rdt,           &                       
1787         & gdept_1d,       &             
1788#if defined key_vvl 
1789         & gdept_n,       &
1790#else 
1791         & gdept_1d,      &
1792#endif                                       
1793         & tmask, umask, vmask                           
1794      USE phycst, ONLY : &              ! Physical constants
1795         & rday                         
1796      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
1797         & tsn,  &             
1798         & un, vn,  &
1799         & sshn
1800#if defined  key_lim3
1801      USE ice, ONLY : &                     ! LIM Ice model variables
1802         & frld
1803#endif
1804#if defined key_lim2
1805      USE ice_2, ONLY : &                     ! LIM Ice model variables
1806         & frld
1807#endif
1808#if defined key_hadocc
1809      USE trc, ONLY :  &                ! HadOCC chlorophyll, fCO2 and pCO2
1810         & HADOCC_CHL, &
1811         & HADOCC_FCO2, &
1812         & HADOCC_PCO2, &
1813         & HADOCC_FILL_FLT
1814#elif defined key_medusa && defined key_foam_medusa
1815      USE trc, ONLY :  &                ! MEDUSA chlorophyll, fCO2 and pCO2
1816         & MEDUSA_CHL, &
1817         & MEDUSA_FCO2, &
1818         & MEDUSA_PCO2, &
1819         & MEDUSA_FILL_FLT
1820#elif defined key_fabm
1821      USE fabm
1822      USE par_fabm
1823#endif
1824#if defined key_spm
1825      USE par_spm, ONLY: &              ! ERSEM/SPM sediments
1826         & jp_spm
1827#endif
1828      IMPLICIT NONE
1829
1830      !! * Arguments
1831      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
1832      !! * Local declarations
1833      INTEGER :: idaystp                ! Number of timesteps per day
1834      INTEGER :: jprofset               ! Profile data set loop variable
1835      INTEGER :: jslaset                ! SLA data set loop variable
1836      INTEGER :: jsstset                ! SST data set loop variable
1837      INTEGER :: jseaiceset             ! sea ice data set loop variable
1838      INTEGER :: jveloset               ! velocity profile data loop variable
1839      INTEGER :: jlogchlset             ! logchl data set loop variable
1840      INTEGER :: jspmset                ! spm data set loop variable
1841      INTEGER :: jfco2set               ! fco2 data set loop variable
1842      INTEGER :: jpco2set               ! pco2 data set loop variable
1843      INTEGER :: jvar                   ! Variable number   
1844#if ! defined key_lim2 && ! defined key_lim3
1845      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
1846#endif
1847      REAL(wp) :: tiny                  ! small number
1848      REAL(wp), DIMENSION(jpi,jpj) :: &
1849         logchl                         ! array for log chlorophyll
1850
1851
1852
1853
1854
1855
1856
1857! JOZEF ADDITION
1858
1859
1860      REAL(wp), DIMENSION(jpi,jpj,nn_logchlpftscc) :: &
1861         logchlpfts
1862
1863
1864! END
1865
1866
1867
1868
1869
1870
1871
1872                         ! array for log chlorophyll pfts
1873      REAL(wp), DIMENSION(jpi,jpj) :: &
1874         maskchl                        ! array for special chlorophyll mask
1875      REAL(wp), DIMENSION(jpi,jpj) :: &
1876         spm                            ! array for spm
1877      REAL(wp), DIMENSION(jpi,jpj) :: &
1878         fco2                           ! array for fco2
1879      REAL(wp), DIMENSION(jpi,jpj) :: &
1880         maskfco2                       ! array for special fco2 mask
1881      REAL(wp), DIMENSION(jpi,jpj) :: &
1882         pco2                           ! array for pco2
1883      REAL(wp), DIMENSION(jpi,jpj) :: &
1884         maskpco2                       ! array for special pco2 mask
1885      INTEGER :: jn                     ! loop index
1886#if defined key_fabm
1887      REAL(wp), DIMENSION(jpi,jpj,jpk) :: logchl_3d
1888
1889
1890
1891
1892
1893
1894     
1895      ! JOZEF ADDITION
1896
1897      REAL(wp), DIMENSION(jpi,jpj,jpk,nn_logchlpftscc) :: logchlpfts_3d
1898#endif
1899
1900      INTEGER :: pft
1901      ! END
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1912 
1913#if ! defined key_lim2 && ! defined key_lim3
1914      CALL wrk_alloc(jpi,jpj,frld) 
1915#endif
1916
1917      IF(lwp) THEN
1918         WRITE(numout,*)
1919         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1920         WRITE(numout,*) '~~~~~~~'
1921      ENDIF
1922
1923      idaystp = NINT( rday / rdt )
1924
1925      !-----------------------------------------------------------------------
1926      ! No LIM => frld == 0.0_wp
1927      !-----------------------------------------------------------------------
1928#if ! defined key_lim2 && ! defined key_lim3
1929      frld(:,:) = 0.0_wp
1930#endif
1931      !-----------------------------------------------------------------------
1932      ! Depending on switches call various observation operators
1933      !-----------------------------------------------------------------------
1934
1935      !  - Temperature/salinity profiles
1936      IF ( ln_t3d .OR. ln_s3d ) THEN
1937         DO jprofset = 1, nprofsets
1938            IF( (.NOT. lk_vvl) .AND. (ln_zco .OR. ln_zps) ) THEN
1939               IF(lwp) THEN
1940                  WRITE(numout,*) 'dia_obs : calling obs_pro_opt'
1941               ENDIF
1942               IF ( ld_enact(jprofset) ) THEN
1943                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1944                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1945                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1946                     &              gdept_1d, tmask, n1dint, n2dint,        & 
1947                     &              kdailyavtypes = endailyavtypes ) 
1948               ELSE
1949                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1950                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1951                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1952                     &              gdept_1d, tmask, n1dint, n2dint               ) 
1953               ENDIF
1954            ELSE
1955               IF(lwp) THEN
1956                  WRITE(numout,*) 'dia_obs : calling obs_pro_sco_opt'
1957               ENDIF
1958               IF ( ld_enact(jprofset) ) THEN
1959                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1960                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1961                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1962                     &              fsdept(:,:,:), fsdepw(:,:,:),           & 
1963                     &              tmask, n1dint, n2dint,                  & 
1964                     &              kdailyavtypes = endailyavtypes ) 
1965               ELSE
1966                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1967                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1968                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1969                     &              fsdept(:,:,:), fsdepw(:,:,:),           &
1970                     &              tmask, n1dint, n2dint ) 
1971               ENDIF
1972            ENDIF
1973         END DO
1974      ENDIF
1975
1976      !  - Sea surface anomaly
1977      IF ( ln_sla ) THEN
1978         DO jslaset = 1, nslasets
1979            CALL obs_sla_opt( sladatqc(jslaset),            &
1980               &              kstp, jpi, jpj, nit000, sshn, &
1981               &              tmask(:,:,1), n2dint )
1982         END DO         
1983      ENDIF
1984
1985      !  - Sea surface temperature
1986      IF ( ln_sst ) THEN
1987         DO jsstset = 1, nsstsets
1988            CALL obs_sst_opt( sstdatqc(jsstset),                &
1989               &              kstp, jpi, jpj, nit000, idaystp,  &
1990               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
1991               &              n2dint, ld_sstnight(jsstset) )
1992         END DO
1993      ENDIF
1994
1995      !  - Sea surface salinity
1996      IF ( ln_sss ) THEN
1997         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1998      ENDIF
1999
2000#if defined key_lim2 || defined key_lim3
2001      IF ( ln_seaice ) THEN
2002         DO jseaiceset = 1, nseaicesets
2003            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
2004               &              kstp, jpi, jpj, nit000, 1.-frld, &
2005               &              tmask(:,:,1), n2dint )
2006         END DO
2007      ENDIF     
2008#endif
2009
2010      !  - Velocity profiles
2011      IF ( ln_vel3d ) THEN
2012         DO jveloset = 1, nvelosets
2013           ! zonal component of velocity
2014           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
2015              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
2016                             n1dint, n2dint, ld_velav(jveloset) )
2017         END DO
2018      ENDIF
2019
2020      IF ( ln_logchl ) THEN
2021
2022#if defined key_hadocc
2023         logchl(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2024#elif defined key_medusa && defined key_foam_medusa
2025         logchl(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2026#elif defined key_fabm
2027         logchl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
2028         logchl(:,:) = logchl_3d(:,:,1)
2029#else
2030         CALL ctl_stop( ' Trying to run logchl observation operator', &
2031            &           ' but no biogeochemical model appears to have been defined' )
2032#endif
2033
2034         maskchl(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
2035
2036         ! Take the log10 where we can, otherwise exclude
2037         tiny = 1.0e-20
2038         WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
2039            logchl(:,:)  = LOG10(logchl(:,:))
2040         ELSEWHERE
2041            logchl(:,:)  = obfillflt
2042            maskchl(:,:) = 0
2043         END WHERE
2044
2045         DO jlogchlset = 1, nlogchlsets
2046             CALL obs_logchl_opt( logchldatqc(jlogchlset),             &
2047               &                  kstp, jpi, jpj, nit000, logchl(:,:), &
2048               &                  maskchl(:,:), n2dint )
2049         END DO         
2050      ENDIF 
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060! JOZEF ADDITION
2061
2062
2063
2064IF ( ln_logchlpft ) THEN
2065
2066#if defined key_hadocc
2067         logchl(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2068#elif defined key_medusa && defined key_foam_medusa
2069         logchl(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2070#elif defined key_fabm
2071         logchlpfts_3d(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1)
2072         logchlpfts(:,:,1) = logchlpfts_3d(:,:,1,1)
2073         logchlpfts_3d(:,:,:,2) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl2)
2074         logchlpfts(:,:,2) = logchlpfts_3d(:,:,1,2)
2075         logchlpfts_3d(:,:,:,3) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl3)
2076         logchlpfts(:,:,3) = logchlpfts_3d(:,:,1,3)
2077         logchlpfts_3d(:,:,:,4) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
2078         logchlpfts(:,:,4) = logchlpfts_3d(:,:,1,4)
2079#else
2080         CALL ctl_stop( ' Trying to run logchl observation operator', &
2081            &           ' but no biogeochemical model appears to have been defined' )
2082#endif
2083
2084                  ! create a special mask to exclude certain things
2085
2086         ! Take the log10 where we can, otherwise exclude
2087         tiny = 1.0e-20
2088
2089         DO pft=1,nn_logchlpftscc
2090
2091         maskchl(:,:) = tmask(:,:,1)
2092
2093         WRITE(numout,*) 'Running pft: ',pft
2094
2095         WHERE(logchlpfts(:,:,pft) > tiny .AND. logchlpfts(:,:,pft) /= obfillflt )
2096            logchlpfts(:,:,pft)  = LOG10(logchlpfts(:,:,pft))
2097         ELSEWHERE
2098            logchlpfts(:,:,pft)  = obfillflt
2099            maskchl(:,:) = 0
2100         END WHERE
2101
2102         WRITE(numout,*) 'Max value of Pft: ', MINVAL(logchlpfts(:,:,pft))
2103
2104         DO jlogchlset = 1, nlogchlpftsets
2105             WRITE(numout,*) 'nlogchlpftsets, nn_logchlpftscc, obfillflt = ', nlogchlpftsets, nn_logchlpftscc, obfillflt
2106             WRITE(numout,*) 'kstp, jpi, jpj, nit000, n2dint = ', kstp, jpi, jpj, nit000, n2dint
2107             !WRITE(numout,*) 'ALLOCATED(logchlpftdatqc) = ', ALLOCATED(logchlpftdatqc)
2108             WRITE(numout,*) 'SHAPE(logchlpftdatqc) = ', SHAPE(logchlpftdatqc)
2109             WRITE(numout,*) 'SIZE(logchlpftdatqc) = ', SIZE(logchlpftdatqc)
2110             !WRITE(numout,*) 'ALLOCATED(logchlpfts) = ', ALLOCATED(logchlpfts)
2111             WRITE(numout,*) 'SHAPE(logchlpfts) = ', SHAPE(logchlpfts)
2112             WRITE(numout,*) 'SIZE(logchlpfts) = ', SIZE(logchlpfts)
2113             !WRITE(numout,*) 'ALLOCATED(maskchl) = ', ALLOCATED(maskchl)
2114             WRITE(numout,*) 'SHAPE(maskchl) = ', SHAPE(maskchl)
2115             WRITE(numout,*) 'SIZE(maskchl) = ', SIZE(maskchl)
2116             CALL flush(numout)
2117             
2118             CALL obs_logchl_opt( logchlpftdatqc(pft),             &
2119               &                  kstp, jpi, jpj, nit000, logchlpfts(:,:,pft), &
2120               &                  maskchl(:,:), n2dint )
2121
2122             WRITE(numout,*) 'Max value of Pft: ', MINVAL(logchlpfts(:,:,pft))
2123         END DO         
2124     
2125         
2126         END DO
2127
2128        ENDIF
2129
2130! END
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147      IF ( ln_spm ) THEN
2148#if defined key_spm
2149         spm(:,:) = 0.0
2150         DO jn = 1, jp_spm
2151            spm(:,:) = spm(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
2152         END DO
2153#else
2154         CALL ctl_stop( ' Trying to run spm observation operator', &
2155            &           ' but no spm model appears to have been defined' )
2156#endif
2157
2158         DO jspmset = 1, nspmsets
2159             CALL obs_spm_opt( spmdatqc(jspmset),                &
2160               &               kstp, jpi, jpj, nit000, spm(:,:), &
2161               &               tmask(:,:,1), n2dint )
2162         END DO         
2163      ENDIF
2164
2165      IF ( ln_fco2 ) THEN
2166         maskfco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
2167#if defined key_hadocc
2168         fco2(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
2169         IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
2170            fco2(:,:) = obfillflt
2171            maskfco2(:,:) = 0
2172            CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
2173               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2174               &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
2175         ENDIF
2176#elif defined key_medusa && defined key_foam_medusa
2177         fco2(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA
2178         IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
2179            fco2(:,:) = obfillflt
2180            maskfco2(:,:) = 0
2181            CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
2182               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2183               &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
2184         ENDIF
2185#elif defined key_fabm
2186         ! First, get pCO2 from FABM
2187         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
2188         pco2(:,:) = pco2_3d(:,:,1)
2189         ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
2190         ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
2191         ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
2192         ! and
2193         ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
2194         ! Marine Chemistry, 2: 203-215.
2195         ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
2196         ! not explicitly included - atmospheric pressure is not necessarily available so this is
2197         ! the best assumption.
2198         ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
2199         ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
2200         ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
2201         ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
2202         fco2(:,:) = pco2(:,:) * EXP((-1636.75                                                                               + &
2203            &                         12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
2204            &                         0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
2205            &                         0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
2206            &                         2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
2207            &                        (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
2208#else
2209         CALL ctl_stop( ' Trying to run fco2 observation operator', &
2210            &           ' but no biogeochemical model appears to have been defined' )
2211#endif
2212
2213         DO jfco2set = 1, nfco2sets
2214             CALL obs_fco2_opt( fco2datqc(jfco2set),                      &
2215               &                kstp, jpi, jpj, nit000, fco2(:,:), &
2216               &                maskfco2(:,:), n2dint )
2217         END DO
2218      ENDIF
2219
2220      IF ( ln_pco2 ) THEN
2221         maskpco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
2222#if defined key_hadocc
2223         pco2(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
2224         IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
2225            pco2(:,:) = obfillflt
2226            maskpco2(:,:) = 0
2227            CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
2228               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2229               &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
2230         ENDIF
2231#elif defined key_medusa && defined key_foam_medusa
2232         pco2(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA
2233         IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
2234            pco2(:,:) = obfillflt
2235            maskpco2(:,:) = 0
2236            CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
2237               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2238               &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
2239         ENDIF
2240#elif defined key_fabm
2241         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
2242         pco2(:,:) = pco2_3d(:,:,1)
2243#else
2244         CALL ctl_stop( ' Trying to run pCO2 observation operator', &
2245            &           ' but no biogeochemical model appears to have been defined' )
2246#endif
2247
2248         DO jpco2set = 1, npco2sets
2249             CALL obs_pco2_opt( pco2datqc(jpco2set),                      &
2250               &                kstp, jpi, jpj, nit000, pco2(:,:), &
2251               &                maskpco2(:,:), n2dint )
2252         END DO
2253      ENDIF
2254
2255#if ! defined key_lim2 && ! defined key_lim3
2256      CALL wrk_dealloc(jpi,jpj,frld) 
2257#endif
2258
2259   END SUBROUTINE dia_obs
2260 
2261   SUBROUTINE dia_obs_wri 
2262      !!----------------------------------------------------------------------
2263      !!                    ***  ROUTINE dia_obs_wri  ***
2264      !!         
2265      !! ** Purpose : Call observation diagnostic output routines
2266      !!
2267      !! ** Method  : Call observation diagnostic output routines
2268      !!
2269      !! ** Action  :
2270      !!
2271      !! History :
2272      !!        !  06-03  (K. Mogensen) Original code
2273      !!        !  06-05  (K. Mogensen) Reformatted
2274      !!        !  06-10  (A. Weaver) Cleaning
2275      !!        !  07-03  (K. Mogensen) General handling of profiles
2276      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
2277      !!----------------------------------------------------------------------
2278      IMPLICIT NONE
2279
2280      !! * Local declarations
2281
2282      INTEGER :: jprofset                 ! Profile data set loop variable
2283      INTEGER :: jveloset                 ! Velocity data set loop variable
2284      INTEGER :: jslaset                  ! SLA data set loop variable
2285      INTEGER :: jsstset                  ! SST data set loop variable
2286      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
2287      INTEGER :: jlogchlset               ! logchl data set loop variable
2288      INTEGER :: jspmset                  ! spm data set loop variable
2289      INTEGER :: jfco2set                 ! fco2 data set loop variable
2290      INTEGER :: jpco2set                 ! pco2 data set loop variable
2291      INTEGER :: jset
2292      INTEGER :: jfbini
2293      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
2294      CHARACTER(LEN=20) :: cdtmp
2295      INTEGER :: pft
2296      !-----------------------------------------------------------------------
2297      ! Depending on switches call various observation output routines
2298      !-----------------------------------------------------------------------
2299
2300      !  - Temperature/salinity profiles
2301
2302      IF( ln_t3d .OR. ln_s3d ) THEN
2303
2304         ! Copy data from prodatqc to profdata structures
2305         DO jprofset = 1, nprofsets
2306
2307            CALL obs_prof_decompress( prodatqc(jprofset), &
2308                 &                    profdata(jprofset), .TRUE., numout )
2309
2310         END DO
2311
2312         ! Write the profiles.
2313
2314         jprofset = 0
2315
2316         ! ENACT insitu data
2317
2318         IF ( ln_ena ) THEN
2319           
2320            jprofset = jprofset + 1
2321
2322            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
2323
2324         ENDIF
2325
2326         ! Coriolis insitu data
2327
2328         IF ( ln_cor ) THEN
2329           
2330            jprofset = jprofset + 1
2331
2332            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
2333           
2334         ENDIF
2335         
2336         ! Feedback insitu data
2337
2338         IF ( ln_profb ) THEN
2339
2340            jfbini = jprofset + 1
2341
2342            DO jprofset = jfbini, nprofsets
2343               
2344               jset = jprofset - jfbini + 1
2345               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
2346               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
2347
2348            END DO
2349
2350         ENDIF
2351
2352      ENDIF
2353
2354      !  - Sea surface anomaly
2355      IF ( ln_sla ) THEN
2356
2357         ! Copy data from sladatqc to sladata structures
2358         DO jslaset = 1, nslasets
2359
2360              CALL obs_surf_decompress( sladatqc(jslaset), &
2361                 &                    sladata(jslaset), .TRUE., numout )
2362
2363         END DO
2364
2365         jslaset = 0 
2366
2367         ! Write the AVISO SLA data
2368
2369         IF ( ln_sladt ) THEN
2370           
2371            jslaset = 1
2372            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
2373            jslaset = 2
2374            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
2375
2376         ENDIF
2377
2378         IF ( ln_slafb ) THEN
2379           
2380            jfbini = jslaset + 1
2381
2382            DO jslaset = jfbini, nslasets
2383               
2384               jset = jslaset - jfbini + 1
2385               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
2386               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
2387
2388            END DO
2389
2390         ENDIF
2391
2392      ENDIF
2393
2394      !  - Sea surface temperature
2395      IF ( ln_sst ) THEN
2396
2397         ! Copy data from sstdatqc to sstdata structures
2398         DO jsstset = 1, nsstsets
2399     
2400              CALL obs_surf_decompress( sstdatqc(jsstset), &
2401                 &                    sstdata(jsstset), .TRUE., numout )
2402
2403         END DO
2404
2405         jsstset = 0 
2406
2407         ! Write the AVISO SST data
2408
2409         IF ( ln_reysst ) THEN
2410           
2411            jsstset = jsstset + 1
2412            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
2413
2414         ENDIF
2415
2416         IF ( ln_ghrsst ) THEN
2417           
2418            jsstset = jsstset + 1
2419            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
2420
2421         ENDIF
2422
2423         IF ( ln_sstfb ) THEN
2424           
2425            jfbini = jsstset + 1
2426
2427            DO jsstset = jfbini, nsstsets
2428               
2429               jset = jsstset - jfbini + 1
2430               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
2431               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
2432
2433            END DO
2434
2435         ENDIF
2436
2437      ENDIF
2438
2439      !  - Sea surface salinity
2440      IF ( ln_sss ) THEN
2441         IF(lwp) WRITE(numout,*) ' SSS currently not available'
2442      ENDIF
2443
2444      !  - Sea Ice Concentration
2445      IF ( ln_seaice ) THEN
2446
2447         ! Copy data from seaicedatqc to seaicedata structures
2448         DO jseaiceset = 1, nseaicesets
2449
2450              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
2451                 &                    seaicedata(jseaiceset), .TRUE., numout )
2452
2453         END DO
2454
2455         ! Write the Sea Ice data
2456         DO jseaiceset = 1, nseaicesets
2457     
2458            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
2459            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
2460
2461         END DO
2462
2463      ENDIF
2464     
2465      ! Velocity data
2466      IF( ln_vel3d ) THEN
2467
2468         ! Copy data from veldatqc to velodata structures
2469         DO jveloset = 1, nvelosets
2470
2471            CALL obs_prof_decompress( veldatqc(jveloset), &
2472                 &                    velodata(jveloset), .TRUE., numout )
2473
2474         END DO
2475
2476         ! Write the profiles.
2477
2478         jveloset = 0
2479
2480         ! Daily averaged data
2481
2482         IF ( ln_velavcur ) THEN
2483           
2484            jveloset = jveloset + 1
2485
2486            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
2487
2488         ENDIF
2489
2490         ! High frequency data
2491
2492         IF ( ln_velhrcur ) THEN
2493           
2494            jveloset = jveloset + 1
2495
2496            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
2497
2498         ENDIF
2499
2500         ! Daily averaged data
2501
2502         IF ( ln_velavadcp ) THEN
2503           
2504            jveloset = jveloset + 1
2505
2506            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
2507
2508         ENDIF
2509
2510         ! High frequency data
2511
2512         IF ( ln_velhradcp ) THEN
2513           
2514            jveloset = jveloset + 1
2515           
2516            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
2517               
2518         ENDIF
2519
2520         ! Feedback velocity data
2521
2522         IF ( ln_velfb ) THEN
2523
2524            jfbini = jveloset + 1
2525
2526            DO jveloset = jfbini, nvelosets
2527               
2528               jset = jveloset - jfbini + 1
2529               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
2530               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
2531
2532            END DO
2533
2534         ENDIF
2535         
2536      ENDIF
2537
2538      !  - log10(chlorophyll)
2539      IF ( ln_logchl ) THEN
2540
2541         ! Copy data from logchldatqc to logchldata structures
2542         DO jlogchlset = 1, nlogchlsets
2543
2544            CALL obs_surf_decompress( logchldatqc(jlogchlset), &
2545                 &                    logchldata(jlogchlset), .TRUE., numout )
2546
2547         END DO
2548         
2549         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2550         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2551         DO jlogchlset = 1, nlogchlsets
2552            WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
2553               logchldata(jlogchlset)%nqc(:)    = 1
2554               logchldata(jlogchlset)%robs(:,1) = obfillflt
2555            END WHERE
2556         END DO
2557
2558         ! Write the logchl data
2559         DO jlogchlset = 1, nlogchlsets
2560     
2561            WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
2562            CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
2563
2564         END DO
2565
2566      ENDIF
2567
2568
2569
2570
2571
2572
2573
2574     
2575
2576
2577
2578
2579
2580
2581
2582
2583! JOZEF ADDITION
2584
2585 IF ( ln_logchlpft ) THEN
2586
2587         WRITE(numout,*)'Pfts loop started'
2588         WRITE(numout,*)'nlogchlpftsets', nlogchlpftsets
2589
2590         DO pft = 1, nn_logchlpftscc         
2591
2592         ! Copy data from logchldatqc to logchldata structures
2593         DO jlogchlset = 1, nlogchlpftsets
2594
2595            CALL obs_surf_decompress( logchlpftdatqc(pft), &
2596                 &                    logchlpftdata(pft), .TRUE., numout )
2597
2598         END DO
2599         
2600         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2601         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2602         DO jlogchlset = 1, nlogchlpftsets
2603            WHERE ( logchlpftdata(pft)%rmod(:,1) == obfillflt )
2604               logchlpftdata(pft)%nqc(:)    = 1
2605               logchlpftdata(pft)%robs(:,1) = obfillflt
2606            END WHERE
2607         END DO
2608
2609         ! Write the logchl data
2610         DO jlogchlset = 1, nlogchlpftsets
2611
2612            WRITE(numout,*)'Code Pfts here!'
2613     
2614           ! WRITE(cdtmp,'(A,I2.2)')'logchlpftfb_',jlogchlset
2615            WRITE(cdtmp,'(A,I2.2,A,I2.2)')'logchlpft',pft,'fb_',jlogchlset
2616            !CALL obs_wri_logchl( cdtmp, logchlpftdata(jlogchlset,pft) )
2617            CALL obs_wri_logchl( cdtmp, logchlpftdata(pft), pft_num=jlogchlset )
2618
2619         END DO
2620
2621         END DO
2622
2623      ENDIF
2624
2625! END
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636
2637
2638
2639
2640
2641
2642
2643
2644
2645      !  - spm
2646      IF ( ln_spm ) THEN
2647
2648         ! Copy data from spmdatqc to spmdata structures
2649         DO jspmset = 1, nspmsets
2650
2651            CALL obs_surf_decompress( spmdatqc(jspmset), &
2652                 &                    spmdata(jspmset), .TRUE., numout )
2653
2654         END DO
2655
2656         ! Write the spm data
2657         DO jspmset = 1, nspmsets
2658     
2659            WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset
2660            CALL obs_wri_spm( cdtmp, spmdata(jspmset) )
2661
2662         END DO
2663
2664      ENDIF
2665
2666      !  - fco2
2667      IF ( ln_fco2 ) THEN
2668
2669         ! Copy data from fco2datqc to fco2data structures
2670         DO jfco2set = 1, nfco2sets
2671
2672            CALL obs_surf_decompress( fco2datqc(jfco2set), &
2673                 &                    fco2data(jfco2set), .TRUE., numout )
2674
2675         END DO
2676         
2677         ! Mark as bad observations with no valid model counterpart due to fCO2 not being in the restart
2678         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2679         DO jfco2set = 1, nfco2sets
2680            WHERE ( fco2data(jfco2set)%rmod(:,1) == obfillflt )
2681               fco2data(jfco2set)%nqc(:)    = 1
2682               fco2data(jfco2set)%robs(:,1) = obfillflt
2683            END WHERE
2684         END DO
2685
2686         ! Write the fco2 data
2687         DO jfco2set = 1, nfco2sets
2688     
2689            WRITE(cdtmp,'(A,I2.2)')'fco2fb_',jfco2set
2690            CALL obs_wri_fco2( cdtmp, fco2data(jfco2set) )
2691
2692         END DO
2693
2694      ENDIF
2695
2696      !  - pco2
2697      IF ( ln_pco2 ) THEN
2698
2699         ! Copy data from pco2datqc to pco2data structures
2700         DO jpco2set = 1, npco2sets
2701
2702            CALL obs_surf_decompress( pco2datqc(jpco2set), &
2703                 &                    pco2data(jpco2set), .TRUE., numout )
2704
2705         END DO
2706         
2707         ! Mark as bad observations with no valid model counterpart due to pco2 not being in the restart
2708         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2709         DO jpco2set = 1, npco2sets
2710            WHERE ( pco2data(jpco2set)%rmod(:,1) == obfillflt )
2711               pco2data(jpco2set)%nqc(:)    = 1
2712               pco2data(jpco2set)%robs(:,1) = obfillflt
2713            END WHERE
2714         END DO
2715
2716         ! Write the pco2 data
2717         DO jpco2set = 1, npco2sets
2718     
2719            WRITE(cdtmp,'(A,I2.2)')'pco2fb_',jpco2set
2720            CALL obs_wri_pco2( cdtmp, pco2data(jpco2set) )
2721
2722         END DO
2723
2724      ENDIF
2725
2726   END SUBROUTINE dia_obs_wri
2727
2728   SUBROUTINE dia_obs_dealloc
2729      IMPLICIT NONE
2730      !!----------------------------------------------------------------------
2731      !!                    *** ROUTINE dia_obs_dealloc ***
2732      !!
2733      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
2734      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
2735      !!
2736      !!  ** Method : Clean up various arrays left behind by the obs_oper.
2737      !!
2738      !!  ** Action :
2739      !!
2740      !!----------------------------------------------------------------------
2741      !! obs_grid deallocation
2742      CALL obs_grid_deallocate
2743
2744      !! diaobs deallocation
2745      IF ( nprofsets > 0 ) THEN
2746          DEALLOCATE(ld_enact, &
2747                  &  profdata, &
2748                  &  prodatqc)
2749      END IF
2750      IF ( ln_sla ) THEN
2751          DEALLOCATE(sladata, &
2752                  &  sladatqc)
2753      END IF
2754      IF ( ln_seaice ) THEN
2755          DEALLOCATE(sladata, &
2756                  &  sladatqc)
2757      END IF
2758      IF ( ln_sst ) THEN
2759          DEALLOCATE(sstdata, &
2760                  &  sstdatqc)
2761      END IF
2762      IF ( ln_vel3d ) THEN
2763          DEALLOCATE(ld_velav, &
2764                  &  velodata, &
2765                  &  veldatqc)
2766      END IF
2767   END SUBROUTINE dia_obs_dealloc
2768
2769   SUBROUTINE ini_date( ddobsini )
2770      !!----------------------------------------------------------------------
2771      !!                    ***  ROUTINE ini_date  ***
2772      !!         
2773      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
2774      !!
2775      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
2776      !!
2777      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
2778      !!
2779      !! History :
2780      !!        !  06-03  (K. Mogensen)  Original code
2781      !!        !  06-05  (K. Mogensen)  Reformatted
2782      !!        !  06-10  (A. Weaver) Cleaning
2783      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
2784      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
2785      !!----------------------------------------------------------------------
2786      USE phycst, ONLY : &            ! Physical constants
2787         & rday
2788!      USE daymod, ONLY : &            ! Time variables
2789!         & nmonth_len           
2790      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
2791         & rdt
2792
2793      IMPLICIT NONE
2794
2795      !! * Arguments
2796      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
2797
2798      !! * Local declarations
2799      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
2800      INTEGER :: imon
2801      INTEGER :: iday
2802      INTEGER :: ihou
2803      INTEGER :: imin
2804      INTEGER :: imday         ! Number of days in month.
2805      REAL(KIND=wp) :: zdayfrc ! Fraction of day
2806
2807      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
2808
2809      !!----------------------------------------------------------------------
2810      !! Initial date initialization (year, month, day, hour, minute)
2811      !! (This assumes that the initial date is for 00z))
2812      !!----------------------------------------------------------------------
2813      iyea =   ndate0 / 10000
2814      imon = ( ndate0 - iyea * 10000 ) / 100
2815      iday =   ndate0 - iyea * 10000 - imon * 100
2816      ihou = 0
2817      imin = 0
2818
2819      !!----------------------------------------------------------------------
2820      !! Compute number of days + number of hours + min since initial time
2821      !!----------------------------------------------------------------------
2822      iday = iday + ( nit000 -1 ) * rdt / rday
2823      zdayfrc = ( nit000 -1 ) * rdt / rday
2824      zdayfrc = zdayfrc - aint(zdayfrc)
2825      ihou = int( zdayfrc * 24 )
2826      imin = int( (zdayfrc * 24 - ihou) * 60 )
2827
2828      !!-----------------------------------------------------------------------
2829      !! Convert number of days (iday) into a real date
2830      !!----------------------------------------------------------------------
2831
2832      CALL calc_month_len( iyea, imonth_len )
2833     
2834      DO WHILE ( iday > imonth_len(imon) )
2835         iday = iday - imonth_len(imon)
2836         imon = imon + 1 
2837         IF ( imon > 12 ) THEN
2838            imon = 1
2839            iyea = iyea + 1
2840            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
2841         ENDIF
2842      END DO
2843
2844      !!----------------------------------------------------------------------
2845      !! Convert it into YYYYMMDD.HHMMSS format.
2846      !!----------------------------------------------------------------------
2847      ddobsini = iyea * 10000_dp + imon * 100_dp + &
2848         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
2849
2850
2851   END SUBROUTINE ini_date
2852
2853   SUBROUTINE fin_date( ddobsfin )
2854      !!----------------------------------------------------------------------
2855      !!                    ***  ROUTINE fin_date  ***
2856      !!         
2857      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
2858      !!
2859      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
2860      !!
2861      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
2862      !!
2863      !! History :
2864      !!        !  06-03  (K. Mogensen)  Original code
2865      !!        !  06-05  (K. Mogensen)  Reformatted
2866      !!        !  06-10  (A. Weaver) Cleaning
2867      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
2868      !!----------------------------------------------------------------------
2869      USE phycst, ONLY : &            ! Physical constants
2870         & rday
2871!      USE daymod, ONLY : &            ! Time variables
2872!         & nmonth_len               
2873      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
2874         & rdt
2875
2876      IMPLICIT NONE
2877
2878      !! * Arguments
2879      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
2880
2881      !! * Local declarations
2882      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
2883      INTEGER :: imon
2884      INTEGER :: iday
2885      INTEGER :: ihou
2886      INTEGER :: imin
2887      INTEGER :: imday         ! Number of days in month.
2888      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
2889         
2890      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
2891           
2892      !-----------------------------------------------------------------------
2893      ! Initial date initialization (year, month, day, hour, minute)
2894      ! (This assumes that the initial date is for 00z)
2895      !-----------------------------------------------------------------------
2896      iyea =   ndate0 / 10000
2897      imon = ( ndate0 - iyea * 10000 ) / 100
2898      iday =   ndate0 - iyea * 10000 - imon * 100
2899      ihou = 0
2900      imin = 0
2901     
2902      !-----------------------------------------------------------------------
2903      ! Compute number of days + number of hours + min since initial time
2904      !-----------------------------------------------------------------------
2905      iday    = iday +  nitend  * rdt / rday
2906      zdayfrc =  nitend  * rdt / rday
2907      zdayfrc = zdayfrc - AINT( zdayfrc )
2908      ihou    = INT( zdayfrc * 24 )
2909      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
2910
2911      !-----------------------------------------------------------------------
2912      ! Convert number of days (iday) into a real date
2913      !----------------------------------------------------------------------
2914
2915      CALL calc_month_len( iyea, imonth_len )
2916     
2917      DO WHILE ( iday > imonth_len(imon) )
2918         iday = iday - imonth_len(imon)
2919         imon = imon + 1 
2920         IF ( imon > 12 ) THEN
2921            imon = 1
2922            iyea = iyea + 1
2923            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
2924         ENDIF
2925      END DO
2926
2927      !-----------------------------------------------------------------------
2928      ! Convert it into YYYYMMDD.HHMMSS format
2929      !-----------------------------------------------------------------------
2930      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
2931         &     + ihou * 0.01_dp  + imin * 0.0001_dp
2932
2933    END SUBROUTINE fin_date
2934   
2935END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.