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

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6854

Last change on this file since 6854 was 6854, checked in by dford, 8 years ago

Initial implementation of observation operator for LogChl?.

File size: 69.8 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_prep                 ! Preparation of obs. (grid search etc).
34   USE obs_oper                 ! Observation operators
35   USE obs_write                ! Writing of observation related diagnostics
36   USE obs_grid                 ! Grid searching
37   USE obs_read_altbias         ! Bias treatment for altimeter
38   USE obs_profiles_def         ! Profile data definitions
39   USE obs_profiles             ! Profile data storage
40   USE obs_surf_def             ! Surface data definitions
41   USE obs_sla                  ! SLA data storage
42   USE obs_sst                  ! SST data storage
43   USE obs_seaice               ! Sea Ice data storage
44   USE obs_logchl               ! logchl data storage
45   USE obs_types                ! Definitions for observation types
46   USE mpp_map                  ! MPP mapping
47   USE lib_mpp                  ! For ctl_warn/stop
48
49   IMPLICIT NONE
50
51   !! * Routine accessibility
52   PRIVATE
53   PUBLIC dia_obs_init, &  ! Initialize and read observations
54      &   dia_obs,      &  ! Compute model equivalent to observations
55      &   dia_obs_wri,  &  ! Write model equivalent to observations
56      &   dia_obs_dealloc  ! Deallocate dia_obs data
57
58   !! * Shared Module variables
59   LOGICAL, PUBLIC, PARAMETER :: &
60#if defined key_diaobs
61      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
62#else
63      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
64#endif
65
66   !! * Module variables
67   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
68   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
69   LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set
70   LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set
71   LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles
72   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies
73   LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files
74   LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files
75   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
76   LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature
77   LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data
78   LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files
79   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
80   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
81   LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data
82   LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data
83   LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data
84   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data
85   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
86   LOGICAL, PUBLIC :: ln_logchl      !: Logical switch for log10(chlorophyll)
87   LOGICAL, PUBLIC :: ln_logchlfb    !: Logical switch for logchl from feedback files
88   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
89   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
90   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
91   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
92   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias 
93   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
94   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
95   LOGICAL, PUBLIC :: ln_sstbias     !: Logical switch for bias corection of SST
96
97   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
98   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
99 
100   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
101   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method
102
103   INTEGER, DIMENSION(imaxavtypes) :: &
104      & endailyavtypes !: ENACT data types which are daily average
105
106   INTEGER, PARAMETER :: MaxNumFiles = 1000
107   LOGICAL, DIMENSION(MaxNumFiles) :: &
108      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
109   !                    !: If so use endailyavtypes
110      & ln_profb_enatim !: Change tim for 820 enact data set.
111   
112   INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type
113
114   LOGICAL, DIMENSION(MaxNumFiles) :: &
115      & ln_velfb_av   !: Is the velocity feedback files daily average?
116   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
117      & ld_enact     !: Profile data is ENACT so use endailyavtypes
118   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
119      & ld_velav     !: Velocity data is daily averaged
120   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
121      & ld_sstnight  !: SST observation corresponds to night mean
122
123   !!----------------------------------------------------------------------
124   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
125   !! $Id$
126   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
127   !!----------------------------------------------------------------------
128
129   !! * Substitutions
130#  include "domzgr_substitute.h90"
131CONTAINS
132
133   SUBROUTINE dia_obs_init
134      !!----------------------------------------------------------------------
135      !!                    ***  ROUTINE dia_obs_init  ***
136      !!         
137      !! ** Purpose : Initialize and read observations
138      !!
139      !! ** Method  : Read the namelist and call reading routines
140      !!
141      !! ** Action  : Read the namelist and call reading routines
142      !!
143      !! History :
144      !!        !  06-03  (K. Mogensen) Original code
145      !!        !  06-05  (A. Weaver) Reformatted
146      !!        !  06-10  (A. Weaver) Cleaning and add controls
147      !!        !  07-03  (K. Mogensen) General handling of profiles
148      !!        !  14-08  (J.While) Incorporated SST bias correction
149      !!----------------------------------------------------------------------
150
151      IMPLICIT NONE
152
153      !! * Local declarations
154      CHARACTER(len=128) :: enactfiles(MaxNumFiles)
155      CHARACTER(len=128) :: coriofiles(MaxNumFiles)
156      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
157      CHARACTER(len=128) :: sstfiles(MaxNumFiles)     
158      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
159      CHARACTER(len=128) :: sstbias_files(MaxNumFiles) 
160      CHARACTER(len=128) :: slafilesact(MaxNumFiles)     
161      CHARACTER(len=128) :: slafilespas(MaxNumFiles)     
162      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
163      CHARACTER(len=128) :: seaicefiles(MaxNumFiles)           
164      CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 
165      CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)   
166      CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
167      CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
168      CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
169      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
170      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
171      CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
172      CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
173      CHARACTER(LEN=128) :: reysstname
174      CHARACTER(LEN=12)  :: reysstfmt
175      CHARACTER(LEN=128) :: bias_file
176      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
177      NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       &
178         &            ln_sla, ln_sladt, ln_slafb,                     &
179         &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       &
180         &            enactfiles, coriofiles, profbfiles,             &
181         &            slafilesact, slafilespas, slafbfiles,           &
182         &            sstfiles, sstfbfiles,                           &
183         &            ln_seaice, seaicefiles,                         &
184         &            dobsini, dobsend, n1dint, n2dint,               &
185         &            nmsshc, mdtcorr, mdtcutoff,                     &
186         &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    &
187         &            ln_sstnight,                                    &
188         &            ln_grid_search_lookup,                          &
189         &            grid_search_file, grid_search_res,              &
190         &            ln_grid_global, bias_file, ln_altbias,          &
191         &            endailyavtypes, ln_s_at_t, ln_profb_ena,        &
192         &            ln_vel3d, ln_velavcur, velavcurfiles,           &
193         &            ln_velhrcur, velhrcurfiles,                     &
194         &            ln_velavadcp, velavadcpfiles,                   &
195         &            ln_velhradcp, velhradcpfiles,                   &
196         &            ln_velfb, velfbfiles, ln_velfb_av,              &
197         &            ln_logchl, ln_logchlfb,                         &
198         &            logchlfiles, logchlfbfiles,                     &
199         &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
200         &            ln_sstbias, sstbias_files
201
202      INTEGER :: jprofset
203      INTEGER :: jveloset
204      INTEGER :: jvar
205      INTEGER :: jnumenact
206      INTEGER :: jnumcorio
207      INTEGER :: jnumprofb
208      INTEGER :: jnumslaact
209      INTEGER :: jnumslapas
210      INTEGER :: jnumslafb
211      INTEGER :: jnumsst
212      INTEGER :: jnumsstfb
213      INTEGER :: jnumsstbias
214      INTEGER :: jnumseaice
215      INTEGER :: jnumvelavcur
216      INTEGER :: jnumvelhrcur 
217      INTEGER :: jnumvelavadcp
218      INTEGER :: jnumvelhradcp   
219      INTEGER :: jnumvelfb
220      INTEGER :: jnumlogchl
221      INTEGER :: jnumlogchlfb
222      INTEGER :: ji
223      INTEGER :: jset
224      INTEGER :: ios                 ! Local integer output status for namelist read
225      LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
226
227      !-----------------------------------------------------------------------
228      ! Read namelist parameters
229      !-----------------------------------------------------------------------
230
231      ln_logchl   = .FALSE.
232      ln_logchlfb = .FALSE.
233     
234      !Initalise all values in namelist arrays
235      enactfiles(:) = ''
236      coriofiles(:) = ''
237      profbfiles(:) = ''
238      slafilesact(:) = ''
239      slafilespas(:) = ''
240      slafbfiles(:) = ''
241      sstfiles(:)   = ''
242      sstfbfiles(:) = ''
243      seaicefiles(:) = ''
244      velcurfiles(:) = ''
245      veladcpfiles(:) = ''
246      velavcurfiles(:) = ''
247      velhrcurfiles(:) = ''
248      velavadcpfiles(:) = ''
249      velhradcpfiles(:) = ''
250      velfbfiles(:) = ''
251      velcurfiles(:) = ''
252      veladcpfiles(:) = ''
253      logchlfiles(:) = ''
254      logchlfbfiles(:) = ''
255      sstbias_files(:) = ''
256      endailyavtypes(:) = -1
257      endailyavtypes(1) = 820
258      ln_profb_ena(:) = .FALSE.
259      ln_profb_enatim(:) = .TRUE.
260      ln_velfb_av(:) = .FALSE.
261      ln_ignmis = .FALSE.
262
263      ! Read Namelist namobs : control observation diagnostics
264      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
265      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
266901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
267
268      CALL ini_date( dobsini )
269      CALL fin_date( dobsend )
270 
271      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
272      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
273902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
274      IF(lwm) WRITE ( numond, namobs )
275
276      ! Count number of files for each type
277      IF (ln_ena) THEN
278         lmask(:) = .FALSE.
279         WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
280         jnumenact = COUNT(lmask)
281      ENDIF
282      IF (ln_cor) THEN
283         lmask(:) = .FALSE.
284         WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
285         jnumcorio = COUNT(lmask)
286      ENDIF
287      IF (ln_profb) THEN
288         lmask(:) = .FALSE.
289         WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
290         jnumprofb = COUNT(lmask)
291      ENDIF
292      IF (ln_sladt) THEN
293         lmask(:) = .FALSE.
294         WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
295         jnumslaact = COUNT(lmask)
296         lmask(:) = .FALSE.
297         WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
298         jnumslapas = COUNT(lmask)
299      ENDIF
300      IF (ln_slafb) THEN
301         lmask(:) = .FALSE.
302         WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
303         jnumslafb = COUNT(lmask)
304         lmask(:) = .FALSE.
305      ENDIF
306      IF (ln_ghrsst) THEN
307         lmask(:) = .FALSE.
308         WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
309         jnumsst = COUNT(lmask)
310      ENDIF     
311      IF (ln_sstfb) THEN
312         lmask(:) = .FALSE.
313         WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
314         jnumsstfb = COUNT(lmask)
315         lmask(:) = .FALSE.
316      ENDIF
317      IF (ln_sstbias) THEN 
318         lmask(:) = .FALSE. 
319         WHERE (sstbias_files(:) /= '') lmask(:) = .TRUE. 
320         jnumsstbias = COUNT(lmask) 
321         lmask(:) = .FALSE. 
322      ENDIF       
323      IF (ln_seaice) THEN
324         lmask(:) = .FALSE.
325         WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
326         jnumseaice = COUNT(lmask)
327      ENDIF
328      IF (ln_velavcur) THEN
329         lmask(:) = .FALSE.
330         WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
331         jnumvelavcur = COUNT(lmask)
332      ENDIF
333      IF (ln_velhrcur) THEN
334         lmask(:) = .FALSE.
335         WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
336         jnumvelhrcur = COUNT(lmask)
337      ENDIF
338      IF (ln_velavadcp) THEN
339         lmask(:) = .FALSE.
340         WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
341         jnumvelavadcp = COUNT(lmask)
342      ENDIF
343      IF (ln_velhradcp) THEN
344         lmask(:) = .FALSE.
345         WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
346         jnumvelhradcp = COUNT(lmask)
347      ENDIF
348      IF (ln_velfb) THEN
349         lmask(:) = .FALSE.
350         WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
351         jnumvelfb = COUNT(lmask)
352         lmask(:) = .FALSE.
353      ENDIF
354      IF (ln_logchl) THEN
355         lmask(:) = .FALSE.
356         WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE.
357         jnumlogchl = COUNT(lmask)
358      ENDIF
359      IF (ln_logchlfb) THEN
360         lmask(:) = .FALSE.
361         WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE.
362         jnumlogchlfb = COUNT(lmask)
363      ENDIF
364     
365      ! Control print
366      IF(lwp) THEN
367         WRITE(numout,*)
368         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
369         WRITE(numout,*) '~~~~~~~~~~~~'
370         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
371         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
372         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
373         WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena
374         WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor
375         WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb
376         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
377         WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt
378         WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb
379         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
380         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
381         WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst   
382         WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst
383         WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb
384         WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias
385         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
386         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
387         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
388         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
389         WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur
390         WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur
391         WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp
392         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
393         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
394         WRITE(numout,*) '             Logical switch for logchl observations          ln_logchl = ', ln_logchl
395         WRITE(numout,*) '             Logical switch for feedback logchl data       ln_logchlfb = ', ln_logchlfb
396         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
397         WRITE(numout,*) &
398   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
399         IF (ln_grid_search_lookup) &
400            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file
401         IF (ln_ena) THEN
402            DO ji = 1, jnumenact
403               WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', &
404                  TRIM(enactfiles(ji))
405            END DO
406         ENDIF
407         IF (ln_cor) THEN
408            DO ji = 1, jnumcorio
409               WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', &
410                  TRIM(coriofiles(ji))
411            END DO
412         ENDIF
413         IF (ln_profb) THEN
414            DO ji = 1, jnumprofb
415               IF (ln_profb_ena(ji)) THEN
416                  WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', &
417                     TRIM(profbfiles(ji))
418               ELSE
419                  WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', &
420                     TRIM(profbfiles(ji))
421               ENDIF
422               WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
423            END DO
424         ENDIF
425         IF (ln_sladt) THEN
426            DO ji = 1, jnumslaact
427               WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', &
428                  TRIM(slafilesact(ji))
429            END DO
430            DO ji = 1, jnumslapas
431               WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', &
432                  TRIM(slafilespas(ji))
433            END DO
434         ENDIF
435         IF (ln_slafb) THEN
436            DO ji = 1, jnumslafb
437               WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', &
438                  TRIM(slafbfiles(ji))
439            END DO
440         ENDIF
441         IF (ln_ghrsst) THEN
442            DO ji = 1, jnumsst
443               WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', &
444                  TRIM(sstfiles(ji))
445            END DO
446         ENDIF
447         IF (ln_sstfb) THEN
448            DO ji = 1, jnumsstfb
449               WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', &
450                  TRIM(sstfbfiles(ji))
451            END DO
452         ENDIF
453         IF (ln_seaice) THEN
454            DO ji = 1, jnumseaice
455               WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', &
456                  TRIM(seaicefiles(ji))
457            END DO
458         ENDIF
459         IF (ln_velavcur) THEN
460            DO ji = 1, jnumvelavcur
461               WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', &
462                  TRIM(velavcurfiles(ji))
463            END DO
464         ENDIF
465         IF (ln_velhrcur) THEN
466            DO ji = 1, jnumvelhrcur
467               WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', &
468                  TRIM(velhrcurfiles(ji))
469            END DO
470         ENDIF
471         IF (ln_velavadcp) THEN
472            DO ji = 1, jnumvelavadcp
473               WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', &
474                  TRIM(velavadcpfiles(ji))
475            END DO
476         ENDIF
477         IF (ln_velhradcp) THEN
478            DO ji = 1, jnumvelhradcp
479               WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', &
480                  TRIM(velhradcpfiles(ji))
481            END DO
482         ENDIF
483         IF (ln_velfb) THEN
484            DO ji = 1, jnumvelfb
485               IF (ln_velfb_av(ji)) THEN
486                  WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', &
487                     TRIM(velfbfiles(ji))
488               ELSE
489                  WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', &
490                     TRIM(velfbfiles(ji))
491               ENDIF
492            END DO
493         ENDIF
494         IF (ln_logchl) THEN
495            DO ji = 1, jnumlogchl
496               WRITE(numout,'(1X,2A)') '             logchl input observation file name        logchlfiles = ', &
497                  TRIM(logchlfiles(ji))
498            END DO
499         ENDIF
500         IF (ln_logchlfb) THEN
501            DO ji = 1, jnumlogchlfb
502               WRITE(numout,'(1X,2A)') '        Feedback logchl input observation file name  logchlfbfiles = ', &
503                  TRIM(logchlfbfiles(ji))
504            END DO
505         ENDIF
506         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini
507         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
508         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
509         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
510         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
511         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
512         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
513         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
514         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
515         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
516         WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes
517
518      ENDIF
519     
520      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
521         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
522         RETURN
523      ENDIF
524
525      IF ( ln_grid_global ) THEN
526         CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
527      ENDIF
528
529      CALL obs_typ_init
530     
531      IF ( ln_grid_global ) THEN     
532         CALL mppmap_init
533      ENDIF
534     
535      ! Parameter control
536#if defined key_diaobs
537      IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
538         & ( .NOT. ln_vel3d ).AND.                                         &
539         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
540         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ) ) THEN
541         IF(lwp) WRITE(numout,cform_war)
542         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
543            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
544            &                    ' ln_logchl are all set to .FALSE.'
545         nwarn = nwarn + 1
546      ENDIF
547#endif
548
549      CALL obs_grid_setup( )
550      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
551         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
552            &                    ' is not available')
553      ENDIF
554      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
555         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
556            &                    ' is not available')
557      ENDIF
558
559      !-----------------------------------------------------------------------
560      ! Depending on switches read the various observation types
561      !-----------------------------------------------------------------------
562      !  - Temperature/salinity profiles
563
564      IF ( ln_t3d .OR. ln_s3d ) THEN
565
566         ! Set the number of variables for profiles to 2 (T and S)
567         nprofvars = 2
568         ! Set the number of extra variables for profiles to 1 (insitu temp).
569         nprofextr = 1
570
571         ! Count how may insitu data sets we have and allocate data.
572         jprofset = 0
573         IF ( ln_ena ) jprofset = jprofset + 1
574         IF ( ln_cor ) jprofset = jprofset + 1
575         IF ( ln_profb ) jprofset = jprofset + jnumprofb
576         nprofsets = jprofset
577         IF ( nprofsets > 0 ) THEN
578            ALLOCATE(ld_enact(nprofsets))
579            ALLOCATE(profdata(nprofsets))
580            ALLOCATE(prodatqc(nprofsets))
581         ENDIF
582
583         jprofset = 0
584         
585         ! ENACT insitu data
586
587         IF ( ln_ena ) THEN
588
589            jprofset = jprofset + 1
590           
591            ld_enact(jprofset) = .TRUE.
592
593            CALL obs_rea_pro_dri( 1, profdata(jprofset),          &
594               &                  jnumenact, enactfiles(1:jnumenact), &
595               &                  nprofvars, nprofextr,        &
596               &                  nitend-nit000+2,             &
597               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
598               &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
599               &                  kdailyavtypes = endailyavtypes )
600
601            DO jvar = 1, 2
602
603               CALL obs_prof_staend( profdata(jprofset), jvar )
604
605            END DO
606
607            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
608               &              ln_t3d, ln_s3d, ln_nea, &
609               &              kdailyavtypes=endailyavtypes )
610           
611         ENDIF
612
613         ! Coriolis insitu data
614
615         IF ( ln_cor ) THEN
616           
617            jprofset = jprofset + 1
618
619            ld_enact(jprofset) = .FALSE.
620
621            CALL obs_rea_pro_dri( 2, profdata(jprofset),          &
622               &                  jnumcorio, coriofiles(1:jnumcorio), &
623               &                  nprofvars, nprofextr,        &
624               &                  nitend-nit000+2,             &
625               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
626               &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
627
628            DO jvar = 1, 2
629
630               CALL obs_prof_staend( profdata(jprofset), jvar )
631
632            END DO
633
634            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
635                 &            ln_t3d, ln_s3d, ln_nea )
636           
637         ENDIF
638 
639         ! Feedback insitu data
640
641         IF ( ln_profb ) THEN
642           
643            DO jset = 1, jnumprofb
644               
645               jprofset = jprofset + 1
646               ld_enact (jprofset) = ln_profb_ena(jset)
647
648               CALL obs_rea_pro_dri( 0, profdata(jprofset),          &
649                  &                  1, profbfiles(jset:jset), &
650                  &                  nprofvars, nprofextr,        &
651                  &                  nitend-nit000+2,             &
652                  &                  dobsini, dobsend, ln_t3d, ln_s3d, &
653                  &                  ln_ignmis, ln_s_at_t, &
654                  &                  ld_enact(jprofset).AND.&
655                  &                  ln_profb_enatim(jset), &
656                  &                  .FALSE., kdailyavtypes = endailyavtypes )
657               
658               DO jvar = 1, 2
659                 
660                  CALL obs_prof_staend( profdata(jprofset), jvar )
661                 
662               END DO
663               
664               IF ( ld_enact(jprofset) ) THEN
665                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
666                     &              ln_t3d, ln_s3d, ln_nea, &
667                     &              kdailyavtypes = endailyavtypes )
668               ELSE
669                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
670                     &              ln_t3d, ln_s3d, ln_nea )
671               ENDIF
672               
673            END DO
674
675         ENDIF
676
677      ENDIF
678
679      !  - Sea level anomalies
680      IF ( ln_sla ) THEN
681        ! Set the number of variables for sla to 1
682         nslavars = 1
683
684         ! Set the number of extra variables for sla to 2
685         nslaextr = 2
686         
687         ! Set the number of sla data sets to 2
688         nslasets = 0
689         IF ( ln_sladt ) THEN
690            nslasets = nslasets + 2
691         ENDIF
692         IF ( ln_slafb ) THEN
693            nslasets = nslasets + jnumslafb
694         ENDIF
695         
696         ALLOCATE(sladata(nslasets))
697         ALLOCATE(sladatqc(nslasets))
698         sladata(:)%nsurf=0
699         sladatqc(:)%nsurf=0
700
701         nslasets = 0
702
703         ! AVISO SLA data
704
705         IF ( ln_sladt ) THEN
706
707            ! Active SLA observations
708           
709            nslasets = nslasets + 1
710           
711            CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
712               &              slafilesact(1:jnumslaact), &
713               &              nslavars, nslaextr, nitend-nit000+2, &
714               &              dobsini, dobsend, ln_ignmis, .FALSE. )
715            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
716               &              ln_sla, ln_nea )
717           
718            ! Passive SLA observations
719           
720            nslasets = nslasets + 1
721           
722            CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
723               &              slafilespas(1:jnumslapas), &
724               &              nslavars, nslaextr, nitend-nit000+2, &
725               &              dobsini, dobsend, ln_ignmis, .FALSE. )
726           
727            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
728               &              ln_sla, ln_nea )
729
730         ENDIF
731         
732         ! Feedback SLA data
733
734         IF ( ln_slafb ) THEN
735
736            DO jset = 1, jnumslafb
737           
738               nslasets = nslasets + 1
739           
740               CALL obs_rea_sla( 0, sladata(nslasets), 1, &
741                  &              slafbfiles(jset:jset), &
742                  &              nslavars, nslaextr, nitend-nit000+2, &
743                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
744               CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
745                  &              ln_sla, ln_nea )
746
747            END DO               
748
749         ENDIF
750         
751         CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
752           
753         ! read in altimeter bias
754         
755         IF ( ln_altbias ) THEN     
756            CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
757         ENDIF
758     
759      ENDIF
760
761      !  - Sea surface height
762      IF ( ln_ssh ) THEN
763         IF(lwp) WRITE(numout,*) ' SSH currently not available'
764      ENDIF
765
766      !  - Sea surface temperature
767      IF ( ln_sst ) THEN
768
769         ! Set the number of variables for sst to 1
770         nsstvars = 1
771
772         ! Set the number of extra variables for sst to 0
773         nsstextr = 0
774
775         nsstsets = 0
776
777         IF (ln_reysst) nsstsets = nsstsets + 1
778         IF (ln_ghrsst) nsstsets = nsstsets + 1
779         IF ( ln_sstfb ) THEN
780            nsstsets = nsstsets + jnumsstfb
781         ENDIF
782
783         ALLOCATE(sstdata(nsstsets))
784         ALLOCATE(sstdatqc(nsstsets))
785         ALLOCATE(ld_sstnight(nsstsets))
786         sstdata(:)%nsurf=0
787         sstdatqc(:)%nsurf=0   
788         ld_sstnight(:)=.false.
789
790         nsstsets = 0
791
792         IF (ln_reysst) THEN
793
794            nsstsets = nsstsets + 1
795
796            ld_sstnight(nsstsets) = ln_sstnight
797
798            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
799               &                  nsstvars, nsstextr, &
800               &                  nitend-nit000+2, dobsini, dobsend )
801            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
802               &              ln_nea )
803
804        ENDIF
805       
806        IF (ln_ghrsst) THEN
807       
808            nsstsets = nsstsets + 1
809
810            ld_sstnight(nsstsets) = ln_sstnight
811         
812            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
813               &              sstfiles(1:jnumsst), &
814               &              nsstvars, nsstextr, nitend-nit000+2, &
815               &              dobsini, dobsend, ln_ignmis, .FALSE. )
816            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
817               &              ln_nea )
818
819        ENDIF
820               
821         ! Feedback SST data
822
823         IF ( ln_sstfb ) THEN
824
825            DO jset = 1, jnumsstfb
826           
827               nsstsets = nsstsets + 1
828
829               ld_sstnight(nsstsets) = ln_sstnight
830           
831               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
832                  &              sstfbfiles(jset:jset), &
833                  &              nsstvars, nsstextr, nitend-nit000+2, &
834                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
835               CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
836                  &              ln_sst, ln_nea )
837
838            END DO               
839
840         ENDIF
841         
842         !Read in bias field and correct SST.
843         IF ( ln_sstbias ) THEN
844            IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
845                                             "  but no bias"// & 
846                                             " files to read in")   
847            CALL obs_app_sstbias( nsstsets, sstdatqc, n2dint, & 
848                                  jnumsstbias, & 
849                                  sstbias_files(1:jnumsstbias) ) 
850         ENDIF
851
852      ENDIF
853
854      !  - Sea surface salinity
855      IF ( ln_sss ) THEN
856         IF(lwp) WRITE(numout,*) ' SSS currently not available'
857      ENDIF
858
859      !  - Sea Ice Concentration
860     
861      IF ( ln_seaice ) THEN
862
863         ! Set the number of variables for seaice to 1
864         nseaicevars = 1
865
866         ! Set the number of extra variables for seaice to 0
867         nseaiceextr = 0
868         
869         ! Set the number of data sets to 1
870         nseaicesets = 1
871
872         ALLOCATE(seaicedata(nseaicesets))
873         ALLOCATE(seaicedatqc(nseaicesets))
874         seaicedata(:)%nsurf=0
875         seaicedatqc(:)%nsurf=0
876
877         CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
878            &                 seaicefiles(1:jnumseaice), &
879            &                 nseaicevars, nseaiceextr, nitend-nit000+2, &
880            &                 dobsini, dobsend, ln_ignmis, .FALSE. )
881
882         CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
883            &                 ln_seaice, ln_nea )
884 
885      ENDIF
886
887      IF (ln_vel3d) THEN
888
889         ! Set the number of variables for profiles to 2 (U and V)
890         nvelovars = 2
891
892         ! Set the number of extra variables for profiles to 2 to store
893         ! rotation parameters
894         nveloextr = 2
895
896         jveloset = 0
897         
898         IF ( ln_velavcur ) jveloset = jveloset + 1
899         IF ( ln_velhrcur ) jveloset = jveloset + 1
900         IF ( ln_velavadcp ) jveloset = jveloset + 1
901         IF ( ln_velhradcp ) jveloset = jveloset + 1
902         IF (ln_velfb) jveloset = jveloset + jnumvelfb
903
904         nvelosets = jveloset
905         IF ( nvelosets > 0 ) THEN
906            ALLOCATE( velodata(nvelosets) )
907            ALLOCATE( veldatqc(nvelosets) )
908            ALLOCATE( ld_velav(nvelosets) )
909         ENDIF
910         
911         jveloset = 0
912         
913         ! Daily averaged data
914
915         IF ( ln_velavcur ) THEN
916           
917            jveloset = jveloset + 1
918           
919            ld_velav(jveloset) = .TRUE.
920           
921            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
922               &                  velavcurfiles(1:jnumvelavcur), &
923               &                  nvelovars, nveloextr, &
924               &                  nitend-nit000+2,              &
925               &                  dobsini, dobsend, ln_ignmis, &
926               &                  ld_velav(jveloset), &
927               &                  .FALSE. )
928           
929            DO jvar = 1, 2
930               CALL obs_prof_staend( velodata(jveloset), jvar )
931            END DO
932           
933            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
934               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
935           
936         ENDIF
937
938         ! High frequency data
939
940         IF ( ln_velhrcur ) THEN
941           
942            jveloset = jveloset + 1
943           
944            ld_velav(jveloset) = .FALSE.
945               
946            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
947               &                  velhrcurfiles(1:jnumvelhrcur), &
948               &                  nvelovars, nveloextr, &
949               &                  nitend-nit000+2,              &
950               &                  dobsini, dobsend, ln_ignmis, &
951               &                  ld_velav(jveloset), &
952               &                  .FALSE. )
953           
954            DO jvar = 1, 2
955               CALL obs_prof_staend( velodata(jveloset), jvar )
956            END DO
957           
958            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
959               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
960           
961         ENDIF
962
963         ! Daily averaged data
964
965         IF ( ln_velavadcp ) THEN
966           
967            jveloset = jveloset + 1
968           
969            ld_velav(jveloset) = .TRUE.
970           
971            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
972               &                  velavadcpfiles(1:jnumvelavadcp), &
973               &                  nvelovars, nveloextr, &
974               &                  nitend-nit000+2,              &
975               &                  dobsini, dobsend, ln_ignmis, &
976               &                  ld_velav(jveloset), &
977               &                  .FALSE. )
978           
979            DO jvar = 1, 2
980               CALL obs_prof_staend( velodata(jveloset), jvar )
981            END DO
982           
983            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
984               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
985           
986         ENDIF
987
988         ! High frequency data
989
990         IF ( ln_velhradcp ) THEN
991           
992            jveloset = jveloset + 1
993           
994            ld_velav(jveloset) = .FALSE.
995               
996            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
997               &                  velhradcpfiles(1:jnumvelhradcp), &
998               &                  nvelovars, nveloextr, &
999               &                  nitend-nit000+2,              &
1000               &                  dobsini, dobsend, ln_ignmis, &
1001               &                  ld_velav(jveloset), &
1002               &                  .FALSE. )
1003           
1004            DO jvar = 1, 2
1005               CALL obs_prof_staend( velodata(jveloset), jvar )
1006            END DO
1007           
1008            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1009               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1010           
1011         ENDIF
1012
1013         IF ( ln_velfb ) THEN
1014
1015            DO jset = 1, jnumvelfb
1016           
1017               jveloset = jveloset + 1
1018
1019               ld_velav(jveloset) = ln_velfb_av(jset)
1020               
1021               CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
1022                  &                  velfbfiles(jset:jset), &
1023                  &                  nvelovars, nveloextr, &
1024                  &                  nitend-nit000+2,              &
1025                  &                  dobsini, dobsend, ln_ignmis, &
1026                  &                  ld_velav(jveloset), &
1027                  &                  .FALSE. )
1028               
1029               DO jvar = 1, 2
1030                  CALL obs_prof_staend( velodata(jveloset), jvar )
1031               END DO
1032               
1033               CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1034                  &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1035
1036
1037            END DO
1038           
1039         ENDIF
1040
1041      ENDIF
1042
1043      !  - log10(chlorophyll)
1044     
1045      IF ( ln_logchl ) THEN
1046
1047         ! Set the number of variables for logchl to 1
1048         nlogchlvars = 1
1049
1050         ! Set the number of extra variables for logchl to 0
1051         nlogchlextr = 0
1052         
1053         IF ( ln_logchlfb ) THEN
1054            nlogchlsets = jnumlogchlfb
1055         ELSE
1056            nlogchlsets = 1
1057         ENDIF
1058
1059         ALLOCATE(logchldata(nlogchlsets))
1060         ALLOCATE(logchldatqc(nlogchlsets))
1061         logchldata(:)%nsurf=0
1062         logchldatqc(:)%nsurf=0
1063
1064         nlogchlsets = 0
1065
1066         IF ( ln_logchlfb ) THEN             ! Feedback file format
1067
1068            DO jset = 1, jnumlogchlfb
1069           
1070               nlogchlsets = nlogchlsets + 1
1071
1072               CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, &
1073                  &                 logchlfbfiles(jset:jset), &
1074                  &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
1075                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1076
1077               CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1078                  &                 ln_logchl, ln_nea )
1079           
1080            ENDDO
1081
1082         ELSE                              ! Original file format
1083
1084            nlogchlsets = nlogchlsets + 1
1085
1086            CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, &
1087               &                 logchlfiles(1:jnumlogchl), &
1088               &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
1089               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1090
1091            CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1092               &                 ln_logchl, ln_nea )
1093
1094         ENDIF
1095 
1096      ENDIF
1097     
1098   END SUBROUTINE dia_obs_init
1099
1100   SUBROUTINE dia_obs( kstp )
1101      !!----------------------------------------------------------------------
1102      !!                    ***  ROUTINE dia_obs  ***
1103      !!         
1104      !! ** Purpose : Call the observation operators on each time step
1105      !!
1106      !! ** Method  : Call the observation operators on each time step to
1107      !!              compute the model equivalent of the following date:
1108      !!               - T profiles
1109      !!               - S profiles
1110      !!               - Sea surface height (referenced to a mean)
1111      !!               - Sea surface temperature
1112      !!               - Sea surface salinity
1113      !!               - Velocity component (U,V) profiles
1114      !!               - Sea surface log10(chlorophyll)
1115      !!
1116      !! ** Action  :
1117      !!
1118      !! History :
1119      !!        !  06-03  (K. Mogensen) Original code
1120      !!        !  06-05  (K. Mogensen) Reformatted
1121      !!        !  06-10  (A. Weaver) Cleaning
1122      !!        !  07-03  (K. Mogensen) General handling of profiles
1123      !!        !  07-04  (G. Smith) Generalized surface operators
1124      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
1125      !!        !  14-08  (J. While) observation operator for profiles in
1126      !!                             generalised vertical coordinates
1127      !!----------------------------------------------------------------------
1128      !! * Modules used
1129      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
1130         & rdt,           &                       
1131         & gdept_1d,       &             
1132#if defined key_vvl 
1133         & gdept_n,       &
1134#else 
1135         & gdept_1d,      &
1136#endif                                       
1137         & tmask, umask, vmask                           
1138      USE phycst, ONLY : &              ! Physical constants
1139         & rday                         
1140      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
1141         & tsn,  &             
1142         & un, vn,  &
1143         & sshn
1144#if defined  key_lim3
1145      USE ice, ONLY : &                     ! LIM Ice model variables
1146         & frld
1147#endif
1148#if defined key_lim2
1149      USE ice_2, ONLY : &                     ! LIM Ice model variables
1150         & frld
1151#endif
1152#if defined key_hadocc
1153      USE trc, ONLY :  &                ! HadOCC chlorophyll
1154         & HADOCC_CHL, &
1155         & HADOCC_FILL_FLT
1156#elif defined key_medusa && defined key_foam_medusa
1157      USE trc, ONLY :  &                ! MEDUSA chlorophyll
1158         & MEDUSA_CHL, &
1159         & MEDUSA_FILL_FLT
1160#elif defined key_fabm
1161      !USE ???                           ! ERSEM chlorophyll
1162#endif
1163      IMPLICIT NONE
1164
1165      !! * Arguments
1166      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
1167      !! * Local declarations
1168      INTEGER :: idaystp                ! Number of timesteps per day
1169      INTEGER :: jprofset               ! Profile data set loop variable
1170      INTEGER :: jslaset                ! SLA data set loop variable
1171      INTEGER :: jsstset                ! SST data set loop variable
1172      INTEGER :: jseaiceset             ! sea ice data set loop variable
1173      INTEGER :: jveloset               ! velocity profile data loop variable
1174      INTEGER :: jlogchlset             ! logchl data set loop variable
1175      INTEGER :: jvar                   ! Variable number   
1176#if ! defined key_lim2 && ! defined key_lim3
1177      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
1178#endif
1179      REAL(wp) :: tiny                  ! small number
1180      REAL(wp), DIMENSION(jpi,jpj) :: &
1181         logchl                         ! array for log chlorophyll
1182      REAL(wp), DIMENSION(jpi,jpj) :: &
1183         maskchl                        ! array for special chlorophyll mask
1184      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1185 
1186#if ! defined key_lim2 && ! defined key_lim3
1187      CALL wrk_alloc(jpi,jpj,frld) 
1188#endif
1189
1190      IF(lwp) THEN
1191         WRITE(numout,*)
1192         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1193         WRITE(numout,*) '~~~~~~~'
1194      ENDIF
1195
1196      idaystp = NINT( rday / rdt )
1197
1198      !-----------------------------------------------------------------------
1199      ! No LIM => frld == 0.0_wp
1200      !-----------------------------------------------------------------------
1201#if ! defined key_lim2 && ! defined key_lim3
1202      frld(:,:) = 0.0_wp
1203#endif
1204      !-----------------------------------------------------------------------
1205      ! Depending on switches call various observation operators
1206      !-----------------------------------------------------------------------
1207
1208      !  - Temperature/salinity profiles
1209      IF ( ln_t3d .OR. ln_s3d ) THEN
1210         DO jprofset = 1, nprofsets
1211            IF( (.NOT. lk_vvl) .AND. (ln_zco .OR. ln_zps) ) THEN
1212               IF(lwp) THEN
1213                  WRITE(numout,*) 'dia_obs : calling obs_pro_opt'
1214               ENDIF
1215               IF ( ld_enact(jprofset) ) THEN
1216                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1217                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1218                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1219                     &              gdept_1d, tmask, n1dint, n2dint,        & 
1220                     &              kdailyavtypes = endailyavtypes ) 
1221               ELSE
1222                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1223                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1224                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1225                     &              gdept_1d, tmask, n1dint, n2dint               ) 
1226               ENDIF
1227            ELSE
1228               IF(lwp) THEN
1229                  WRITE(numout,*) 'dia_obs : calling obs_pro_sco_opt'
1230               ENDIF
1231               IF ( ld_enact(jprofset) ) THEN
1232                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1233                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1234                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1235                     &              fsdept(:,:,:), fsdepw(:,:,:),           & 
1236                     &              tmask, n1dint, n2dint,                  & 
1237                     &              kdailyavtypes = endailyavtypes ) 
1238               ELSE
1239                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1240                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1241                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1242                     &              fsdept(:,:,:), fsdepw(:,:,:),           &
1243                     &              tmask, n1dint, n2dint ) 
1244               ENDIF
1245            ENDIF
1246         END DO
1247      ENDIF
1248
1249      !  - Sea surface anomaly
1250      IF ( ln_sla ) THEN
1251         DO jslaset = 1, nslasets
1252            CALL obs_sla_opt( sladatqc(jslaset),            &
1253               &              kstp, jpi, jpj, nit000, sshn, &
1254               &              tmask(:,:,1), n2dint )
1255         END DO         
1256      ENDIF
1257
1258      !  - Sea surface temperature
1259      IF ( ln_sst ) THEN
1260         DO jsstset = 1, nsstsets
1261            CALL obs_sst_opt( sstdatqc(jsstset),                &
1262               &              kstp, jpi, jpj, nit000, idaystp,  &
1263               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
1264               &              n2dint, ld_sstnight(jsstset) )
1265         END DO
1266      ENDIF
1267
1268      !  - Sea surface salinity
1269      IF ( ln_sss ) THEN
1270         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1271      ENDIF
1272
1273#if defined key_lim2 || defined key_lim3
1274      IF ( ln_seaice ) THEN
1275         DO jseaiceset = 1, nseaicesets
1276            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
1277               &              kstp, jpi, jpj, nit000, 1.-frld, &
1278               &              tmask(:,:,1), n2dint )
1279         END DO
1280      ENDIF     
1281#endif
1282
1283      !  - Velocity profiles
1284      IF ( ln_vel3d ) THEN
1285         DO jveloset = 1, nvelosets
1286           ! zonal component of velocity
1287           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
1288              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
1289                             n1dint, n2dint, ld_velav(jveloset) )
1290         END DO
1291      ENDIF
1292
1293      IF ( ln_logchl ) THEN
1294
1295#if defined key_hadocc
1296         logchl(:,:)  = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
1297#elif defined key_medusa && defined key_foam_medusa
1298         logchl(:,:)  = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
1299#elif defined key_fabm
1300         !logchl(:,:)  =  ???                 ! (not log) chlorophyll from ERSEM
1301         CALL ctl_stop( ' Trying to run logchl observation operator', &
1302            &           ' but not properly implemented for FABM-ERSEM yet' )
1303#else
1304         CALL ctl_stop( ' Trying to run logchl observation operator', &
1305            &           ' but no biogeochemical model appears to have been defined' )
1306#endif
1307
1308         maskchl(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
1309
1310         ! Take the log10 where we can, otherwise exclude
1311         tiny = 1.0e-20
1312         WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
1313            logchl(:,:)  = LOG10(logchl(:,:))
1314         ELSEWHERE
1315            logchl(:,:)  = obfillflt
1316            maskchl(:,:) = 0
1317         END WHERE
1318
1319         DO jlogchlset = 1, nlogchlsets
1320             CALL obs_logchl_opt( logchldatqc(jlogchlset),             &
1321               &                  kstp, jpi, jpj, nit000, logchl(:,:), &
1322               &                  maskchl(:,:), n2dint )
1323         END DO         
1324      ENDIF 
1325
1326#if ! defined key_lim2 && ! defined key_lim3
1327      CALL wrk_dealloc(jpi,jpj,frld) 
1328#endif
1329
1330   END SUBROUTINE dia_obs
1331 
1332   SUBROUTINE dia_obs_wri 
1333      !!----------------------------------------------------------------------
1334      !!                    ***  ROUTINE dia_obs_wri  ***
1335      !!         
1336      !! ** Purpose : Call observation diagnostic output routines
1337      !!
1338      !! ** Method  : Call observation diagnostic output routines
1339      !!
1340      !! ** Action  :
1341      !!
1342      !! History :
1343      !!        !  06-03  (K. Mogensen) Original code
1344      !!        !  06-05  (K. Mogensen) Reformatted
1345      !!        !  06-10  (A. Weaver) Cleaning
1346      !!        !  07-03  (K. Mogensen) General handling of profiles
1347      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
1348      !!----------------------------------------------------------------------
1349      IMPLICIT NONE
1350
1351      !! * Local declarations
1352
1353      INTEGER :: jprofset                 ! Profile data set loop variable
1354      INTEGER :: jveloset                 ! Velocity data set loop variable
1355      INTEGER :: jslaset                  ! SLA data set loop variable
1356      INTEGER :: jsstset                  ! SST data set loop variable
1357      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
1358      INTEGER :: jlogchlset               ! logchl data set loop variable
1359      INTEGER :: jset
1360      INTEGER :: jfbini
1361      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1362      CHARACTER(LEN=10) :: cdtmp
1363      !-----------------------------------------------------------------------
1364      ! Depending on switches call various observation output routines
1365      !-----------------------------------------------------------------------
1366
1367      !  - Temperature/salinity profiles
1368
1369      IF( ln_t3d .OR. ln_s3d ) THEN
1370
1371         ! Copy data from prodatqc to profdata structures
1372         DO jprofset = 1, nprofsets
1373
1374            CALL obs_prof_decompress( prodatqc(jprofset), &
1375                 &                    profdata(jprofset), .TRUE., numout )
1376
1377         END DO
1378
1379         ! Write the profiles.
1380
1381         jprofset = 0
1382
1383         ! ENACT insitu data
1384
1385         IF ( ln_ena ) THEN
1386           
1387            jprofset = jprofset + 1
1388
1389            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
1390
1391         ENDIF
1392
1393         ! Coriolis insitu data
1394
1395         IF ( ln_cor ) THEN
1396           
1397            jprofset = jprofset + 1
1398
1399            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
1400           
1401         ENDIF
1402         
1403         ! Feedback insitu data
1404
1405         IF ( ln_profb ) THEN
1406
1407            jfbini = jprofset + 1
1408
1409            DO jprofset = jfbini, nprofsets
1410               
1411               jset = jprofset - jfbini + 1
1412               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
1413               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
1414
1415            END DO
1416
1417         ENDIF
1418
1419      ENDIF
1420
1421      !  - Sea surface anomaly
1422      IF ( ln_sla ) THEN
1423
1424         ! Copy data from sladatqc to sladata structures
1425         DO jslaset = 1, nslasets
1426
1427              CALL obs_surf_decompress( sladatqc(jslaset), &
1428                 &                    sladata(jslaset), .TRUE., numout )
1429
1430         END DO
1431
1432         jslaset = 0 
1433
1434         ! Write the AVISO SLA data
1435
1436         IF ( ln_sladt ) THEN
1437           
1438            jslaset = 1
1439            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
1440            jslaset = 2
1441            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
1442
1443         ENDIF
1444
1445         IF ( ln_slafb ) THEN
1446           
1447            jfbini = jslaset + 1
1448
1449            DO jslaset = jfbini, nslasets
1450               
1451               jset = jslaset - jfbini + 1
1452               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
1453               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
1454
1455            END DO
1456
1457         ENDIF
1458
1459      ENDIF
1460
1461      !  - Sea surface temperature
1462      IF ( ln_sst ) THEN
1463
1464         ! Copy data from sstdatqc to sstdata structures
1465         DO jsstset = 1, nsstsets
1466     
1467              CALL obs_surf_decompress( sstdatqc(jsstset), &
1468                 &                    sstdata(jsstset), .TRUE., numout )
1469
1470         END DO
1471
1472         jsstset = 0 
1473
1474         ! Write the AVISO SST data
1475
1476         IF ( ln_reysst ) THEN
1477           
1478            jsstset = jsstset + 1
1479            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
1480
1481         ENDIF
1482
1483         IF ( ln_ghrsst ) THEN
1484           
1485            jsstset = jsstset + 1
1486            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
1487
1488         ENDIF
1489
1490         IF ( ln_sstfb ) THEN
1491           
1492            jfbini = jsstset + 1
1493
1494            DO jsstset = jfbini, nsstsets
1495               
1496               jset = jsstset - jfbini + 1
1497               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
1498               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
1499
1500            END DO
1501
1502         ENDIF
1503
1504      ENDIF
1505
1506      !  - Sea surface salinity
1507      IF ( ln_sss ) THEN
1508         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1509      ENDIF
1510
1511      !  - Sea Ice Concentration
1512      IF ( ln_seaice ) THEN
1513
1514         ! Copy data from seaicedatqc to seaicedata structures
1515         DO jseaiceset = 1, nseaicesets
1516
1517              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
1518                 &                    seaicedata(jseaiceset), .TRUE., numout )
1519
1520         END DO
1521
1522         ! Write the Sea Ice data
1523         DO jseaiceset = 1, nseaicesets
1524     
1525            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
1526            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
1527
1528         END DO
1529
1530      ENDIF
1531     
1532      ! Velocity data
1533      IF( ln_vel3d ) THEN
1534
1535         ! Copy data from veldatqc to velodata structures
1536         DO jveloset = 1, nvelosets
1537
1538            CALL obs_prof_decompress( veldatqc(jveloset), &
1539                 &                    velodata(jveloset), .TRUE., numout )
1540
1541         END DO
1542
1543         ! Write the profiles.
1544
1545         jveloset = 0
1546
1547         ! Daily averaged data
1548
1549         IF ( ln_velavcur ) THEN
1550           
1551            jveloset = jveloset + 1
1552
1553            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
1554
1555         ENDIF
1556
1557         ! High frequency data
1558
1559         IF ( ln_velhrcur ) THEN
1560           
1561            jveloset = jveloset + 1
1562
1563            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
1564
1565         ENDIF
1566
1567         ! Daily averaged data
1568
1569         IF ( ln_velavadcp ) THEN
1570           
1571            jveloset = jveloset + 1
1572
1573            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
1574
1575         ENDIF
1576
1577         ! High frequency data
1578
1579         IF ( ln_velhradcp ) THEN
1580           
1581            jveloset = jveloset + 1
1582           
1583            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
1584               
1585         ENDIF
1586
1587         ! Feedback velocity data
1588
1589         IF ( ln_velfb ) THEN
1590
1591            jfbini = jveloset + 1
1592
1593            DO jveloset = jfbini, nvelosets
1594               
1595               jset = jveloset - jfbini + 1
1596               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
1597               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
1598
1599            END DO
1600
1601         ENDIF
1602         
1603      ENDIF
1604
1605      !  - log10(chlorophyll)
1606      IF ( ln_logchl ) THEN
1607
1608         ! Copy data from logchldatqc to logchldata structures
1609         DO jlogchlset = 1, nlogchlsets
1610
1611            CALL obs_surf_decompress( logchldatqc(jlogchlset), &
1612                 &                    logchldata(jlogchlset), .TRUE., numout )
1613
1614         END DO
1615         
1616         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
1617         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
1618         DO jlogchlset = 1, nlogchlsets
1619            WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
1620               logchldata(jlogchlset)%nqc(:)    = 1
1621               logchldata(jlogchlset)%robs(:,1) = obfillflt
1622            END WHERE
1623         END DO
1624
1625         ! Write the logchl data
1626         DO jlogchlset = 1, nlogchlsets
1627     
1628            WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
1629            CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
1630
1631         END DO
1632
1633      ENDIF
1634
1635   END SUBROUTINE dia_obs_wri
1636
1637   SUBROUTINE dia_obs_dealloc
1638      IMPLICIT NONE
1639      !!----------------------------------------------------------------------
1640      !!                    *** ROUTINE dia_obs_dealloc ***
1641      !!
1642      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1643      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1644      !!
1645      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1646      !!
1647      !!  ** Action :
1648      !!
1649      !!----------------------------------------------------------------------
1650      !! obs_grid deallocation
1651      CALL obs_grid_deallocate
1652
1653      !! diaobs deallocation
1654      IF ( nprofsets > 0 ) THEN
1655          DEALLOCATE(ld_enact, &
1656                  &  profdata, &
1657                  &  prodatqc)
1658      END IF
1659      IF ( ln_sla ) THEN
1660          DEALLOCATE(sladata, &
1661                  &  sladatqc)
1662      END IF
1663      IF ( ln_seaice ) THEN
1664          DEALLOCATE(sladata, &
1665                  &  sladatqc)
1666      END IF
1667      IF ( ln_sst ) THEN
1668          DEALLOCATE(sstdata, &
1669                  &  sstdatqc)
1670      END IF
1671      IF ( ln_vel3d ) THEN
1672          DEALLOCATE(ld_velav, &
1673                  &  velodata, &
1674                  &  veldatqc)
1675      END IF
1676   END SUBROUTINE dia_obs_dealloc
1677
1678   SUBROUTINE ini_date( ddobsini )
1679      !!----------------------------------------------------------------------
1680      !!                    ***  ROUTINE ini_date  ***
1681      !!         
1682      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
1683      !!
1684      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
1685      !!
1686      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
1687      !!
1688      !! History :
1689      !!        !  06-03  (K. Mogensen)  Original code
1690      !!        !  06-05  (K. Mogensen)  Reformatted
1691      !!        !  06-10  (A. Weaver) Cleaning
1692      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1693      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1694      !!----------------------------------------------------------------------
1695      USE phycst, ONLY : &            ! Physical constants
1696         & rday
1697!      USE daymod, ONLY : &            ! Time variables
1698!         & nmonth_len           
1699      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1700         & rdt
1701
1702      IMPLICIT NONE
1703
1704      !! * Arguments
1705      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
1706
1707      !! * Local declarations
1708      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1709      INTEGER :: imon
1710      INTEGER :: iday
1711      INTEGER :: ihou
1712      INTEGER :: imin
1713      INTEGER :: imday         ! Number of days in month.
1714      REAL(KIND=wp) :: zdayfrc ! Fraction of day
1715
1716      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1717
1718      !!----------------------------------------------------------------------
1719      !! Initial date initialization (year, month, day, hour, minute)
1720      !! (This assumes that the initial date is for 00z))
1721      !!----------------------------------------------------------------------
1722      iyea =   ndate0 / 10000
1723      imon = ( ndate0 - iyea * 10000 ) / 100
1724      iday =   ndate0 - iyea * 10000 - imon * 100
1725      ihou = 0
1726      imin = 0
1727
1728      !!----------------------------------------------------------------------
1729      !! Compute number of days + number of hours + min since initial time
1730      !!----------------------------------------------------------------------
1731      iday = iday + ( nit000 -1 ) * rdt / rday
1732      zdayfrc = ( nit000 -1 ) * rdt / rday
1733      zdayfrc = zdayfrc - aint(zdayfrc)
1734      ihou = int( zdayfrc * 24 )
1735      imin = int( (zdayfrc * 24 - ihou) * 60 )
1736
1737      !!-----------------------------------------------------------------------
1738      !! Convert number of days (iday) into a real date
1739      !!----------------------------------------------------------------------
1740
1741      CALL calc_month_len( iyea, imonth_len )
1742     
1743      DO WHILE ( iday > imonth_len(imon) )
1744         iday = iday - imonth_len(imon)
1745         imon = imon + 1 
1746         IF ( imon > 12 ) THEN
1747            imon = 1
1748            iyea = iyea + 1
1749            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1750         ENDIF
1751      END DO
1752
1753      !!----------------------------------------------------------------------
1754      !! Convert it into YYYYMMDD.HHMMSS format.
1755      !!----------------------------------------------------------------------
1756      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1757         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1758
1759
1760   END SUBROUTINE ini_date
1761
1762   SUBROUTINE fin_date( ddobsfin )
1763      !!----------------------------------------------------------------------
1764      !!                    ***  ROUTINE fin_date  ***
1765      !!         
1766      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
1767      !!
1768      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
1769      !!
1770      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
1771      !!
1772      !! History :
1773      !!        !  06-03  (K. Mogensen)  Original code
1774      !!        !  06-05  (K. Mogensen)  Reformatted
1775      !!        !  06-10  (A. Weaver) Cleaning
1776      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1777      !!----------------------------------------------------------------------
1778      USE phycst, ONLY : &            ! Physical constants
1779         & rday
1780!      USE daymod, ONLY : &            ! Time variables
1781!         & nmonth_len               
1782      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1783         & rdt
1784
1785      IMPLICIT NONE
1786
1787      !! * Arguments
1788      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
1789
1790      !! * Local declarations
1791      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1792      INTEGER :: imon
1793      INTEGER :: iday
1794      INTEGER :: ihou
1795      INTEGER :: imin
1796      INTEGER :: imday         ! Number of days in month.
1797      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
1798         
1799      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1800           
1801      !-----------------------------------------------------------------------
1802      ! Initial date initialization (year, month, day, hour, minute)
1803      ! (This assumes that the initial date is for 00z)
1804      !-----------------------------------------------------------------------
1805      iyea =   ndate0 / 10000
1806      imon = ( ndate0 - iyea * 10000 ) / 100
1807      iday =   ndate0 - iyea * 10000 - imon * 100
1808      ihou = 0
1809      imin = 0
1810     
1811      !-----------------------------------------------------------------------
1812      ! Compute number of days + number of hours + min since initial time
1813      !-----------------------------------------------------------------------
1814      iday    = iday +  nitend  * rdt / rday
1815      zdayfrc =  nitend  * rdt / rday
1816      zdayfrc = zdayfrc - AINT( zdayfrc )
1817      ihou    = INT( zdayfrc * 24 )
1818      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1819
1820      !-----------------------------------------------------------------------
1821      ! Convert number of days (iday) into a real date
1822      !----------------------------------------------------------------------
1823
1824      CALL calc_month_len( iyea, imonth_len )
1825     
1826      DO WHILE ( iday > imonth_len(imon) )
1827         iday = iday - imonth_len(imon)
1828         imon = imon + 1 
1829         IF ( imon > 12 ) THEN
1830            imon = 1
1831            iyea = iyea + 1
1832            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1833         ENDIF
1834      END DO
1835
1836      !-----------------------------------------------------------------------
1837      ! Convert it into YYYYMMDD.HHMMSS format
1838      !-----------------------------------------------------------------------
1839      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1840         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1841
1842    END SUBROUTINE fin_date
1843   
1844END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.