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/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS/diaobs.F90 @ 14986

Last change on this file since 14986 was 14986, checked in by sparonuz, 3 years ago

Merge trunk -r14984:HEAD

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