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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9041

Last change on this file since 9041 was 9041, checked in by timgraham, 6 years ago

Fix in diaobs

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