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 NEMO/trunk/src/OCE/OBS – NEMO

source: NEMO/trunk/src/OCE/OBS/diaobs.F90 @ 9656

Last change on this file since 9656 was 9656, checked in by clem, 6 years ago

remove the remaining references to LIM

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