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

source: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9003

Last change on this file since 9003 was 8992, checked in by timgraham, 7 years ago

Merged Met Office branch in

  • Property svn:keywords set to Id
File size: 48.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   !!   dia_obs_init : Reading and prepare observations
10   !!   dia_obs      : Compute model equivalent to observations
11   !!   dia_obs_wri  : Write observational diagnostics
12   !!   calc_date    : Compute the date of timestep in YYYYMMDD.HHMMSS format
13   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS
14   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
15   !!----------------------------------------------------------------------
16   !! * Modules used
17   USE wrk_nemo                 ! Memory Allocation
18   USE par_kind                 ! Precision variables
19   USE in_out_manager           ! I/O manager
20   USE par_oce
21   USE dom_oce                  ! Ocean space and time domain variables
22   USE obs_read_prof            ! Reading and allocation of profile obs
23   USE obs_read_surf            ! Reading and allocation of surface obs
24   USE obs_sstbias              ! Bias correction routine for SST
25   USE obs_readmdt              ! Reading and allocation of MDT for SLA.
26   USE obs_prep                 ! Preparation of obs. (grid search etc).
27   USE obs_oper                 ! Observation operators
28   USE obs_write                ! Writing of observation related diagnostics
29   USE obs_grid                 ! Grid searching
30   USE obs_read_altbias         ! Bias treatment for altimeter
31   USE obs_profiles_def         ! Profile data definitions
32   USE obs_surf_def             ! Surface data definitions
33   USE obs_types                ! Definitions for observation types
34   USE mpp_map                  ! MPP mapping
35   USE lib_mpp                  ! For ctl_warn/stop
36
37   IMPLICIT NONE
38
39   !! * Routine accessibility
40   PRIVATE
41   PUBLIC dia_obs_init, &  ! Initialize and read observations
42      &   dia_obs,      &  ! Compute model equivalent to observations
43      &   dia_obs_wri,  &  ! Write model equivalent to observations
44      &   dia_obs_dealloc, &  ! Deallocate dia_obs data
45      &   calc_date           ! Compute the date of a timestep
46
47   !! * Module variables
48   LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator
49   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs
50   LOGICAL :: ln_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres
51   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres
52   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres
53   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
54
55   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
56   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
57   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
58   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
59   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
60   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
61   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
62   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
63
64   INTEGER :: nn_1dint       !: Vertical interpolation method
65   INTEGER :: nn_2dint       !: Default horizontal interpolation method
66   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method
67   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method
68   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method
69   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method
70   INTEGER, DIMENSION(imaxavtypes) :: &
71      & nn_profdavtypes      !: Profile data types representing a daily average
72   INTEGER :: nproftypes     !: Number of profile obs types
73   INTEGER :: nsurftypes     !: Number of surface obs types
74   INTEGER, DIMENSION(:), ALLOCATABLE :: &
75      & nvarsprof, &         !: Number of profile variables
76      & nvarssurf            !: Number of surface variables
77   INTEGER, DIMENSION(:), ALLOCATABLE :: &
78      & nextrprof, &         !: Number of profile extra variables
79      & nextrsurf            !: Number of surface extra variables
80   INTEGER, DIMENSION(:), ALLOCATABLE :: &
81      & n2dintsurf           !: Interpolation option for surface variables
82   REAL(wp), DIMENSION(:), ALLOCATABLE :: &
83      & zavglamscl, &        !: E/W diameter of averaging footprint for surface variables
84      & zavgphiscl           !: N/S diameter of averaging footprint for surface variables
85   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
86      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres
87      & llnightav            !: Logical for calculating night-time averages
88
89   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
90      & surfdata, &          !: Initial surface data
91      & surfdataqc           !: Surface data after quality control
92   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
93      & profdata, &          !: Initial profile data
94      & profdataqc           !: Profile data after quality control
95
96   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
97      & cobstypesprof, &     !: Profile obs types
98      & cobstypessurf        !: Surface obs types
99
100   !!----------------------------------------------------------------------
101   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
102   !! $Id$
103   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
104   !!----------------------------------------------------------------------
105
106CONTAINS
107
108   SUBROUTINE dia_obs_init
109      !!----------------------------------------------------------------------
110      !!                    ***  ROUTINE dia_obs_init  ***
111      !!         
112      !! ** Purpose : Initialize and read observations
113      !!
114      !! ** Method  : Read the namelist and call reading routines
115      !!
116      !! ** Action  : Read the namelist and call reading routines
117      !!
118      !! History :
119      !!        !  06-03  (K. Mogensen) Original code
120      !!        !  06-05  (A. Weaver) Reformatted
121      !!        !  06-10  (A. Weaver) Cleaning and add controls
122      !!        !  07-03  (K. Mogensen) General handling of profiles
123      !!        !  14-08  (J.While) Incorporated SST bias correction 
124      !!        !  15-02  (M. Martin) Simplification of namelist and code
125      !!----------------------------------------------------------------------
126
127      IMPLICIT NONE
128
129      !! * Local declarations
130      INTEGER, PARAMETER :: &
131         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
132      INTEGER, DIMENSION(:), ALLOCATABLE :: &
133         & ifilesprof, &         ! Number of profile files
134         & ifilessurf            ! Number of surface files
135      INTEGER :: ios             ! Local integer output status for namelist read
136      INTEGER :: jtype           ! Counter for obs types
137      INTEGER :: jvar            ! Counter for variables
138      INTEGER :: jfile           ! Counter for files
139
140      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
141         & cn_profbfiles, &      ! T/S profile input filenames
142         & cn_sstfbfiles, &      ! Sea surface temperature input filenames
143         & cn_sssfbfiles, &      ! Sea surface salinity input filenames
144         & cn_slafbfiles, &      ! Sea level anomaly input filenames
145         & cn_sicfbfiles, &      ! Seaice concentration input filenames
146         & cn_velfbfiles, &      ! Velocity profile input filenames
147         & cn_sstbiasfiles      ! SST bias input filenames
148      CHARACTER(LEN=128) :: &
149         & cn_altbiasfile        ! Altimeter bias input filename
150      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
151         & clproffiles, &        ! Profile filenames
152         & clsurffiles           ! Surface filenames
153
154      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
155      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
156      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies
157      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
158      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity
159      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
160      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
161      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
162      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
163      LOGICAL :: ln_sstbias      ! Logical switch for bias corection of SST
164      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
165      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
166      LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs.
167      LOGICAL :: llvar1          ! Logical for profile variable 1
168      LOGICAL :: llvar2          ! Logical for profile variable 1
169      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files
170
171      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
172      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
173      REAL(wp), POINTER, DIMENSION(:,:) :: &
174         & zglam1, &             ! Model longitudes for profile variable 1
175         & zglam2                ! Model longitudes for profile variable 2
176      REAL(wp), POINTER, DIMENSION(:,:) :: &
177         & zgphi1, &             ! Model latitudes for profile variable 1
178         & zgphi2                ! Model latitudes for profile variable 2
179      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
180         & zmask1, &             ! Model land/sea mask associated with variable 1
181         & zmask2                ! Model land/sea mask associated with variable 2
182
183      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
184         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
185         &            ln_altbias, ln_sstbias, ln_nea,                 &
186         &            ln_grid_global, ln_grid_search_lookup,          &
187         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
188         &            ln_sstnight,                                    &
189         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
190         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
191         &            cn_profbfiles, cn_slafbfiles,                   &
192         &            cn_sstfbfiles, cn_sicfbfiles,                   &
193         &            cn_velfbfiles, cn_sssfbfiles,                   &
194         &            cn_sstbiasfiles, cn_altbiasfile,                &
195         &            cn_gridsearchfile, rn_gridsearchres,            &
196         &            rn_dobsini, rn_dobsend,                         &
197         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
198         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
199         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
200         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
201         &            nn_1dint, nn_2dint,                             &
202         &            nn_2dint_sla, nn_2dint_sst,                     &
203         &            nn_2dint_sss, nn_2dint_sic,                     &
204         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
205         &            nn_profdavtypes
206
207      INTEGER :: jnumsstbias
208      CALL wrk_alloc( jpi, jpj, zglam1 )
209      CALL wrk_alloc( jpi, jpj, zglam2 )
210      CALL wrk_alloc( jpi, jpj, zgphi1 )
211      CALL wrk_alloc( jpi, jpj, zgphi2 )
212      CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
213      CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
214
215      !-----------------------------------------------------------------------
216      ! Read namelist parameters
217      !-----------------------------------------------------------------------
218      ! Some namelist arrays need initialising
219      cn_profbfiles(:) = ''
220      cn_slafbfiles(:) = ''
221      cn_sstfbfiles(:) = ''
222      cn_sicfbfiles(:) = ''
223      cn_velfbfiles(:) = ''
224      cn_sssfbfiles(:)    = ''
225      cn_sstbiasfiles(:) = ''
226      nn_profdavtypes(:) = -1
227
228      CALL ini_date( rn_dobsini )
229      CALL fin_date( rn_dobsend )
230
231      ! Read namelist namobs : control observation diagnostics
232      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
233      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
234901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
235
236      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
237      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
238902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
239      IF(lwm) WRITE ( numond, namobs )
240
241      IF ( .NOT. ln_diaobs ) THEN
242         IF(lwp) WRITE(numout,cform_war)
243         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs'
244         RETURN
245      ENDIF
246
247      IF(lwp) THEN
248         WRITE(numout,*)
249         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
250         WRITE(numout,*) '~~~~~~~~~~~~'
251         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
252         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
253         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
254         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
255         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
256         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
257         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
258         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
259         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global
260         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
261         IF (ln_grid_search_lookup) &
262            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
263         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
264         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
265         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
266         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint
267         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
268         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
269         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
270         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
271         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
272         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
273         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias
274         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
275         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
276         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
277      ENDIF
278      !-----------------------------------------------------------------------
279      ! Set up list of observation types to be used
280      ! and the files associated with each type
281      !-----------------------------------------------------------------------
282
283      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
284      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) )
285
286      IF (ln_sstbias) THEN
287         lmask(:) = .FALSE. 
288         WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE. 
289         jnumsstbias = COUNT(lmask) 
290         lmask(:) = .FALSE. 
291      ENDIF     
292
293      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
294         IF(lwp) WRITE(numout,cform_war)
295         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
296            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &
297            &                    ' are set to .FALSE. so turning off calls to dia_obs'
298         nwarn = nwarn + 1
299         ln_diaobs = .FALSE.
300         RETURN
301      ENDIF
302
303      IF ( nproftypes > 0 ) THEN
304
305         ALLOCATE( cobstypesprof(nproftypes) )
306         ALLOCATE( ifilesprof(nproftypes) )
307         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
308
309         jtype = 0
310         IF (ln_t3d .OR. ln_s3d) THEN
311            jtype = jtype + 1
312            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', &
313               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
314         ENDIF
315         IF (ln_vel3d) THEN
316            jtype = jtype + 1
317            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
318               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
319         ENDIF
320
321      ENDIF
322
323      IF ( nsurftypes > 0 ) THEN
324
325         ALLOCATE( cobstypessurf(nsurftypes) )
326         ALLOCATE( ifilessurf(nsurftypes) )
327         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
328         ALLOCATE(n2dintsurf(nsurftypes))
329         ALLOCATE(zavglamscl(nsurftypes))
330         ALLOCATE(zavgphiscl(nsurftypes))
331         ALLOCATE(lfpindegs(nsurftypes))
332         ALLOCATE(llnightav(nsurftypes))
333
334         jtype = 0
335         IF (ln_sla) THEN
336            jtype = jtype + 1
337            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', &
338               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
339            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      &
340               &                  nn_2dint, nn_2dint_sla,             &
341               &                  rn_sla_avglamscl, rn_sla_avgphiscl, &
342               &                  ln_sla_fp_indegs, .FALSE.,          &
343               &                  n2dintsurf, zavglamscl, zavgphiscl, &
344               &                  lfpindegs, llnightav )
345         ENDIF
346         IF (ln_sst) THEN
347            jtype = jtype + 1
348            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', &
349               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
350            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      &
351               &                  nn_2dint, nn_2dint_sst,             &
352               &                  rn_sst_avglamscl, rn_sst_avgphiscl, &
353               &                  ln_sst_fp_indegs, ln_sstnight,      &
354               &                  n2dintsurf, zavglamscl, zavgphiscl, &
355               &                  lfpindegs, llnightav )
356         ENDIF
357#if defined key_lim2 || defined key_lim3 || defined key_cice
358         IF (ln_sic) THEN
359            jtype = jtype + 1
360            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', &
361               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
362            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      &
363               &                  nn_2dint, nn_2dint_sic,             &
364               &                  rn_sic_avglamscl, rn_sic_avgphiscl, &
365               &                  ln_sic_fp_indegs, .FALSE.,          &
366               &                  n2dintsurf, zavglamscl, zavgphiscl, &
367               &                  lfpindegs, llnightav )
368         ENDIF
369#endif
370         IF (ln_sss) THEN
371            jtype = jtype + 1
372            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', &
373               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
374            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      &
375               &                  nn_2dint, nn_2dint_sss,             &
376               &                  rn_sss_avglamscl, rn_sss_avgphiscl, &
377               &                  ln_sss_fp_indegs, .FALSE.,          &
378               &                  n2dintsurf, zavglamscl, zavgphiscl, &
379               &                  lfpindegs, llnightav )
380         ENDIF
381
382      ENDIF
383
384
385
386      !-----------------------------------------------------------------------
387      ! Obs operator parameter checking and initialisations
388      !-----------------------------------------------------------------------
389
390      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
391         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
392         RETURN
393      ENDIF
394
395      IF ( ln_grid_global ) THEN
396         CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
397      ENDIF
398
399      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
400         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
401            &                    ' is not available')
402      ENDIF
403
404      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN
405         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
406            &                    ' is not available')
407      ENDIF
408
409      CALL obs_typ_init
410      IF(ln_grid_global) THEN
411         CALL mppmap_init
412      ENDIF
413
414      CALL obs_grid_setup( )
415
416      !-----------------------------------------------------------------------
417      ! Depending on switches read the various observation types
418      !-----------------------------------------------------------------------
419
420      IF ( nproftypes > 0 ) THEN
421
422         ALLOCATE(profdata(nproftypes))
423         ALLOCATE(profdataqc(nproftypes))
424         ALLOCATE(nvarsprof(nproftypes))
425         ALLOCATE(nextrprof(nproftypes))
426
427         DO jtype = 1, nproftypes
428
429            nvarsprof(jtype) = 2
430            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
431               nextrprof(jtype) = 1
432               llvar1 = ln_t3d
433               llvar2 = ln_s3d
434               zglam1 = glamt
435               zgphi1 = gphit
436               zmask1 = tmask
437               zglam2 = glamt
438               zgphi2 = gphit
439               zmask2 = tmask
440            ENDIF
441            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
442               nextrprof(jtype) = 2
443               llvar1 = ln_vel3d
444               llvar2 = ln_vel3d
445               zglam1 = glamu
446               zgphi1 = gphiu
447               zmask1 = umask
448               zglam2 = glamv
449               zgphi2 = gphiv
450               zmask2 = vmask
451            ENDIF
452
453            !Read in profile or profile obs types
454            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
455               &               clproffiles(jtype,1:ifilesprof(jtype)), &
456               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
457               &               rn_dobsini, rn_dobsend, llvar1, llvar2, &
458               &               ln_ignmis, ln_s_at_t, .FALSE., &
459               &               kdailyavtypes = nn_profdavtypes )
460
461            DO jvar = 1, nvarsprof(jtype)
462               CALL obs_prof_staend( profdata(jtype), jvar )
463            END DO
464
465            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
466               &               llvar1, llvar2, &
467               &               jpi, jpj, jpk, &
468               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
469               &               ln_nea, ln_bound_reject, &
470               &               kdailyavtypes = nn_profdavtypes )
471
472         END DO
473
474         DEALLOCATE( ifilesprof, clproffiles )
475
476      ENDIF
477
478      IF ( nsurftypes > 0 ) THEN
479
480         ALLOCATE(surfdata(nsurftypes))
481         ALLOCATE(surfdataqc(nsurftypes))
482         ALLOCATE(nvarssurf(nsurftypes))
483         ALLOCATE(nextrsurf(nsurftypes))
484
485         DO jtype = 1, nsurftypes
486
487            nvarssurf(jtype) = 1
488            nextrsurf(jtype) = 0
489            llnightav(jtype) = .FALSE.
490            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
491            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight
492
493            !Read in surface obs types
494            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
495               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
496               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
497               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
498
499            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
500
501            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
502               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
503               IF ( ln_altbias ) &
504                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
505            ENDIF
506
507            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
508               jnumsstbias = 0
509               DO jfile = 1, jpmaxnfiles
510                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
511                     &  jnumsstbias = jnumsstbias + 1
512               END DO
513               IF ( jnumsstbias == 0 ) THEN
514                  CALL ctl_stop("ln_sstbias set but no bias files to read in")   
515               ENDIF
516
517               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
518                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
519
520            ENDIF
521         END DO
522
523         DEALLOCATE( ifilessurf, clsurffiles )
524
525      ENDIF
526
527      CALL wrk_dealloc( jpi, jpj, zglam1 )
528      CALL wrk_dealloc( jpi, jpj, zglam2 )
529      CALL wrk_dealloc( jpi, jpj, zgphi1 )
530      CALL wrk_dealloc( jpi, jpj, zgphi2 )
531      CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
532      CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
533
534   END SUBROUTINE dia_obs_init
535
536   SUBROUTINE dia_obs( kstp )
537      !!----------------------------------------------------------------------
538      !!                    ***  ROUTINE dia_obs  ***
539      !!         
540      !! ** Purpose : Call the observation operators on each time step
541      !!
542      !! ** Method  : Call the observation operators on each time step to
543      !!              compute the model equivalent of the following data:
544      !!               - Profile data, currently T/S or U/V
545      !!               - Surface data, currently SST, SLA or sea-ice concentration.
546      !!
547      !! ** Action  :
548      !!
549      !! History :
550      !!        !  06-03  (K. Mogensen) Original code
551      !!        !  06-05  (K. Mogensen) Reformatted
552      !!        !  06-10  (A. Weaver) Cleaning
553      !!        !  07-03  (K. Mogensen) General handling of profiles
554      !!        !  07-04  (G. Smith) Generalized surface operators
555      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
556      !!        !  14-08  (J. While) observation operator for profiles in
557      !!                             generalised vertical coordinates
558      !!        !  15-08  (M. Martin) Combined surface/profile routines.
559      !!----------------------------------------------------------------------
560      !! * Modules used
561      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
562         & gdept_n,       &     
563         & gdept_1d     
564      USE phycst, ONLY : &              ! Physical constants
565         & rday                         
566      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
567         & tsn,  &             
568         & un, vn, &
569         & sshn 
570      USE phycst, ONLY : &         ! Physical constants
571         & rday
572#if defined  key_lim3
573      USE ice, ONLY : &            ! LIM3 Ice model variables
574         & frld
575#endif
576#if defined key_lim2
577      USE ice_2, ONLY : &          ! LIM2 Ice model variables
578         & frld
579#endif
580#if defined key_cice
581      USE sbc_oce, ONLY : fr_i     ! ice fraction
582#endif
583
584      IMPLICIT NONE
585
586      !! * Arguments
587      INTEGER, INTENT(IN) :: kstp  ! Current timestep
588      !! * Local declarations
589      INTEGER :: idaystp           ! Number of timesteps per day
590      INTEGER :: jtype             ! Data loop variable
591      INTEGER :: jvar              ! Variable number
592      INTEGER :: ji, jj            ! Loop counters
593      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
594         & zprofvar1, &            ! Model values for 1st variable in a prof ob
595         & zprofvar2               ! Model values for 2nd variable in a prof ob
596      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
597         & zprofmask1, &           ! Mask associated with zprofvar1
598         & zprofmask2              ! Mask associated with zprofvar2
599      REAL(wp), POINTER, DIMENSION(:,:) :: &
600         & zsurfvar, &             ! Model values equivalent to surface ob.
601         & zsurfmask               ! Mask associated with surface variable
602      REAL(wp), POINTER, DIMENSION(:,:) :: &
603         & zglam1,    &            ! Model longitudes for prof variable 1
604         & zglam2,    &            ! Model longitudes for prof variable 2
605         & zgphi1,    &            ! Model latitudes for prof variable 1
606         & zgphi2                  ! Model latitudes for prof variable 2
607
608      !Allocate local work arrays
609      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
610      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
611      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
612      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
613      CALL wrk_alloc( jpi, jpj, zsurfvar )
614      CALL wrk_alloc( jpi, jpj, zsurfmask )
615      CALL wrk_alloc( jpi, jpj, zglam1 )
616      CALL wrk_alloc( jpi, jpj, zglam2 )
617      CALL wrk_alloc( jpi, jpj, zgphi1 )
618      CALL wrk_alloc( jpi, jpj, zgphi2 )
619
620      IF(lwp) THEN
621         WRITE(numout,*)
622         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
623         WRITE(numout,*) '~~~~~~~'
624      ENDIF
625
626      idaystp = NINT( rday / rdt )
627
628      !-----------------------------------------------------------------------
629      ! Call the profile and surface observation operators
630      !-----------------------------------------------------------------------
631
632      IF ( nproftypes > 0 ) THEN
633
634         DO jtype = 1, nproftypes
635
636            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
637            CASE('prof')
638               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
639               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
640               zprofmask1(:,:,:) = tmask(:,:,:)
641               zprofmask2(:,:,:) = tmask(:,:,:)
642               zglam1(:,:) = glamt(:,:)
643               zglam2(:,:) = glamt(:,:)
644               zgphi1(:,:) = gphit(:,:)
645               zgphi2(:,:) = gphit(:,:)
646            CASE('vel')
647               zprofvar1(:,:,:) = un(:,:,:)
648               zprofvar2(:,:,:) = vn(:,:,:)
649               zprofmask1(:,:,:) = umask(:,:,:)
650               zprofmask2(:,:,:) = vmask(:,:,:)
651               zglam1(:,:) = glamu(:,:)
652               zglam2(:,:) = glamv(:,:)
653               zgphi1(:,:) = gphiu(:,:)
654               zgphi2(:,:) = gphiv(:,:)
655            CASE DEFAULT
656               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
657            END SELECT
658
659            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
660               &               nit000, idaystp,                         &
661               &               zprofvar1, zprofvar2,                    &
662               &               gdept_n(:,:,:), gdepw_n(:,:,:),            & 
663               &               zprofmask1, zprofmask2,                  &
664               &               zglam1, zglam2, zgphi1, zgphi2,          &
665               &               nn_1dint, nn_2dint,                      &
666               &               kdailyavtypes = nn_profdavtypes )
667
668         END DO
669
670      ENDIF
671
672      IF ( nsurftypes > 0 ) THEN
673
674         DO jtype = 1, nsurftypes
675
676            !Defaults which might be changed
677            zsurfmask(:,:) = tmask(:,:,1)
678
679            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
680            CASE('sst')
681               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
682            CASE('sla')
683               zsurfvar(:,:) = sshn(:,:)
684            CASE('sss')
685               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
686            CASE('sic')
687               IF ( kstp == 0 ) THEN
688                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
689                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
690                        &           'time-step but some obs are valid then.' )
691                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
692                        &           ' sea-ice obs will be missed'
693                  ENDIF
694                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
695                     &                        surfdataqc(jtype)%nsstp(1)
696                  CYCLE
697               ELSE
698#if defined key_cice
699                  zsurfvar(:,:) = fr_i(:,:)
700#elif defined key_lim2 || defined key_lim3
701                  zsurfvar(:,:) = 1._wp - frld(:,:)
702#else
703               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
704                  &           ' but no sea-ice model appears to have been defined' )
705#endif
706               ENDIF
707
708            END SELECT
709
710            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
711               &               nit000, idaystp, zsurfvar, zsurfmask,    &
712               &               n2dintsurf(jtype), llnightav(jtype),     &
713               &               zavglamscl(jtype), zavgphiscl(jtype),     &
714               &               lfpindegs(jtype) )
715
716         END DO
717
718      ENDIF
719
720      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
721      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
722      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
723      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
724      CALL wrk_dealloc( jpi, jpj, zsurfvar )
725      CALL wrk_dealloc( jpi, jpj, zsurfmask )
726      CALL wrk_dealloc( jpi, jpj, zglam1 )
727      CALL wrk_dealloc( jpi, jpj, zglam2 )
728      CALL wrk_dealloc( jpi, jpj, zgphi1 )
729      CALL wrk_dealloc( jpi, jpj, zgphi2 )
730
731   END SUBROUTINE dia_obs
732
733   SUBROUTINE dia_obs_wri
734      !!----------------------------------------------------------------------
735      !!                    ***  ROUTINE dia_obs_wri  ***
736      !!         
737      !! ** Purpose : Call observation diagnostic output routines
738      !!
739      !! ** Method  : Call observation diagnostic output routines
740      !!
741      !! ** Action  :
742      !!
743      !! History :
744      !!        !  06-03  (K. Mogensen) Original code
745      !!        !  06-05  (K. Mogensen) Reformatted
746      !!        !  06-10  (A. Weaver) Cleaning
747      !!        !  07-03  (K. Mogensen) General handling of profiles
748      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
749      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
750      !!----------------------------------------------------------------------
751      !! * Modules used
752      USE obs_rot_vel          ! Rotation of velocities
753
754      IMPLICIT NONE
755
756      !! * Local declarations
757      INTEGER :: jtype                    ! Data set loop variable
758      INTEGER :: jo, jvar, jk
759      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
760         & zu, &
761         & zv
762
763      !-----------------------------------------------------------------------
764      ! Depending on switches call various observation output routines
765      !-----------------------------------------------------------------------
766
767      IF ( nproftypes > 0 ) THEN
768
769         DO jtype = 1, nproftypes
770
771            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
772
773               ! For velocity data, rotate the model velocities to N/S, E/W
774               ! using the compressed data structure.
775               ALLOCATE( &
776                  & zu(profdataqc(jtype)%nvprot(1)), &
777                  & zv(profdataqc(jtype)%nvprot(2))  &
778                  & )
779
780               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
781
782               DO jo = 1, profdataqc(jtype)%nprof
783                  DO jvar = 1, 2
784                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
785
786                        IF ( jvar == 1 ) THEN
787                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
788                        ELSE
789                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
790                        ENDIF
791
792                     END DO
793                  END DO
794               END DO
795
796               DEALLOCATE( zu )
797               DEALLOCATE( zv )
798
799            END IF
800
801            CALL obs_prof_decompress( profdataqc(jtype), &
802               &                      profdata(jtype), .TRUE., numout )
803
804            CALL obs_wri_prof( profdata(jtype) )
805
806         END DO
807
808      ENDIF
809
810      IF ( nsurftypes > 0 ) THEN
811
812         DO jtype = 1, nsurftypes
813
814            CALL obs_surf_decompress( surfdataqc(jtype), &
815               &                      surfdata(jtype), .TRUE., numout )
816
817            CALL obs_wri_surf( surfdata(jtype) )
818
819         END DO
820
821      ENDIF
822
823   END SUBROUTINE dia_obs_wri
824
825   SUBROUTINE dia_obs_dealloc
826      IMPLICIT NONE
827      !!----------------------------------------------------------------------
828      !!                    *** ROUTINE dia_obs_dealloc ***
829      !!
830      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
831      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
832      !!
833      !!  ** Method : Clean up various arrays left behind by the obs_oper.
834      !!
835      !!  ** Action :
836      !!
837      !!----------------------------------------------------------------------
838      ! obs_grid deallocation
839      CALL obs_grid_deallocate
840
841      ! diaobs deallocation
842      IF ( nproftypes > 0 ) &
843         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
844
845      IF ( nsurftypes > 0 ) &
846         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
847         &               n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav )
848
849   END SUBROUTINE dia_obs_dealloc
850
851   SUBROUTINE calc_date( kstp, ddobs )
852      !!----------------------------------------------------------------------
853      !!                    ***  ROUTINE calc_date  ***
854      !!         
855      !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format
856      !!
857      !! ** Method  : Get date in double precision YYYYMMDD.HHMMSS format
858      !!
859      !! ** Action  : Get date in double precision YYYYMMDD.HHMMSS format
860      !!
861      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
862      !!
863      !! History :
864      !!        !  06-03  (K. Mogensen)  Original code
865      !!        !  06-05  (K. Mogensen)  Reformatted
866      !!        !  06-10  (A. Weaver) Cleaning
867      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
868      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
869      !!        !  2014-09  (D. Lea) New generic routine now deals with arbitrary initial time of day
870      !!----------------------------------------------------------------------
871      USE phycst, ONLY : &            ! Physical constants
872         & rday
873      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
874         & rdt
875
876      IMPLICIT NONE
877
878      !! * Arguments
879      REAL(KIND=dp), INTENT(OUT) :: ddobs                        ! Date in YYYYMMDD.HHMMSS
880      INTEGER :: kstp
881
882      !! * Local declarations
883      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
884      INTEGER :: imon
885      INTEGER :: iday
886      INTEGER :: ihou
887      INTEGER :: imin
888      INTEGER :: imday       ! Number of days in month.
889      REAL(wp) :: zdayfrc    ! Fraction of day
890
891      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
892
893      !!----------------------------------------------------------------------
894      !! Initial date initialization (year, month, day, hour, minute)
895      !!----------------------------------------------------------------------
896      iyea =   ndate0 / 10000
897      imon = ( ndate0 - iyea * 10000 ) / 100
898      iday =   ndate0 - iyea * 10000 - imon * 100
899      ihou =   nn_time0 / 100
900      imin = ( nn_time0 - ihou * 100 ) 
901
902      !!----------------------------------------------------------------------
903      !! Compute number of days + number of hours + min since initial time
904      !!----------------------------------------------------------------------
905      zdayfrc = kstp * rdt / rday
906      zdayfrc = zdayfrc - aint(zdayfrc)
907      imin = imin + int( zdayfrc * 24 * 60 ) 
908      DO WHILE (imin >= 60) 
909        imin=imin-60
910        ihou=ihou+1
911      END DO
912      DO WHILE (ihou >= 24)
913        ihou=ihou-24
914        iday=iday+1
915      END DO
916      iday = iday + kstp * rdt / rday 
917
918      !-----------------------------------------------------------------------
919      ! Convert number of days (iday) into a real date
920      !----------------------------------------------------------------------
921
922      CALL calc_month_len( iyea, imonth_len )
923
924      DO WHILE ( iday > imonth_len(imon) )
925         iday = iday - imonth_len(imon)
926         imon = imon + 1 
927         IF ( imon > 12 ) THEN
928            imon = 1
929            iyea = iyea + 1
930            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
931         ENDIF
932      END DO
933
934      !----------------------------------------------------------------------
935      ! Convert it into YYYYMMDD.HHMMSS format.
936      !----------------------------------------------------------------------
937      ddobs = iyea * 10000_dp + imon * 100_dp + &
938         &    iday + ihou * 0.01_dp + imin * 0.0001_dp
939
940   END SUBROUTINE calc_date
941
942   SUBROUTINE ini_date( ddobsini )
943      !!----------------------------------------------------------------------
944      !!                    ***  ROUTINE ini_date  ***
945      !!         
946      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
947      !!
948      !! ** Method  :
949      !!
950      !! ** Action  :
951      !!
952      !! History :
953      !!        !  06-03  (K. Mogensen)  Original code
954      !!        !  06-05  (K. Mogensen)  Reformatted
955      !!        !  06-10  (A. Weaver) Cleaning
956      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
957      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
958      !!----------------------------------------------------------------------
959
960      IMPLICIT NONE
961
962      !! * Arguments
963      REAL(KIND=dp), INTENT(OUT) :: ddobsini                   ! Initial date in YYYYMMDD.HHMMSS
964
965      CALL calc_date( nit000 - 1, ddobsini )
966
967   END SUBROUTINE ini_date
968
969   SUBROUTINE fin_date( ddobsfin )
970      !!----------------------------------------------------------------------
971      !!                    ***  ROUTINE fin_date  ***
972      !!         
973      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
974      !!
975      !! ** Method  :
976      !!
977      !! ** Action  :
978      !!
979      !! History :
980      !!        !  06-03  (K. Mogensen)  Original code
981      !!        !  06-05  (K. Mogensen)  Reformatted
982      !!        !  06-10  (A. Weaver) Cleaning
983      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
984      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
985      !!----------------------------------------------------------------------
986
987      IMPLICIT NONE
988
989      !! * Arguments
990      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
991
992      CALL calc_date( nitend, ddobsfin )
993
994   END SUBROUTINE fin_date
995   
996    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
997       &                         cfilestype, ifiles, cobstypes, cfiles )
998
999    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
1000    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
1001    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs
1002    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1003       &                   ifiles      ! Out appended number of files for this type
1004
1005    CHARACTER(len=6), INTENT(IN) :: ctypein 
1006    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
1007       &                   cfilestype  ! In list of files for this obs type
1008    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
1009       &                   cobstypes   ! Out appended list of obs types
1010    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
1011       &                   cfiles      ! Out appended list of files for all types
1012
1013    !Local variables
1014    INTEGER :: jfile
1015
1016    cfiles(jtype,:) = cfilestype(:)
1017    cobstypes(jtype) = ctypein
1018    ifiles(jtype) = 0
1019    DO jfile = 1, jpmaxnfiles
1020       IF ( trim(cfiles(jtype,jfile)) /= '' ) &
1021                 ifiles(jtype) = ifiles(jtype) + 1
1022    END DO
1023
1024    IF ( ifiles(jtype) == 0 ) THEN
1025         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   &
1026            &           ' set to true but no files available to read' )
1027    ENDIF
1028
1029    IF(lwp) THEN   
1030       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
1031       DO jfile = 1, ifiles(jtype)
1032          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
1033       END DO
1034    ENDIF
1035
1036    END SUBROUTINE obs_settypefiles
1037
1038    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
1039               &                  n2dint_default, n2dint_type,        &
1040               &                  zavglamscl_type, zavgphiscl_type,   &
1041               &                  lfp_indegs_type, lavnight_type,     &
1042               &                  n2dint, zavglamscl, zavgphiscl,     &
1043               &                  lfpindegs, lavnight )
1044
1045    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
1046    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
1047    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
1048    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
1049    REAL(wp), INTENT(IN) :: &
1050       &                    zavglamscl_type, & !E/W diameter of obs footprint for this type
1051       &                    zavgphiscl_type    !N/S diameter of obs footprint for this type
1052    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
1053    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
1054    CHARACTER(len=6), INTENT(IN) :: ctypein 
1055
1056    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1057       &                    n2dint 
1058    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1059       &                    zavglamscl, zavgphiscl
1060    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1061       &                    lfpindegs, lavnight
1062
1063    lavnight(jtype) = lavnight_type
1064
1065    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
1066       n2dint(jtype) = n2dint_type
1067    ELSE
1068       n2dint(jtype) = n2dint_default
1069    ENDIF
1070
1071    ! For averaging observation footprints set options for size of footprint
1072    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
1073       IF ( zavglamscl_type > 0._wp ) THEN
1074          zavglamscl(jtype) = zavglamscl_type
1075       ELSE
1076          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1077                         'scale (zavglamscl) for observation type '//TRIM(ctypein) )     
1078       ENDIF
1079
1080       IF ( zavgphiscl_type > 0._wp ) THEN
1081          zavgphiscl(jtype) = zavgphiscl_type
1082       ELSE
1083          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1084                         'scale (zavgphiscl) for observation type '//TRIM(ctypein) )     
1085       ENDIF
1086
1087       lfpindegs(jtype) = lfp_indegs_type 
1088
1089    ENDIF
1090
1091    ! Write out info
1092    IF(lwp) THEN
1093       IF ( n2dint(jtype) <= 4 ) THEN
1094          WRITE(numout,*) '             '//TRIM(ctypein)// &
1095             &            ' model counterparts will be interpolated horizontally'
1096       ELSE IF ( n2dint(jtype) <= 6 ) THEN
1097          WRITE(numout,*) '             '//TRIM(ctypein)// &
1098             &            ' model counterparts will be averaged horizontally'
1099          WRITE(numout,*) '             '//'    with E/W scale: ',zavglamscl(jtype)
1100          WRITE(numout,*) '             '//'    with N/S scale: ',zavgphiscl(jtype)
1101          IF ( lfpindegs(jtype) ) THEN
1102              WRITE(numout,*) '             '//'    (in degrees)'
1103          ELSE
1104              WRITE(numout,*) '             '//'    (in metres)'
1105          ENDIF
1106       ENDIF
1107    ENDIF
1108
1109    END SUBROUTINE obs_setinterpopts
1110
1111END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.