New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaobs.F90 in branches/UKMO/dev_r4650_general_vert_coord_obsoper_pfts/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

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

Last change on this file was 9054, checked in by dford, 6 years ago

Merge in fix for reading/write logchl STD values.

File size: 99.7 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 1 (STD)
1423         nlogchlextr = 1
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 1 (STD)
1493         nlogchlpftextr = 1
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         & rt0
1797      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
1798         & tsn,  &             
1799         & un, vn,  &
1800         & sshn
1801#if defined  key_lim3
1802      USE ice, ONLY : &                     ! LIM Ice model variables
1803         & frld
1804#endif
1805#if defined key_lim2
1806      USE ice_2, ONLY : &                     ! LIM Ice model variables
1807         & frld
1808#endif
1809#if defined key_hadocc
1810      USE trc, ONLY :  &                ! HadOCC chlorophyll, fCO2 and pCO2
1811         & HADOCC_CHL, &
1812         & HADOCC_FCO2, &
1813         & HADOCC_PCO2, &
1814         & HADOCC_FILL_FLT
1815#elif defined key_medusa && defined key_foam_medusa
1816      USE trc, ONLY :  &                ! MEDUSA chlorophyll, fCO2 and pCO2
1817         & MEDUSA_CHL, &
1818         & MEDUSA_FCO2, &
1819         & MEDUSA_PCO2, &
1820         & MEDUSA_FILL_FLT
1821#elif defined key_fabm
1822      USE fabm
1823      USE par_fabm
1824#endif
1825#if defined key_spm
1826      USE par_spm, ONLY: &              ! ERSEM/SPM sediments
1827         & jp_spm
1828#endif
1829      IMPLICIT NONE
1830
1831      !! * Arguments
1832      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
1833      !! * Local declarations
1834      INTEGER :: idaystp                ! Number of timesteps per day
1835      INTEGER :: jprofset               ! Profile data set loop variable
1836      INTEGER :: jslaset                ! SLA data set loop variable
1837      INTEGER :: jsstset                ! SST data set loop variable
1838      INTEGER :: jseaiceset             ! sea ice data set loop variable
1839      INTEGER :: jveloset               ! velocity profile data loop variable
1840      INTEGER :: jlogchlset             ! logchl data set loop variable
1841      INTEGER :: jspmset                ! spm data set loop variable
1842      INTEGER :: jfco2set               ! fco2 data set loop variable
1843      INTEGER :: jpco2set               ! pco2 data set loop variable
1844      INTEGER :: jvar                   ! Variable number   
1845#if ! defined key_lim2 && ! defined key_lim3
1846      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
1847#endif
1848      REAL(wp) :: tiny                  ! small number
1849      REAL(wp), DIMENSION(jpi,jpj) :: &
1850         logchl                         ! array for log chlorophyll
1851
1852
1853
1854
1855
1856
1857
1858! JOZEF ADDITION
1859
1860
1861      REAL(wp), DIMENSION(jpi,jpj,nn_logchlpftscc) :: &
1862         logchlpfts
1863
1864
1865! END
1866
1867
1868
1869
1870
1871
1872
1873                         ! array for log chlorophyll pfts
1874      REAL(wp), DIMENSION(jpi,jpj) :: &
1875         maskchl                        ! array for special chlorophyll mask
1876      REAL(wp), DIMENSION(jpi,jpj) :: &
1877         spm                            ! array for spm
1878      REAL(wp), DIMENSION(jpi,jpj) :: &
1879         fco2                           ! array for fco2
1880      REAL(wp), DIMENSION(jpi,jpj) :: &
1881         maskfco2                       ! array for special fco2 mask
1882      REAL(wp), DIMENSION(jpi,jpj) :: &
1883         pco2                           ! array for pco2
1884      REAL(wp), DIMENSION(jpi,jpj) :: &
1885         maskpco2                       ! array for special pco2 mask
1886      INTEGER :: jn                     ! loop index
1887#if defined key_fabm
1888      REAL(wp), DIMENSION(jpi,jpj,jpk) :: logchl_3d
1889      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pco2_3d
1890
1891
1892
1893
1894
1895
1896     
1897      ! JOZEF ADDITION
1898
1899      REAL(wp), DIMENSION(jpi,jpj,jpk,nn_logchlpftscc) :: logchlpfts_3d
1900#endif
1901
1902      INTEGER :: pft
1903      ! END
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1914 
1915#if ! defined key_lim2 && ! defined key_lim3
1916      CALL wrk_alloc(jpi,jpj,frld) 
1917#endif
1918
1919      IF(lwp) THEN
1920         WRITE(numout,*)
1921         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1922         WRITE(numout,*) '~~~~~~~'
1923      ENDIF
1924
1925      idaystp = NINT( rday / rdt )
1926
1927      !-----------------------------------------------------------------------
1928      ! No LIM => frld == 0.0_wp
1929      !-----------------------------------------------------------------------
1930#if ! defined key_lim2 && ! defined key_lim3
1931      frld(:,:) = 0.0_wp
1932#endif
1933      !-----------------------------------------------------------------------
1934      ! Depending on switches call various observation operators
1935      !-----------------------------------------------------------------------
1936
1937      !  - Temperature/salinity profiles
1938      IF ( ln_t3d .OR. ln_s3d ) THEN
1939         DO jprofset = 1, nprofsets
1940            IF( (.NOT. lk_vvl) .AND. (ln_zco .OR. ln_zps) ) THEN
1941               IF(lwp) THEN
1942                  WRITE(numout,*) 'dia_obs : calling obs_pro_opt'
1943               ENDIF
1944               IF ( ld_enact(jprofset) ) THEN
1945                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1946                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1947                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1948                     &              gdept_1d, tmask, n1dint, n2dint,        & 
1949                     &              kdailyavtypes = endailyavtypes ) 
1950               ELSE
1951                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1952                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1953                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1954                     &              gdept_1d, tmask, n1dint, n2dint               ) 
1955               ENDIF
1956            ELSE
1957               IF(lwp) THEN
1958                  WRITE(numout,*) 'dia_obs : calling obs_pro_sco_opt'
1959               ENDIF
1960               IF ( ld_enact(jprofset) ) THEN
1961                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1962                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1963                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1964                     &              fsdept(:,:,:), fsdepw(:,:,:),           & 
1965                     &              tmask, n1dint, n2dint,                  & 
1966                     &              kdailyavtypes = endailyavtypes ) 
1967               ELSE
1968                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1969                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1970                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1971                     &              fsdept(:,:,:), fsdepw(:,:,:),           &
1972                     &              tmask, n1dint, n2dint ) 
1973               ENDIF
1974            ENDIF
1975         END DO
1976      ENDIF
1977
1978      !  - Sea surface anomaly
1979      IF ( ln_sla ) THEN
1980         DO jslaset = 1, nslasets
1981            CALL obs_sla_opt( sladatqc(jslaset),            &
1982               &              kstp, jpi, jpj, nit000, sshn, &
1983               &              tmask(:,:,1), n2dint )
1984         END DO         
1985      ENDIF
1986
1987      !  - Sea surface temperature
1988      IF ( ln_sst ) THEN
1989         DO jsstset = 1, nsstsets
1990            CALL obs_sst_opt( sstdatqc(jsstset),                &
1991               &              kstp, jpi, jpj, nit000, idaystp,  &
1992               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
1993               &              n2dint, ld_sstnight(jsstset) )
1994         END DO
1995      ENDIF
1996
1997      !  - Sea surface salinity
1998      IF ( ln_sss ) THEN
1999         IF(lwp) WRITE(numout,*) ' SSS currently not available'
2000      ENDIF
2001
2002#if defined key_lim2 || defined key_lim3
2003      IF ( ln_seaice ) THEN
2004         DO jseaiceset = 1, nseaicesets
2005            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
2006               &              kstp, jpi, jpj, nit000, 1.-frld, &
2007               &              tmask(:,:,1), n2dint )
2008         END DO
2009      ENDIF     
2010#endif
2011
2012      !  - Velocity profiles
2013      IF ( ln_vel3d ) THEN
2014         DO jveloset = 1, nvelosets
2015           ! zonal component of velocity
2016           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
2017              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
2018                             n1dint, n2dint, ld_velav(jveloset) )
2019         END DO
2020      ENDIF
2021
2022      IF ( ln_logchl ) THEN
2023
2024#if defined key_hadocc
2025         logchl(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2026#elif defined key_medusa && defined key_foam_medusa
2027         logchl(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2028#elif defined key_fabm
2029         logchl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
2030         logchl(:,:) = logchl_3d(:,:,1)
2031#else
2032         CALL ctl_stop( ' Trying to run logchl observation operator', &
2033            &           ' but no biogeochemical model appears to have been defined' )
2034#endif
2035
2036         maskchl(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
2037
2038         ! Take the log10 where we can, otherwise exclude
2039         tiny = 1.0e-20
2040         WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
2041            logchl(:,:)  = LOG10(logchl(:,:))
2042         ELSEWHERE
2043            logchl(:,:)  = obfillflt
2044            maskchl(:,:) = 0
2045         END WHERE
2046
2047         DO jlogchlset = 1, nlogchlsets
2048             CALL obs_logchl_opt( logchldatqc(jlogchlset),             &
2049               &                  kstp, jpi, jpj, nit000, logchl(:,:), &
2050               &                  maskchl(:,:), n2dint )
2051         END DO         
2052      ENDIF 
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062! JOZEF ADDITION
2063
2064
2065
2066IF ( ln_logchlpft ) THEN
2067
2068#if defined key_hadocc
2069         logchl(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2070#elif defined key_medusa && defined key_foam_medusa
2071         logchl(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
2072#elif defined key_fabm
2073         logchlpfts_3d(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1)
2074         logchlpfts(:,:,1) = logchlpfts_3d(:,:,1,1)
2075         logchlpfts_3d(:,:,:,2) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl2)
2076         logchlpfts(:,:,2) = logchlpfts_3d(:,:,1,2)
2077         logchlpfts_3d(:,:,:,3) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl3)
2078         logchlpfts(:,:,3) = logchlpfts_3d(:,:,1,3)
2079         logchlpfts_3d(:,:,:,4) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
2080         logchlpfts(:,:,4) = logchlpfts_3d(:,:,1,4)
2081#else
2082         CALL ctl_stop( ' Trying to run logchl observation operator', &
2083            &           ' but no biogeochemical model appears to have been defined' )
2084#endif
2085
2086                  ! create a special mask to exclude certain things
2087
2088         ! Take the log10 where we can, otherwise exclude
2089         tiny = 1.0e-20
2090
2091         DO pft=1,nn_logchlpftscc
2092
2093         maskchl(:,:) = tmask(:,:,1)
2094
2095         WRITE(numout,*) 'Running pft: ',pft
2096
2097         WHERE(logchlpfts(:,:,pft) > tiny .AND. logchlpfts(:,:,pft) /= obfillflt )
2098            logchlpfts(:,:,pft)  = LOG10(logchlpfts(:,:,pft))
2099         ELSEWHERE
2100            logchlpfts(:,:,pft)  = obfillflt
2101            maskchl(:,:) = 0
2102         END WHERE
2103
2104         WRITE(numout,*) 'Max value of Pft: ', MINVAL(logchlpfts(:,:,pft))
2105
2106         DO jlogchlset = 1, nlogchlpftsets
2107             WRITE(numout,*) 'nlogchlpftsets, nn_logchlpftscc, obfillflt = ', nlogchlpftsets, nn_logchlpftscc, obfillflt
2108             WRITE(numout,*) 'kstp, jpi, jpj, nit000, n2dint = ', kstp, jpi, jpj, nit000, n2dint
2109             !WRITE(numout,*) 'ALLOCATED(logchlpftdatqc) = ', ALLOCATED(logchlpftdatqc)
2110             WRITE(numout,*) 'SHAPE(logchlpftdatqc) = ', SHAPE(logchlpftdatqc)
2111             WRITE(numout,*) 'SIZE(logchlpftdatqc) = ', SIZE(logchlpftdatqc)
2112             !WRITE(numout,*) 'ALLOCATED(logchlpfts) = ', ALLOCATED(logchlpfts)
2113             WRITE(numout,*) 'SHAPE(logchlpfts) = ', SHAPE(logchlpfts)
2114             WRITE(numout,*) 'SIZE(logchlpfts) = ', SIZE(logchlpfts)
2115             !WRITE(numout,*) 'ALLOCATED(maskchl) = ', ALLOCATED(maskchl)
2116             WRITE(numout,*) 'SHAPE(maskchl) = ', SHAPE(maskchl)
2117             WRITE(numout,*) 'SIZE(maskchl) = ', SIZE(maskchl)
2118             CALL flush(numout)
2119             
2120             CALL obs_logchl_opt( logchlpftdatqc(pft),             &
2121               &                  kstp, jpi, jpj, nit000, logchlpfts(:,:,pft), &
2122               &                  maskchl(:,:), n2dint )
2123
2124             WRITE(numout,*) 'Max value of Pft: ', MINVAL(logchlpfts(:,:,pft))
2125         END DO         
2126     
2127         
2128         END DO
2129
2130        ENDIF
2131
2132! END
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149      IF ( ln_spm ) THEN
2150#if defined key_spm
2151         spm(:,:) = 0.0
2152         DO jn = 1, jp_spm
2153            spm(:,:) = spm(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
2154         END DO
2155#else
2156         CALL ctl_stop( ' Trying to run spm observation operator', &
2157            &           ' but no spm model appears to have been defined' )
2158#endif
2159
2160         DO jspmset = 1, nspmsets
2161             CALL obs_spm_opt( spmdatqc(jspmset),                &
2162               &               kstp, jpi, jpj, nit000, spm(:,:), &
2163               &               tmask(:,:,1), n2dint )
2164         END DO         
2165      ENDIF
2166
2167      IF ( ln_fco2 ) THEN
2168         maskfco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
2169#if defined key_hadocc
2170         fco2(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
2171         IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
2172            fco2(:,:) = obfillflt
2173            maskfco2(:,:) = 0
2174            CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
2175               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2176               &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
2177         ENDIF
2178#elif defined key_medusa && defined key_foam_medusa
2179         fco2(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA
2180         IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
2181            fco2(:,:) = obfillflt
2182            maskfco2(:,:) = 0
2183            CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
2184               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2185               &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
2186         ENDIF
2187#elif defined key_fabm
2188         ! First, get pCO2 from FABM
2189         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
2190         pco2(:,:) = pco2_3d(:,:,1)
2191         ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
2192         ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
2193         ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
2194         ! and
2195         ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
2196         ! Marine Chemistry, 2: 203-215.
2197         ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
2198         ! not explicitly included - atmospheric pressure is not necessarily available so this is
2199         ! the best assumption.
2200         ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
2201         ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
2202         ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
2203         ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
2204         fco2(:,:) = pco2(:,:) * EXP((-1636.75                                                                               + &
2205            &                         12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
2206            &                         0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
2207            &                         0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
2208            &                         2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
2209            &                        (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
2210#else
2211         CALL ctl_stop( ' Trying to run fco2 observation operator', &
2212            &           ' but no biogeochemical model appears to have been defined' )
2213#endif
2214
2215         DO jfco2set = 1, nfco2sets
2216             CALL obs_fco2_opt( fco2datqc(jfco2set),                      &
2217               &                kstp, jpi, jpj, nit000, fco2(:,:), &
2218               &                maskfco2(:,:), n2dint )
2219         END DO
2220      ENDIF
2221
2222      IF ( ln_pco2 ) THEN
2223         maskpco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
2224#if defined key_hadocc
2225         pco2(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
2226         IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
2227            pco2(:,:) = obfillflt
2228            maskpco2(:,:) = 0
2229            CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
2230               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2231               &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
2232         ENDIF
2233#elif defined key_medusa && defined key_foam_medusa
2234         pco2(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA
2235         IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
2236            pco2(:,:) = obfillflt
2237            maskpco2(:,:) = 0
2238            CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
2239               &           ' on timestep ' // TRIM(STR(kstp)),                              &
2240               &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
2241         ENDIF
2242#elif defined key_fabm
2243         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
2244         pco2(:,:) = pco2_3d(:,:,1)
2245#else
2246         CALL ctl_stop( ' Trying to run pCO2 observation operator', &
2247            &           ' but no biogeochemical model appears to have been defined' )
2248#endif
2249
2250         DO jpco2set = 1, npco2sets
2251             CALL obs_pco2_opt( pco2datqc(jpco2set),                      &
2252               &                kstp, jpi, jpj, nit000, pco2(:,:), &
2253               &                maskpco2(:,:), n2dint )
2254         END DO
2255      ENDIF
2256
2257#if ! defined key_lim2 && ! defined key_lim3
2258      CALL wrk_dealloc(jpi,jpj,frld) 
2259#endif
2260
2261   END SUBROUTINE dia_obs
2262 
2263   SUBROUTINE dia_obs_wri 
2264      !!----------------------------------------------------------------------
2265      !!                    ***  ROUTINE dia_obs_wri  ***
2266      !!         
2267      !! ** Purpose : Call observation diagnostic output routines
2268      !!
2269      !! ** Method  : Call observation diagnostic output routines
2270      !!
2271      !! ** Action  :
2272      !!
2273      !! History :
2274      !!        !  06-03  (K. Mogensen) Original code
2275      !!        !  06-05  (K. Mogensen) Reformatted
2276      !!        !  06-10  (A. Weaver) Cleaning
2277      !!        !  07-03  (K. Mogensen) General handling of profiles
2278      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
2279      !!----------------------------------------------------------------------
2280      IMPLICIT NONE
2281
2282      !! * Local declarations
2283
2284      INTEGER :: jprofset                 ! Profile data set loop variable
2285      INTEGER :: jveloset                 ! Velocity data set loop variable
2286      INTEGER :: jslaset                  ! SLA data set loop variable
2287      INTEGER :: jsstset                  ! SST data set loop variable
2288      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
2289      INTEGER :: jlogchlset               ! logchl data set loop variable
2290      INTEGER :: jspmset                  ! spm data set loop variable
2291      INTEGER :: jfco2set                 ! fco2 data set loop variable
2292      INTEGER :: jpco2set                 ! pco2 data set loop variable
2293      INTEGER :: jset
2294      INTEGER :: jfbini
2295      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
2296      CHARACTER(LEN=20) :: cdtmp
2297      INTEGER :: pft
2298      !-----------------------------------------------------------------------
2299      ! Depending on switches call various observation output routines
2300      !-----------------------------------------------------------------------
2301
2302      !  - Temperature/salinity profiles
2303
2304      IF( ln_t3d .OR. ln_s3d ) THEN
2305
2306         ! Copy data from prodatqc to profdata structures
2307         DO jprofset = 1, nprofsets
2308
2309            CALL obs_prof_decompress( prodatqc(jprofset), &
2310                 &                    profdata(jprofset), .TRUE., numout )
2311
2312         END DO
2313
2314         ! Write the profiles.
2315
2316         jprofset = 0
2317
2318         ! ENACT insitu data
2319
2320         IF ( ln_ena ) THEN
2321           
2322            jprofset = jprofset + 1
2323
2324            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
2325
2326         ENDIF
2327
2328         ! Coriolis insitu data
2329
2330         IF ( ln_cor ) THEN
2331           
2332            jprofset = jprofset + 1
2333
2334            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
2335           
2336         ENDIF
2337         
2338         ! Feedback insitu data
2339
2340         IF ( ln_profb ) THEN
2341
2342            jfbini = jprofset + 1
2343
2344            DO jprofset = jfbini, nprofsets
2345               
2346               jset = jprofset - jfbini + 1
2347               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
2348               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
2349
2350            END DO
2351
2352         ENDIF
2353
2354      ENDIF
2355
2356      !  - Sea surface anomaly
2357      IF ( ln_sla ) THEN
2358
2359         ! Copy data from sladatqc to sladata structures
2360         DO jslaset = 1, nslasets
2361
2362              CALL obs_surf_decompress( sladatqc(jslaset), &
2363                 &                    sladata(jslaset), .TRUE., numout )
2364
2365         END DO
2366
2367         jslaset = 0 
2368
2369         ! Write the AVISO SLA data
2370
2371         IF ( ln_sladt ) THEN
2372           
2373            jslaset = 1
2374            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
2375            jslaset = 2
2376            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
2377
2378         ENDIF
2379
2380         IF ( ln_slafb ) THEN
2381           
2382            jfbini = jslaset + 1
2383
2384            DO jslaset = jfbini, nslasets
2385               
2386               jset = jslaset - jfbini + 1
2387               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
2388               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
2389
2390            END DO
2391
2392         ENDIF
2393
2394      ENDIF
2395
2396      !  - Sea surface temperature
2397      IF ( ln_sst ) THEN
2398
2399         ! Copy data from sstdatqc to sstdata structures
2400         DO jsstset = 1, nsstsets
2401     
2402              CALL obs_surf_decompress( sstdatqc(jsstset), &
2403                 &                    sstdata(jsstset), .TRUE., numout )
2404
2405         END DO
2406
2407         jsstset = 0 
2408
2409         ! Write the AVISO SST data
2410
2411         IF ( ln_reysst ) THEN
2412           
2413            jsstset = jsstset + 1
2414            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
2415
2416         ENDIF
2417
2418         IF ( ln_ghrsst ) THEN
2419           
2420            jsstset = jsstset + 1
2421            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
2422
2423         ENDIF
2424
2425         IF ( ln_sstfb ) THEN
2426           
2427            jfbini = jsstset + 1
2428
2429            DO jsstset = jfbini, nsstsets
2430               
2431               jset = jsstset - jfbini + 1
2432               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
2433               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
2434
2435            END DO
2436
2437         ENDIF
2438
2439      ENDIF
2440
2441      !  - Sea surface salinity
2442      IF ( ln_sss ) THEN
2443         IF(lwp) WRITE(numout,*) ' SSS currently not available'
2444      ENDIF
2445
2446      !  - Sea Ice Concentration
2447      IF ( ln_seaice ) THEN
2448
2449         ! Copy data from seaicedatqc to seaicedata structures
2450         DO jseaiceset = 1, nseaicesets
2451
2452              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
2453                 &                    seaicedata(jseaiceset), .TRUE., numout )
2454
2455         END DO
2456
2457         ! Write the Sea Ice data
2458         DO jseaiceset = 1, nseaicesets
2459     
2460            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
2461            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
2462
2463         END DO
2464
2465      ENDIF
2466     
2467      ! Velocity data
2468      IF( ln_vel3d ) THEN
2469
2470         ! Copy data from veldatqc to velodata structures
2471         DO jveloset = 1, nvelosets
2472
2473            CALL obs_prof_decompress( veldatqc(jveloset), &
2474                 &                    velodata(jveloset), .TRUE., numout )
2475
2476         END DO
2477
2478         ! Write the profiles.
2479
2480         jveloset = 0
2481
2482         ! Daily averaged data
2483
2484         IF ( ln_velavcur ) THEN
2485           
2486            jveloset = jveloset + 1
2487
2488            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
2489
2490         ENDIF
2491
2492         ! High frequency data
2493
2494         IF ( ln_velhrcur ) THEN
2495           
2496            jveloset = jveloset + 1
2497
2498            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
2499
2500         ENDIF
2501
2502         ! Daily averaged data
2503
2504         IF ( ln_velavadcp ) THEN
2505           
2506            jveloset = jveloset + 1
2507
2508            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
2509
2510         ENDIF
2511
2512         ! High frequency data
2513
2514         IF ( ln_velhradcp ) THEN
2515           
2516            jveloset = jveloset + 1
2517           
2518            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
2519               
2520         ENDIF
2521
2522         ! Feedback velocity data
2523
2524         IF ( ln_velfb ) THEN
2525
2526            jfbini = jveloset + 1
2527
2528            DO jveloset = jfbini, nvelosets
2529               
2530               jset = jveloset - jfbini + 1
2531               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
2532               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
2533
2534            END DO
2535
2536         ENDIF
2537         
2538      ENDIF
2539
2540      !  - log10(chlorophyll)
2541      IF ( ln_logchl ) THEN
2542
2543         ! Copy data from logchldatqc to logchldata structures
2544         DO jlogchlset = 1, nlogchlsets
2545
2546            CALL obs_surf_decompress( logchldatqc(jlogchlset), &
2547                 &                    logchldata(jlogchlset), .TRUE., numout )
2548
2549         END DO
2550         
2551         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2552         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2553         DO jlogchlset = 1, nlogchlsets
2554            WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
2555               logchldata(jlogchlset)%nqc(:)    = 1
2556               logchldata(jlogchlset)%robs(:,1) = obfillflt
2557            END WHERE
2558         END DO
2559
2560         ! Write the logchl data
2561         ! If padd is ever defined and added as an optional argument
2562         ! to the obs_wri_logchl call, then the observed STD values
2563         ! will need adding to padd in order to be written out
2564         DO jlogchlset = 1, nlogchlsets
2565     
2566            WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
2567            CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
2568
2569         END DO
2570
2571      ENDIF
2572
2573
2574
2575
2576
2577
2578
2579     
2580
2581
2582
2583
2584
2585
2586
2587
2588! JOZEF ADDITION
2589
2590 IF ( ln_logchlpft ) THEN
2591
2592         WRITE(numout,*)'Pfts loop started'
2593         WRITE(numout,*)'nlogchlpftsets', nlogchlpftsets
2594
2595         DO pft = 1, nn_logchlpftscc         
2596
2597         ! Copy data from logchldatqc to logchldata structures
2598         DO jlogchlset = 1, nlogchlpftsets
2599
2600            CALL obs_surf_decompress( logchlpftdatqc(pft), &
2601                 &                    logchlpftdata(pft), .TRUE., numout )
2602
2603         END DO
2604         
2605         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2606         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2607         DO jlogchlset = 1, nlogchlpftsets
2608            WHERE ( logchlpftdata(pft)%rmod(:,1) == obfillflt )
2609               logchlpftdata(pft)%nqc(:)    = 1
2610               logchlpftdata(pft)%robs(:,1) = obfillflt
2611            END WHERE
2612         END DO
2613
2614         ! Write the logchl data
2615         DO jlogchlset = 1, nlogchlpftsets
2616
2617            WRITE(numout,*)'Code Pfts here!'
2618     
2619           ! WRITE(cdtmp,'(A,I2.2)')'logchlpftfb_',jlogchlset
2620            WRITE(cdtmp,'(A,I2.2,A,I2.2)')'logchlpft',pft,'fb_',jlogchlset
2621            !CALL obs_wri_logchl( cdtmp, logchlpftdata(jlogchlset,pft) )
2622            CALL obs_wri_logchl( cdtmp, logchlpftdata(pft), pft_num=jlogchlset )
2623
2624         END DO
2625
2626         END DO
2627
2628      ENDIF
2629
2630! END
2631
2632
2633
2634
2635
2636
2637
2638
2639
2640
2641
2642
2643
2644
2645
2646
2647
2648
2649
2650      !  - spm
2651      IF ( ln_spm ) THEN
2652
2653         ! Copy data from spmdatqc to spmdata structures
2654         DO jspmset = 1, nspmsets
2655
2656            CALL obs_surf_decompress( spmdatqc(jspmset), &
2657                 &                    spmdata(jspmset), .TRUE., numout )
2658
2659         END DO
2660
2661         ! Write the spm data
2662         DO jspmset = 1, nspmsets
2663     
2664            WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset
2665            CALL obs_wri_spm( cdtmp, spmdata(jspmset) )
2666
2667         END DO
2668
2669      ENDIF
2670
2671      !  - fco2
2672      IF ( ln_fco2 ) THEN
2673
2674         ! Copy data from fco2datqc to fco2data structures
2675         DO jfco2set = 1, nfco2sets
2676
2677            CALL obs_surf_decompress( fco2datqc(jfco2set), &
2678                 &                    fco2data(jfco2set), .TRUE., numout )
2679
2680         END DO
2681         
2682         ! Mark as bad observations with no valid model counterpart due to fCO2 not being in the restart
2683         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2684         DO jfco2set = 1, nfco2sets
2685            WHERE ( fco2data(jfco2set)%rmod(:,1) == obfillflt )
2686               fco2data(jfco2set)%nqc(:)    = 1
2687               fco2data(jfco2set)%robs(:,1) = obfillflt
2688            END WHERE
2689         END DO
2690
2691         ! Write the fco2 data
2692         DO jfco2set = 1, nfco2sets
2693     
2694            WRITE(cdtmp,'(A,I2.2)')'fco2fb_',jfco2set
2695            CALL obs_wri_fco2( cdtmp, fco2data(jfco2set) )
2696
2697         END DO
2698
2699      ENDIF
2700
2701      !  - pco2
2702      IF ( ln_pco2 ) THEN
2703
2704         ! Copy data from pco2datqc to pco2data structures
2705         DO jpco2set = 1, npco2sets
2706
2707            CALL obs_surf_decompress( pco2datqc(jpco2set), &
2708                 &                    pco2data(jpco2set), .TRUE., numout )
2709
2710         END DO
2711         
2712         ! Mark as bad observations with no valid model counterpart due to pco2 not being in the restart
2713         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2714         DO jpco2set = 1, npco2sets
2715            WHERE ( pco2data(jpco2set)%rmod(:,1) == obfillflt )
2716               pco2data(jpco2set)%nqc(:)    = 1
2717               pco2data(jpco2set)%robs(:,1) = obfillflt
2718            END WHERE
2719         END DO
2720
2721         ! Write the pco2 data
2722         DO jpco2set = 1, npco2sets
2723     
2724            WRITE(cdtmp,'(A,I2.2)')'pco2fb_',jpco2set
2725            CALL obs_wri_pco2( cdtmp, pco2data(jpco2set) )
2726
2727         END DO
2728
2729      ENDIF
2730
2731   END SUBROUTINE dia_obs_wri
2732
2733   SUBROUTINE dia_obs_dealloc
2734      IMPLICIT NONE
2735      !!----------------------------------------------------------------------
2736      !!                    *** ROUTINE dia_obs_dealloc ***
2737      !!
2738      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
2739      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
2740      !!
2741      !!  ** Method : Clean up various arrays left behind by the obs_oper.
2742      !!
2743      !!  ** Action :
2744      !!
2745      !!----------------------------------------------------------------------
2746      !! obs_grid deallocation
2747      CALL obs_grid_deallocate
2748
2749      !! diaobs deallocation
2750      IF ( nprofsets > 0 ) THEN
2751          DEALLOCATE(ld_enact, &
2752                  &  profdata, &
2753                  &  prodatqc)
2754      END IF
2755      IF ( ln_sla ) THEN
2756          DEALLOCATE(sladata, &
2757                  &  sladatqc)
2758      END IF
2759      IF ( ln_seaice ) THEN
2760          DEALLOCATE(sladata, &
2761                  &  sladatqc)
2762      END IF
2763      IF ( ln_sst ) THEN
2764          DEALLOCATE(sstdata, &
2765                  &  sstdatqc)
2766      END IF
2767      IF ( ln_vel3d ) THEN
2768          DEALLOCATE(ld_velav, &
2769                  &  velodata, &
2770                  &  veldatqc)
2771      END IF
2772   END SUBROUTINE dia_obs_dealloc
2773
2774   SUBROUTINE ini_date( ddobsini )
2775      !!----------------------------------------------------------------------
2776      !!                    ***  ROUTINE ini_date  ***
2777      !!         
2778      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
2779      !!
2780      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
2781      !!
2782      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
2783      !!
2784      !! History :
2785      !!        !  06-03  (K. Mogensen)  Original code
2786      !!        !  06-05  (K. Mogensen)  Reformatted
2787      !!        !  06-10  (A. Weaver) Cleaning
2788      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
2789      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
2790      !!----------------------------------------------------------------------
2791      USE phycst, ONLY : &            ! Physical constants
2792         & rday
2793!      USE daymod, ONLY : &            ! Time variables
2794!         & nmonth_len           
2795      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
2796         & rdt
2797
2798      IMPLICIT NONE
2799
2800      !! * Arguments
2801      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
2802
2803      !! * Local declarations
2804      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
2805      INTEGER :: imon
2806      INTEGER :: iday
2807      INTEGER :: ihou
2808      INTEGER :: imin
2809      INTEGER :: imday         ! Number of days in month.
2810      REAL(KIND=wp) :: zdayfrc ! Fraction of day
2811
2812      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
2813
2814      !!----------------------------------------------------------------------
2815      !! Initial date initialization (year, month, day, hour, minute)
2816      !! (This assumes that the initial date is for 00z))
2817      !!----------------------------------------------------------------------
2818      iyea =   ndate0 / 10000
2819      imon = ( ndate0 - iyea * 10000 ) / 100
2820      iday =   ndate0 - iyea * 10000 - imon * 100
2821      ihou = 0
2822      imin = 0
2823
2824      !!----------------------------------------------------------------------
2825      !! Compute number of days + number of hours + min since initial time
2826      !!----------------------------------------------------------------------
2827      iday = iday + ( nit000 -1 ) * rdt / rday
2828      zdayfrc = ( nit000 -1 ) * rdt / rday
2829      zdayfrc = zdayfrc - aint(zdayfrc)
2830      ihou = int( zdayfrc * 24 )
2831      imin = int( (zdayfrc * 24 - ihou) * 60 )
2832
2833      !!-----------------------------------------------------------------------
2834      !! Convert number of days (iday) into a real date
2835      !!----------------------------------------------------------------------
2836
2837      CALL calc_month_len( iyea, imonth_len )
2838     
2839      DO WHILE ( iday > imonth_len(imon) )
2840         iday = iday - imonth_len(imon)
2841         imon = imon + 1 
2842         IF ( imon > 12 ) THEN
2843            imon = 1
2844            iyea = iyea + 1
2845            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
2846         ENDIF
2847      END DO
2848
2849      !!----------------------------------------------------------------------
2850      !! Convert it into YYYYMMDD.HHMMSS format.
2851      !!----------------------------------------------------------------------
2852      ddobsini = iyea * 10000_dp + imon * 100_dp + &
2853         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
2854
2855
2856   END SUBROUTINE ini_date
2857
2858   SUBROUTINE fin_date( ddobsfin )
2859      !!----------------------------------------------------------------------
2860      !!                    ***  ROUTINE fin_date  ***
2861      !!         
2862      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
2863      !!
2864      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
2865      !!
2866      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
2867      !!
2868      !! History :
2869      !!        !  06-03  (K. Mogensen)  Original code
2870      !!        !  06-05  (K. Mogensen)  Reformatted
2871      !!        !  06-10  (A. Weaver) Cleaning
2872      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
2873      !!----------------------------------------------------------------------
2874      USE phycst, ONLY : &            ! Physical constants
2875         & rday
2876!      USE daymod, ONLY : &            ! Time variables
2877!         & nmonth_len               
2878      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
2879         & rdt
2880
2881      IMPLICIT NONE
2882
2883      !! * Arguments
2884      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
2885
2886      !! * Local declarations
2887      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
2888      INTEGER :: imon
2889      INTEGER :: iday
2890      INTEGER :: ihou
2891      INTEGER :: imin
2892      INTEGER :: imday         ! Number of days in month.
2893      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
2894         
2895      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
2896           
2897      !-----------------------------------------------------------------------
2898      ! Initial date initialization (year, month, day, hour, minute)
2899      ! (This assumes that the initial date is for 00z)
2900      !-----------------------------------------------------------------------
2901      iyea =   ndate0 / 10000
2902      imon = ( ndate0 - iyea * 10000 ) / 100
2903      iday =   ndate0 - iyea * 10000 - imon * 100
2904      ihou = 0
2905      imin = 0
2906     
2907      !-----------------------------------------------------------------------
2908      ! Compute number of days + number of hours + min since initial time
2909      !-----------------------------------------------------------------------
2910      iday    = iday +  nitend  * rdt / rday
2911      zdayfrc =  nitend  * rdt / rday
2912      zdayfrc = zdayfrc - AINT( zdayfrc )
2913      ihou    = INT( zdayfrc * 24 )
2914      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
2915
2916      !-----------------------------------------------------------------------
2917      ! Convert number of days (iday) into a real date
2918      !----------------------------------------------------------------------
2919
2920      CALL calc_month_len( iyea, imonth_len )
2921     
2922      DO WHILE ( iday > imonth_len(imon) )
2923         iday = iday - imonth_len(imon)
2924         imon = imon + 1 
2925         IF ( imon > 12 ) THEN
2926            imon = 1
2927            iyea = iyea + 1
2928            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
2929         ENDIF
2930      END DO
2931
2932      !-----------------------------------------------------------------------
2933      ! Convert it into YYYYMMDD.HHMMSS format
2934      !-----------------------------------------------------------------------
2935      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
2936         &     + ihou * 0.01_dp  + imin * 0.0001_dp
2937
2938    END SUBROUTINE fin_date
2939   
2940END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.