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 @ 9125

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

Removed wrk_arrays from whole code. No change in SETTE results from this.

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