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

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9013

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

Update to read in and write out logchl obs STD values. See internal Met Office NEMO ticket 736.

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