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/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/diaobs.F90 @ 15395

Last change on this file since 15395 was 15395, checked in by dford, 3 years ago

Fix use of POTM/TEMP and grid type output.
File size: 50.2 KB
Line 
1MODULE diaobs
2   !!======================================================================
3   !!                       ***  MODULE diaobs  ***
4   !! Observation diagnostics: Computation of the misfit between data and
5   !!                          their model equivalent
6   !!======================================================================
7   !! 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 timing            ! Timing
31   USE par_oce           ! ocean parameter
32   USE dom_oce           ! Ocean space and time domain variables
33   USE sbc_oce           ! Sea-ice fraction
34   !
35   USE obs_read_prof     ! Reading and allocation of profile obs
36   USE obs_read_surf     ! Reading and allocation of surface obs
37   USE obs_bias          ! Bias correction routine
38   USE obs_readmdt       ! Reading and allocation of MDT for SLA.
39   USE obs_readsnowdepth ! Get model snow depth for conversion of freeboard to ice thickness
40   USE obs_prep          ! Preparation of obs. (grid search etc).
41   USE obs_oper          ! Observation operators
42   USE obs_write         ! Writing of observation related diagnostics
43   USE obs_grid          ! Grid searching
44   USE obs_read_altbias  ! Bias treatment for altimeter
45   USE obs_profiles_def  ! Profile data definitions
46   USE obs_surf_def      ! Surface data definitions
47   USE obs_types         ! Definitions for observation types
48   USE obs_group_def     ! Definitions for observation groups
49   !
50   USE mpp_map           ! MPP mapping
51   USE lib_mpp           ! For ctl_warn/stop
52
53   IMPLICIT NONE
54   PRIVATE
55
56   PUBLIC dia_obs_init     ! Initialize and read observations
57   PUBLIC dia_obs          ! Compute model equivalent to observations
58   PUBLIC dia_obs_wri      ! Write model equivalent to observations
59   PUBLIC calc_date        ! Compute the date of a timestep
60
61   LOGICAL, PUBLIC :: ln_diaobs            !: Logical switch for the obs operator
62   
63   INTEGER :: nn_obsgroups
64
65   TYPE(obs_group), DIMENSION(:), ALLOCATABLE ::   sobsgroups   ! Obs groups
66
67   !!----------------------------------------------------------------------
68   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
69   !! $Id$
70   !! Software governed by the CeCILL license (see ./LICENSE)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE dia_obs_init
75      !!----------------------------------------------------------------------
76      !!                    ***  ROUTINE dia_obs_init  ***
77      !!         
78      !! ** Purpose : Initialize and read observations
79      !!
80      !! ** Method  : Read the namelist and call reading routines
81      !!
82      !! ** Action  : Read the namelist and call reading routines
83      !!
84      !!----------------------------------------------------------------------
85#if defined key_si3
86      USE ice, ONLY : &     ! Sea ice variables
87         & hm_s             ! Snow depth for freeboard conversion
88#elif defined key_cice
89      USE sbc_oce, ONLY : & ! Sea ice variables
90         & thick_s          ! Snow depth for freeboard conversion
91#endif
92      USE obs_fbm, ONLY : &
93         & fbrmdi           ! Real missing data indicator
94
95      IMPLICIT NONE
96
97      INTEGER, PARAMETER ::   jpmaxngroups = 1000    ! Maximum number of obs groups
98      INTEGER :: ios             ! Local integer output status for namelist read
99      INTEGER :: jtype           ! Counter for obs types
100      INTEGER :: jvar            ! Counter for variables
101      INTEGER :: jfile           ! Counter for files
102      INTEGER :: jenabled
103      INTEGER :: jgroup
104      !
105      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar
106      !
107      REAL(dp) :: rn_dobsini      ! Obs window start date YYYYMMDD.HHMMSS
108      REAL(dp) :: rn_dobsend      ! Obs window end date   YYYYMMDD.HHMMSS
109      !!
110      NAMELIST/namobs/ln_diaobs, nn_obsgroups,               &
111         &            ln_grid_global, ln_grid_search_lookup, &
112         &            cn_gridsearchfile, rn_gridsearchres,   &
113         &            rn_dobsini, rn_dobsend
114      !-----------------------------------------------------------------------
115
116      !-----------------------------------------------------------------------
117      ! Read namelist parameters
118      !-----------------------------------------------------------------------
119      ! Initialise time window
120      CALL ini_date( rn_dobsini )
121      CALL fin_date( rn_dobsend )
122
123      ! Read namelist namobs : control observation diagnostics
124      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
125      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
126901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namobs in reference namelist' )
127      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
128      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
129902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namobs in configuration namelist' )
130      IF(lwm) WRITE ( numond, namobs )
131
132      IF( .NOT.ln_diaobs ) THEN
133         IF(lwp) WRITE(numout,*)
134         IF(lwp) WRITE(numout,*) 'dia_obs_init : NO Observation diagnostic used'
135         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
136         RETURN
137      ENDIF
138
139      IF(lwp) THEN
140         WRITE(numout,*)
141         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
142         WRITE(numout,*) '~~~~~~~~~~~~'
143         WRITE(numout,*) '   Namelist namobs : set observation diagnostic parameters'
144         WRITE(numout,*) '      Number of namobs_dta namelists to read             nn_obsgroups = ', nn_obsgroups
145         WRITE(numout,*) '      Global distribution of observations              ln_grid_global = ', ln_grid_global
146         WRITE(numout,*) '      Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
147         IF (ln_grid_search_lookup) THEN
148            WRITE(numout,*) '      Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
149            WRITE(numout,*) '      Grid search resolution                         rn_gridsearchres = ', rn_gridsearchres
150         ENDIF
151      ENDIF
152
153      IF( ln_grid_global ) THEN
154         CALL ctl_warn( 'dia_obs_init: ln_grid_global=T may cause memory issues when used with a large number of processors' )
155      ENDIF
156
157      !-----------------------------------------------------------------------
158      ! Read namobs_dta namelists and set up observation groups
159      !-----------------------------------------------------------------------
160
161      IF( nn_obsgroups == 0 ) THEN
162         CALL ctl_warn( 'dia_obs_init: ln_diaobs is set to true, but nn_obsgroups == 0',   &
163            &           ' so turning off calls to dia_obs'  )
164         ln_diaobs = .FALSE.
165         RETURN
166      ENDIF
167
168      ALLOCATE( sobsgroups(nn_obsgroups) )
169
170      jenabled = 0
171      DO jgroup = 1, nn_obsgroups
172         CALL obs_group_read ( sobsgroups(jgroup) )
173         CALL obs_group_check( sobsgroups(jgroup), jgroup )
174         IF (sobsgroups(jgroup)%lenabled) THEN
175            jenabled = jenabled + 1
176            IF( sobsgroups(jgroup)%lvel  .AND.  .NOT.ln_grid_global ) THEN
177               CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
178            ENDIF
179         ENDIF
180      END DO
181
182      IF( jenabled == 0 ) THEN
183         CALL ctl_warn( 'dia_obs_init: ln_diaobs is set to true, and nn_obsgroups > 0',   &
184            &           ' but no groups are enabled so turning off calls to dia_obs'  )
185         ln_diaobs = .FALSE.
186         RETURN
187      ENDIF
188
189
190      !-----------------------------------------------------------------------
191      ! Obs operator parameter checking and initialisations
192      !-----------------------------------------------------------------------
193      !
194      CALL obs_typ_init
195      IF( ln_grid_global )   CALL mppmap_init
196      !
197      CALL obs_grid_setup( )
198
199      !-----------------------------------------------------------------------
200      ! Depending on switches read the various observation types
201      !-----------------------------------------------------------------------
202      !
203      DO jgroup = 1, nn_obsgroups
204         IF (sobsgroups(jgroup)%lenabled) THEN
205            IF (sobsgroups(jgroup)%lprof) THEN
206               !
207               ! Read in profile or profile obs types
208               !
209               ALLOCATE( llvar(sobsgroups(jgroup)%nobstypes) )
210               llvar(:) = .TRUE.
211               !
212               CALL obs_rea_prof( sobsgroups(jgroup)%sprofdata,   &
213                  &               sobsgroups(jgroup)%nobsfiles,   &
214                  &               sobsgroups(jgroup)%cobsfiles,   &
215                  &               sobsgroups(jgroup)%nobstypes,   &
216                  &               sobsgroups(jgroup)%naddvars,    &
217                  &               sobsgroups(jgroup)%nextvars,    &
218                  &               nitend-nit000+2,                &
219                  &               rn_dobsini,                     &
220                  &               rn_dobsend,                     &
221                  &               llvar,                          &
222                  &               sobsgroups(jgroup)%lignmis,     &
223                  &               sobsgroups(jgroup)%lall_at_all, &
224                  &               .FALSE.,                        &
225                  &               sobsgroups(jgroup)%cobstypes,   &
226                  &               kdailyavtypes = sobsgroups(jgroup)%nprofdavtypes )
227               !
228               DO jvar = 1, sobsgroups(jgroup)%nobstypes
229                  CALL obs_prof_staend( sobsgroups(jgroup)%sprofdata, jvar )
230               END DO
231               !
232               IF ( sobsgroups(jgroup)%sprofdata%next > 0 ) THEN
233                  CALL obs_prof_staend_ext( sobsgroups(jgroup)%sprofdata )
234               ENDIF
235               !
236               IF( sobsgroups(jgroup)%loutput_clim ) THEN
237                  sobsgroups(jgroup)%sprofdata%caddvars(sobsgroups(jgroup)%nadd_clm)  = 'CLM'
238                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
239                     sobsgroups(jgroup)%sprofdata%var(jvar)%vadd(:,sobsgroups(jgroup)%nadd_clm) = fbrmdi
240                     sobsgroups(jgroup)%sprofdata%caddlong(sobsgroups(jgroup)%nadd_clm,jvar) = 'Climatology'
241                     sobsgroups(jgroup)%sprofdata%caddunit(sobsgroups(jgroup)%nadd_clm,jvar) = sobsgroups(jgroup)%sprofdata%cunit(jvar)
242                  END DO
243               ENDIF
244               !
245               sobsgroups(jgroup)%sprofdata%cgrid = sobsgroups(jgroup)%cgrid
246               !
247               CALL obs_pre_prof( sobsgroups(jgroup)%sprofdata,     &
248                  &               sobsgroups(jgroup)%sprofdataqc,   &
249                  &               llvar,                            &
250                  &               jpi, jpj, jpk,                    &
251                  &               sobsgroups(jgroup)%rmask,         &
252                  &               sobsgroups(jgroup)%rglam,         &
253                  &               sobsgroups(jgroup)%rgphi,         &
254                  &               sobsgroups(jgroup)%lnea,          &
255                  &               sobsgroups(jgroup)%lbound_reject, &
256                  &               kdailyavtypes = sobsgroups(jgroup)%nprofdavtypes )
257               !
258               DEALLOCATE( llvar )
259               !
260            ELSEIF (sobsgroups(jgroup)%lsurf) THEN
261               !
262               ! Read in surface obs types
263               !
264               CALL obs_rea_surf( sobsgroups(jgroup)%ssurfdata,         &
265                  &               sobsgroups(jgroup)%nobsfiles,         &
266                  &               sobsgroups(jgroup)%cobsfiles,         &
267                  &               sobsgroups(jgroup)%nobstypes,         &
268                  &               sobsgroups(jgroup)%naddvars,          &
269                  &               sobsgroups(jgroup)%nextvars,          &
270                  &               nitend-nit000+2,                      &
271                  &               rn_dobsini,                           &
272                  &               rn_dobsend,                           &
273                  &               sobsgroups(jgroup)%rtime_mean_period, &
274                  &               sobsgroups(jgroup)%ltime_mean_bkg,    &
275                  &               sobsgroups(jgroup)%lignmis,           &
276                  &               .FALSE.,                              &
277                  &               sobsgroups(jgroup)%lnight,            &
278                  &               sobsgroups(jgroup)%cobstypes )
279               !
280               IF( sobsgroups(jgroup)%lsla ) THEN
281                  sobsgroups(jgroup)%ssurfdata%cextvars(sobsgroups(jgroup)%next_mdt) = 'MDT'
282                  sobsgroups(jgroup)%ssurfdata%cextlong(sobsgroups(jgroup)%next_mdt) = 'Mean dynamic topography'
283                  sobsgroups(jgroup)%ssurfdata%cextunit(sobsgroups(jgroup)%next_mdt) = 'Metres'
284                  sobsgroups(jgroup)%ssurfdata%caddvars(sobsgroups(jgroup)%nadd_ssh) = 'SSH'
285                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
286                     sobsgroups(jgroup)%ssurfdata%caddlong(sobsgroups(jgroup)%nadd_ssh,jvar) = 'Model Sea surface height'
287                     sobsgroups(jgroup)%ssurfdata%caddunit(sobsgroups(jgroup)%nadd_ssh,jvar) = 'Metres'
288                  END DO
289               ENDIF
290               !
291               IF( sobsgroups(jgroup)%lfbd ) THEN
292                  sobsgroups(jgroup)%ssurfdata%cextvars(sobsgroups(jgroup)%next_snow) = 'SNOW'
293                  sobsgroups(jgroup)%ssurfdata%cextlong(sobsgroups(jgroup)%next_snow) = 'Snow thickness'
294                  sobsgroups(jgroup)%ssurfdata%cextunit(sobsgroups(jgroup)%next_snow) = 'Metres'
295                  sobsgroups(jgroup)%ssurfdata%caddvars(sobsgroups(jgroup)%nadd_fbd)  = 'FBD'
296                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
297                     sobsgroups(jgroup)%ssurfdata%caddlong(sobsgroups(jgroup)%nadd_fbd,jvar) = 'Freeboard'
298                     sobsgroups(jgroup)%ssurfdata%caddunit(sobsgroups(jgroup)%nadd_fbd,jvar) = 'Metres'
299                  END DO
300               ENDIF
301               !
302               IF( sobsgroups(jgroup)%loutput_clim ) THEN
303                  sobsgroups(jgroup)%ssurfdata%caddvars(sobsgroups(jgroup)%nadd_clm)  = 'CLM'
304                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
305                     sobsgroups(jgroup)%ssurfdata%radd(:,:,jvar) = fbrmdi
306                     sobsgroups(jgroup)%ssurfdata%caddlong(sobsgroups(jgroup)%nadd_clm,jvar) = 'Climatology'
307                     sobsgroups(jgroup)%ssurfdata%caddunit(sobsgroups(jgroup)%nadd_clm,jvar) = sobsgroups(jgroup)%ssurfdata%cunit(jvar)
308                  END DO
309               ENDIF
310               !
311               sobsgroups(jgroup)%ssurfdata%cgrid = sobsgroups(jgroup)%cgrid
312               !
313               CALL obs_pre_surf( sobsgroups(jgroup)%ssurfdata,      &
314                  &               sobsgroups(jgroup)%ssurfdataqc,    &
315                  &               jpi, jpj,                          &
316                  &               sobsgroups(jgroup)%rmask(:,:,1,:), &
317                  &               sobsgroups(jgroup)%rglam,          &
318                  &               sobsgroups(jgroup)%rgphi,          &
319                  &               sobsgroups(jgroup)%lnea,           &
320                  &               sobsgroups(jgroup)%lbound_reject )
321               !
322               IF( sobsgroups(jgroup)%lsla ) THEN
323                  CALL obs_rea_mdt( sobsgroups(jgroup)%ssurfdataqc, &
324                     &              sobsgroups(jgroup)%n2dint,      &
325                     &              sobsgroups(jgroup)%next_mdt,    &
326                     &              sobsgroups(jgroup)%nmsshc,      &
327                     &              sobsgroups(jgroup)%rmdtcorr,    &
328                     &              sobsgroups(jgroup)%rmdtcutoff )
329                  IF( sobsgroups(jgroup)%laltbias ) THEN
330                     !CALL obs_rea_altbias( sobsgroups(jgroup)%ssurfdataqc, &
331                     !   &                  sobsgroups(jgroup)%n2dint,      &
332                     !   &                  sobsgroups(jgroup)%caltbiasfile )
333                     CALL obs_app_bias( sobsgroups(jgroup)%ssurfdataqc,   &
334                        &               sobsgroups(jgroup)%next_mdt,      & 
335                        &               sobsgroups(jgroup)%n2dint,        & 
336                        &               1,                                &
337                        &               sobsgroups(jgroup)%caltbiasfile,  &
338                        &               'altbias',                        &
339                        &               ld_extvar=.TRUE. ) 
340                  ENDIF
341               ENDIF
342               !
343#if defined key_si3
344               IF( sobsgroups(jgroup)%lfbd ) THEN
345                  CALL obs_rea_snowdepth( sobsgroups(jgroup)%ssurfdataqc, &
346                     &                    sobsgroups(jgroup)%n2dint,      &
347                     &                    sobsgroups(jgroup)%next_snow,   &
348                     &                    hm_s(:,:) )
349               ENDIF
350#elif defined key_cice
351               IF( sobsgroups(jgroup)%lfbd ) THEN
352                  CALL obs_rea_snowdepth( sobsgroups(jgroup)%ssurfdataqc, &
353                     &                    sobsgroups(jgroup)%n2dint,      &
354                     &                    sobsgroups(jgroup)%next_snow,   &
355                     &                    thick_s(:,:) )
356               ENDIF
357#endif
358               !
359               IF( sobsgroups(jgroup)%lobsbias ) THEN
360                  CALL obs_app_bias( sobsgroups(jgroup)%ssurfdataqc,   &
361                     &               sobsgroups(jgroup)%nbiasvar,      & 
362                     &               sobsgroups(jgroup)%n2dint,        & 
363                     &               sobsgroups(jgroup)%nobsbiasfiles, &
364                     &               sobsgroups(jgroup)%cobsbiasfiles, &
365                     &               sobsgroups(jgroup)%cbiasvarname ) 
366               ENDIF
367               !
368            ENDIF
369         ENDIF
370      END DO
371      !
372   END SUBROUTINE dia_obs_init
373
374
375   SUBROUTINE dia_obs( kstp )
376      !!----------------------------------------------------------------------
377      !!                    ***  ROUTINE dia_obs  ***
378      !!         
379      !! ** Purpose : Call the observation operators on each time step
380      !!
381      !! ** Method  : Call the observation operators on each time step to
382      !!              compute the model equivalent of the following data:
383      !!               - Profile data, currently T/S or U/V
384      !!               - Surface data, currently SST, SLA or sea-ice concentration.
385      !!
386      !! ** Action  :
387      !!----------------------------------------------------------------------
388      USE dom_oce, ONLY : gdept_n, gdept_1d   ! Ocean space and time domain variables
389      USE phycst , ONLY : rday                ! Physical constants
390      USE oce    , ONLY : tsn, un, vn, sshn   ! Ocean dynamics and tracers variables
391      USE phycst , ONLY : rday                ! Physical constants
392#if defined key_si3
393      USE ice    , ONLY : at_i, hm_i          ! SI3 Ice model variables
394#elif defined key_cice
395      USE sbc_oce, ONLY : fr_i, thick_i       ! CICE Ice model variables
396#endif
397      USE tradmp,  ONLY : tclim, sclim        ! T&S climatology
398      USE obs_fbm, ONLY : fbrmdi              ! Real missing data indicator
399
400      IMPLICIT NONE
401
402      !! * Arguments
403      INTEGER, INTENT(IN) :: kstp  ! Current timestep
404      !! * Local declarations
405      INTEGER :: idaystp           ! Number of timesteps per day
406      INTEGER :: imeanstp          ! Number of timesteps for time averaging
407      INTEGER :: jtype             ! Data loop variable
408      INTEGER :: jvar              ! Variable number
409      INTEGER :: jgroup
410      INTEGER :: ji, jj, jobs      ! Loop counters
411      LOGICAL :: lstp0             ! Flag special treatment on zeroth time step
412      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
413         & zprofvar, &             ! Model values for variables in a prof ob
414         & zprofclim               ! Climatology values for variables in a prof ob
415      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
416         & zsurfvar, &             ! Model values for variables in a surf ob
417         & zsurfclim               ! Climatology values for variables in a surf ob
418
419      !-----------------------------------------------------------------------
420
421      IF( ln_timing )   CALL timing_start('dia_obs')
422
423      IF(lwp) THEN
424         WRITE(numout,*)
425         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
426         WRITE(numout,*) '~~~~~~~'
427      ENDIF
428
429      idaystp = NINT( rday / rdt )
430
431      !-----------------------------------------------------------------------
432      ! Call the profile and surface observation operators
433      !-----------------------------------------------------------------------
434
435      ALLOCATE( zprofvar(jpi,jpj,jpk),  &
436         &      zprofclim(jpi,jpj,jpk), &
437         &      zsurfvar(jpi,jpj),      &
438         &      zsurfclim(jpi,jpj) )
439
440      DO jgroup = 1, nn_obsgroups
441         IF (sobsgroups(jgroup)%lenabled) THEN
442
443            IF (sobsgroups(jgroup)%lprof) THEN
444
445               zprofclim(:,:,:) = fbrmdi
446
447               DO jvar = 1, sobsgroups(jgroup)%nobstypes
448
449                  SELECT CASE ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) )
450                  CASE('POTM')
451                     zprofvar(:,:,:) = tsn(:,:,:,jp_tem)
452                     IF (sobsgroups(jgroup)%loutput_clim) THEN
453                        zprofclim(:,:,:) = tclim(:,:,:)
454                     ENDIF
455                  CASE('PSAL')
456                     zprofvar(:,:,:) = tsn(:,:,:,jp_sal)
457                     IF (sobsgroups(jgroup)%loutput_clim) THEN
458                        zprofclim(:,:,:) = sclim(:,:,:)
459                     ENDIF
460                  CASE('UVEL')
461                     zprofvar(:,:,:) = un(:,:,:)
462                  CASE('VVEL')
463                     zprofvar(:,:,:) = vn(:,:,:)
464                  CASE DEFAULT
465                     CALL ctl_stop( 'Unknown profile observation type '//TRIM(sobsgroups(jgroup)%cobstypes(jvar))//' in dia_obs' )
466                  END SELECT
467
468                  CALL obs_prof_opt( sobsgroups(jgroup)%sprofdataqc,       &
469                     &               kstp, jpi, jpj, jpk,                  &
470                     &               nit000, idaystp, jvar,                &
471                     &               zprofvar,                             &
472                     &               sobsgroups(jgroup)%loutput_clim,      &
473                     &               sobsgroups(jgroup)%nadd_clm,          &
474                     &               zprofclim,                            &
475                     &               gdept_n,                              &
476                     &               gdepw_n,                              & 
477                     &               sobsgroups(jgroup)%rmask(:,:,:,jvar), &
478                     &               sobsgroups(jgroup)%rglam(:,:,jvar),   &
479                     &               sobsgroups(jgroup)%rgphi(:,:,jvar),   &
480                     &               sobsgroups(jgroup)%n1dint,            &
481                     &               sobsgroups(jgroup)%n2dint,            &
482                     &               kdailyavtypes = sobsgroups(jgroup)%nprofdavtypes )
483
484               END DO
485
486            ELSEIF (sobsgroups(jgroup)%lsurf) THEN
487
488               zsurfclim(:,:) = fbrmdi
489
490               DO jvar = 1, sobsgroups(jgroup)%nobstypes
491
492                  lstp0 = .FALSE.
493                  SELECT CASE ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) )
494                  CASE('SST')
495                     zsurfvar(:,:) = tsn(:,:,1,jp_tem)
496                     IF (sobsgroups(jgroup)%loutput_clim) THEN
497                        zsurfclim(:,:) = tclim(:,:,1)
498                     ENDIF
499                  CASE('SLA')
500                     zsurfvar(:,:) = sshn(:,:)
501                  CASE('SSS')
502                     zsurfvar(:,:) = tsn(:,:,1,jp_sal)
503                     IF (sobsgroups(jgroup)%loutput_clim) THEN
504                        zsurfclim(:,:) = sclim(:,:,1)
505                     ENDIF
506                  CASE('UVEL')
507                     zsurfvar(:,:) = un(:,:,1)
508                  CASE('VVEL')
509                     zsurfvar(:,:) = vn(:,:,1)
510                  CASE('ICECONC')
511                     IF ( kstp == 0 ) THEN
512                        lstp0 = .TRUE.
513                     ELSE
514#if defined key_si3
515                        zsurfvar(:,:) = at_i(:,:)
516#elif defined key_cice
517                        zsurfvar(:,:) = fr_i(:,:)
518#else
519                        CALL ctl_stop( ' Trying to run sea-ice observation operator', &
520                           &           ' but no sea-ice model appears to have been defined' )
521#endif
522                     ENDIF
523                  CASE('SIT','FBD')
524                     IF ( kstp == 0 ) THEN
525                        lstp0 = .TRUE.
526                     ELSE
527#if defined key_si3
528                        zsurfvar(:,:) = hm_i(:,:)
529#elif defined key_cice
530                        zsurfvar(:,:) = thick_i(:,:)
531#else
532                        CALL ctl_stop( ' Trying to run sea-ice observation operator', &
533                           &           ' but no sea-ice model appears to have been defined' )
534#endif
535                     ENDIF
536                  END SELECT
537
538                  IF ( lstp0 ) THEN
539                     IF ( sobsgroups(jgroup)%ssurfdataqc%nsstpmpp(1) > 0 ) THEN
540                        DO jobs = sobsgroups(jgroup)%ssurfdataqc%nsurfup + 1, &
541                           &      sobsgroups(jgroup)%ssurfdataqc%nsurfup + sobsgroups(jgroup)%ssurfdataqc%nsstp(1)
542                           sobsgroups(jgroup)%ssurfdata%nqc(jobs) = IBSET(sobsgroups(jgroup)%ssurfdata%nqc(jobs),13)
543                        END DO
544                        IF ( lwp ) THEN
545                           CALL ctl_warn( TRIM(sobsgroups(jgroup)%cobstypes(jvar))// &
546                              &           ' not initialised on zeroth '           // &
547                              &           'time-step but some obs are valid then.' )
548                           WRITE(numout,*)sobsgroups(jgroup)%ssurfdataqc%nsstpmpp(1), &
549                              &           TRIM(sobsgroups(jgroup)%cobstypes(jvar)),   &
550                              &           'observations will be flagged as bad'
551                        ENDIF
552                     ENDIF
553                     IF ( jvar == sobsgroups(jgroup)%ssurfdataqc%nvar ) THEN
554                        sobsgroups(jgroup)%ssurfdataqc%nsurfup = sobsgroups(jgroup)%ssurfdataqc%nsurfup + &
555                           &                                     sobsgroups(jgroup)%ssurfdataqc%nsstp(1)
556                     ENDIF
557                  ELSE
558                     IF ( sobsgroups(jgroup)%ltime_mean_bkg ) THEN
559                        ! Number of time-steps in meaning period
560                        imeanstp = NINT( ( sobsgroups(jgroup)%rtime_mean_period * 60.0 * 60.0 ) / rdt )
561                     ENDIF
562                     CALL obs_surf_opt( sobsgroups(jgroup)%ssurfdataqc,       &
563                        &               kstp, jpi, jpj,                       &
564                        &               nit000, idaystp,                      &
565                        &               jvar, zsurfvar,                       &
566                        &               sobsgroups(jgroup)%loutput_clim,      &
567                        &               sobsgroups(jgroup)%nadd_clm,          &
568                        &               zsurfclim,                            &
569                        &               sobsgroups(jgroup)%rmask(:,:,1,jvar), &
570                        &               sobsgroups(jgroup)%n2dint,            &
571                        &               sobsgroups(jgroup)%lnight,            &
572                        &               sobsgroups(jgroup)%ravglamscl,        &
573                        &               sobsgroups(jgroup)%ravgphiscl,        &
574                        &               sobsgroups(jgroup)%lfp_indegs,        &
575                        &               sobsgroups(jgroup)%ltime_mean_bkg,    &
576                        &               imeanstp,                             &
577                        &               kssh=sobsgroups(jgroup)%nadd_ssh,     &
578                        &               kmdt=sobsgroups(jgroup)%next_mdt,     &
579                        &               kfbd=sobsgroups(jgroup)%nadd_fbd,     &
580                        &               ksnow=sobsgroups(jgroup)%next_snow )
581                  ENDIF
582
583               END DO
584
585            ENDIF
586
587         ENDIF
588      END DO
589
590      DEALLOCATE( zprofvar, zprofclim, &
591         &        zsurfvar, zsurfclim )
592
593      IF( ln_timing )   CALL timing_stop('dia_obs')
594
595   END SUBROUTINE dia_obs
596
597   SUBROUTINE dia_obs_wri
598      !!----------------------------------------------------------------------
599      !!                    ***  ROUTINE dia_obs_wri  ***
600      !!         
601      !! ** Purpose : Call observation diagnostic output routines
602      !!
603      !! ** Method  : Call observation diagnostic output routines
604      !!
605      !! ** Action  :
606      !!
607      !! History :
608      !!        !  06-03  (K. Mogensen) Original code
609      !!        !  06-05  (K. Mogensen) Reformatted
610      !!        !  06-10  (A. Weaver) Cleaning
611      !!        !  07-03  (K. Mogensen) General handling of profiles
612      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
613      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
614      !!----------------------------------------------------------------------
615      !! * Modules used
616      USE obs_rot_vel          ! Rotation of velocities
617
618      IMPLICIT NONE
619
620      !! * Local declarations
621      INTEGER :: jgroup                   ! Data set loop variable
622      INTEGER :: jo, jvar, jk, jadd, jext, jadd2, jext2
623      INTEGER :: iuvar, ivvar
624      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
625         & zu, &
626         & zv
627      LOGICAL, DIMENSION(:), ALLOCATABLE :: ll_write
628      TYPE(obswriinfo) :: sladd, slext
629
630      IF( ln_timing )   CALL timing_start('dia_obs_wri')
631
632      !-----------------------------------------------------------------------
633      ! Depending on switches call various observation output routines
634      !-----------------------------------------------------------------------
635
636      DO jgroup = 1, nn_obsgroups
637         IF (sobsgroups(jgroup)%lenabled) THEN
638
639            IF (sobsgroups(jgroup)%lprof) THEN
640
641               IF (sobsgroups(jgroup)%lvel) THEN
642                  iuvar = 0
643                  ivvar = 0
644                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
645                     IF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_uvel ) THEN
646                        iuvar = jvar
647                     ELSEIF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_vvel ) THEN
648                        ivvar = jvar
649                     ENDIF
650                  END DO
651                  IF ( (iuvar > 0) .AND. (ivvar > 0) ) THEN
652
653                     ! For velocity data, rotate the model velocities to N/S, E/W
654                     ! using the compressed data structure.
655                     ALLOCATE( &
656                        & zu(sobsgroups(jgroup)%sprofdataqc%nvprot(iuvar)), &
657                        & zv(sobsgroups(jgroup)%sprofdataqc%nvprot(ivvar))  &
658                        & )
659
660                     CALL obs_rotvel_pro( sobsgroups(jgroup)%sprofdataqc, sobsgroups(jgroup)%n2dint, &
661                        &                 iuvar, ivvar, zu, zv )
662
663                     DO jo = 1, sobsgroups(jgroup)%sprofdataqc%nprof
664                        DO jk = sobsgroups(jgroup)%sprofdataqc%npvsta(jo,iuvar), sobsgroups(jgroup)%sprofdataqc%npvend(jo,iuvar)
665                           sobsgroups(jgroup)%sprofdataqc%var(iuvar)%vmod(jk) = zu(jk)
666                        END DO
667                        DO jk = sobsgroups(jgroup)%sprofdataqc%npvsta(jo,ivvar), sobsgroups(jgroup)%sprofdataqc%npvend(jo,ivvar)
668                           sobsgroups(jgroup)%sprofdataqc%var(ivvar)%vmod(jk) = zv(jk)
669                        END DO
670                     END DO
671
672                     DEALLOCATE( zu )
673                     DEALLOCATE( zv )
674
675                  ELSE
676                     CALL ctl_stop( 'Could not identify velocity observation variables to rotate' )
677                  END IF
678               END IF
679
680               CALL obs_prof_decompress( sobsgroups(jgroup)%sprofdataqc, &
681                  &                      sobsgroups(jgroup)%sprofdata, .TRUE., numout )
682
683               ! Put additional and extra variable information into obswriinfo structure
684               ! used by obs_write.
685               ! add/ext variables generated by the OBS code (1...sobsgroups(jgroup)%naddvars)
686               ! may duplicate ones read in (%naddvars+1...sobsgroups(jgroup)%sprofdata%nadd)
687               ! Check for this, and if so only write out the version generated by the OBS code
688               sladd%inum = sobsgroups(jgroup)%sprofdata%nadd
689               ALLOCATE( ll_write(sobsgroups(jgroup)%sprofdata%nadd) )
690               ll_write(:) = .TRUE.
691               IF ( (sobsgroups(jgroup)%naddvars > 0) .AND. &
692                  & (sobsgroups(jgroup)%sprofdata%nadd > sobsgroups(jgroup)%naddvars) ) THEN
693                  DO jadd = sobsgroups(jgroup)%naddvars + 1, sobsgroups(jgroup)%sprofdata%nadd
694                     DO jadd2 = 1, sobsgroups(jgroup)%naddvars
695                        IF ( TRIM(sobsgroups(jgroup)%sprofdata%caddvars(jadd )) == &
696                           & TRIM(sobsgroups(jgroup)%sprofdata%caddvars(jadd2)) ) THEN
697                           sladd%inum = sladd%inum - 1
698                           ll_write(jadd) = .FALSE.
699                        ENDIF
700                     END DO
701                  END DO
702               ENDIF
703               IF ( sladd%inum > 0 ) THEN
704                  ALLOCATE( sladd%ipoint(sladd%inum),                                   &
705                     &      sladd%cdname(sladd%inum),                                   &
706                     &      sladd%cdlong(sladd%inum,sobsgroups(jgroup)%sprofdata%nvar), &
707                     &      sladd%cdunit(sladd%inum,sobsgroups(jgroup)%sprofdata%nvar) )
708                  jadd2 = 0
709                  DO jadd = 1, sobsgroups(jgroup)%sprofdata%nadd
710                     IF ( ll_write(jadd) ) THEN
711                        jadd2 = jadd2 + 1
712                        sladd%ipoint(jadd2) = jadd
713                        sladd%cdname(jadd2) = sobsgroups(jgroup)%sprofdata%caddvars(jadd)
714                        DO jvar = 1, sobsgroups(jgroup)%sprofdata%nvar
715                           sladd%cdlong(jadd2,jvar) = sobsgroups(jgroup)%sprofdata%caddlong(jadd,jvar)
716                           sladd%cdunit(jadd2,jvar) = sobsgroups(jgroup)%sprofdata%caddunit(jadd,jvar)
717                        END DO
718                     ENDIF
719                  END DO
720               ENDIF
721               DEALLOCATE( ll_write )
722               
723               slext%inum = sobsgroups(jgroup)%sprofdata%next
724               ALLOCATE( ll_write(sobsgroups(jgroup)%sprofdata%next) )
725               ll_write(:) = .TRUE.
726               IF ( (sobsgroups(jgroup)%nextvars > 0) .AND. &
727                  & (sobsgroups(jgroup)%sprofdata%next > sobsgroups(jgroup)%nextvars) ) THEN
728                  DO jext = sobsgroups(jgroup)%nextvars + 1, sobsgroups(jgroup)%sprofdata%next
729                     DO jext2 = 1, sobsgroups(jgroup)%nextvars
730                        IF ( TRIM(sobsgroups(jgroup)%sprofdata%cextvars(jext )) == &
731                           & TRIM(sobsgroups(jgroup)%sprofdata%cextvars(jext2)) ) THEN
732                           slext%inum = slext%inum - 1
733                           ll_write(jext) = .FALSE.
734                        ENDIF
735                     END DO
736                  END DO
737               ENDIF
738               IF ( slext%inum > 0 ) THEN
739                  ALLOCATE( slext%ipoint(slext%inum),   &
740                     &      slext%cdname(slext%inum),   &
741                     &      slext%cdlong(slext%inum,1), &
742                     &      slext%cdunit(slext%inum,1) )
743                  jext2 = 0
744                  DO jext = 1, sobsgroups(jgroup)%sprofdata%next
745                     IF ( ll_write(jext) ) THEN
746                        jext2 = jext2 + 1
747                        slext%ipoint(jext2)   = jext
748                        slext%cdname(jext2)   = sobsgroups(jgroup)%sprofdata%cextvars(jext)
749                        slext%cdlong(jext2,1) = sobsgroups(jgroup)%sprofdata%cextlong(jext)
750                        slext%cdunit(jext2,1) = sobsgroups(jgroup)%sprofdata%cextunit(jext)
751                     ENDIF
752                  END DO
753               ENDIF
754               DEALLOCATE( ll_write )
755
756               CALL obs_wri_prof( sobsgroups(jgroup)%sprofdata, sobsgroups(jgroup)%cgroupname, sladd, slext )
757
758               IF ( sladd%inum > 0 ) THEN
759                  DEALLOCATE( sladd%ipoint, sladd%cdname, sladd%cdlong, sladd%cdunit )
760               ENDIF
761               IF ( slext%inum > 0 ) THEN
762                  DEALLOCATE( slext%ipoint, slext%cdname, slext%cdlong, slext%cdunit )
763               ENDIF
764
765            ELSEIF (sobsgroups(jgroup)%lsurf) THEN
766
767               IF (sobsgroups(jgroup)%lvel) THEN
768                  iuvar = 0
769                  ivvar = 0
770                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
771                     IF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_uvel ) THEN
772                        iuvar = jvar
773                     ELSEIF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_vvel ) THEN
774                        ivvar = jvar
775                     ENDIF
776                  END DO
777                  IF ( (iuvar > 0) .AND. (ivvar > 0) ) THEN
778
779                     ! For velocity data, rotate the model velocities to N/S, E/W
780                     ! using the compressed data structure.
781                     ALLOCATE( &
782                        & zu(sobsgroups(jgroup)%ssurfdataqc%nsurf), &
783                        & zv(sobsgroups(jgroup)%ssurfdataqc%nsurf)  &
784                        & )
785
786                     CALL obs_rotvel_surf( sobsgroups(jgroup)%ssurfdataqc, sobsgroups(jgroup)%n2dint, &
787                        &                  iuvar, ivvar, zu, zv )
788
789                     DO jo = 1, sobsgroups(jgroup)%ssurfdataqc%nsurf
790                        sobsgroups(jgroup)%ssurfdataqc%rmod(jo,iuvar) = zu(jo)
791                        sobsgroups(jgroup)%ssurfdataqc%rmod(jo,ivvar) = zv(jo)
792                     END DO
793
794                     DEALLOCATE( zu )
795                     DEALLOCATE( zv )
796
797                  ELSE
798                     CALL ctl_stop( 'Could not identify velocity observation variables to rotate' )
799                  END IF
800               END IF
801
802               CALL obs_surf_decompress( sobsgroups(jgroup)%ssurfdataqc, &
803                  &                      sobsgroups(jgroup)%ssurfdata, .TRUE., numout )
804
805               IF (sobsgroups(jgroup)%lfbd) THEN
806                  ! Input observations were freeboard, but we're outputting ice thickness
807                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
808                     IF ( sobsgroups(jgroup)%cobstypes(jvar) == cobsname_fbd ) THEN
809                        sobsgroups(jgroup)%ssurfdata%cvars(jvar) = 'SIT'
810                        sobsgroups(jgroup)%ssurfdata%clong(jvar) = 'Sea ice thickness'
811                        sobsgroups(jgroup)%ssurfdata%cunit(jvar) = 'm'
812                        EXIT
813                     ENDIF
814                  END DO
815               ENDIF
816
817               ! Put additional and extra variable information into obswriinfo structure
818               ! used by obs_write.
819               ! add/ext variables generated by the OBS code (1...sobsgroups(jgroup)%naddvars)
820               ! may duplicate ones read in (%naddvars+1...sobsgroups(jgroup)%ssurfdata%nadd)
821               ! Check for this, and if so only write out the version generated by the OBS code
822               sladd%inum = sobsgroups(jgroup)%ssurfdata%nadd
823               ALLOCATE( ll_write(sobsgroups(jgroup)%ssurfdata%nadd) )
824               ll_write(:) = .TRUE.
825               IF ( (sobsgroups(jgroup)%naddvars > 0) .AND. &
826                  & (sobsgroups(jgroup)%ssurfdata%nadd > sobsgroups(jgroup)%naddvars) ) THEN
827                  DO jadd = sobsgroups(jgroup)%naddvars + 1, sobsgroups(jgroup)%ssurfdata%nadd
828                     DO jadd2 = 1, sobsgroups(jgroup)%naddvars
829                        IF ( TRIM(sobsgroups(jgroup)%ssurfdata%caddvars(jadd )) == &
830                           & TRIM(sobsgroups(jgroup)%ssurfdata%caddvars(jadd2)) ) THEN
831                           sladd%inum = sladd%inum - 1
832                           ll_write(jadd) = .FALSE.
833                        ENDIF
834                     END DO
835                  END DO
836               ENDIF
837               IF ( sladd%inum > 0 ) THEN
838                  ALLOCATE( sladd%ipoint(sladd%inum),                                   &
839                     &      sladd%cdname(sladd%inum),                                   &
840                     &      sladd%cdlong(sladd%inum,sobsgroups(jgroup)%ssurfdata%nvar), &
841                     &      sladd%cdunit(sladd%inum,sobsgroups(jgroup)%ssurfdata%nvar) )
842                  jadd2 = 0
843                  DO jadd = 1, sobsgroups(jgroup)%ssurfdata%nadd
844                     IF ( ll_write(jadd) ) THEN
845                        jadd2 = jadd2 + 1
846                        sladd%ipoint(jadd2) = jadd
847                        sladd%cdname(jadd2) = sobsgroups(jgroup)%ssurfdata%caddvars(jadd)
848                        DO jvar = 1, sobsgroups(jgroup)%ssurfdata%nvar
849                           sladd%cdlong(jadd2,jvar) = sobsgroups(jgroup)%ssurfdata%caddlong(jadd,jvar)
850                           sladd%cdunit(jadd2,jvar) = sobsgroups(jgroup)%ssurfdata%caddunit(jadd,jvar)
851                        END DO
852                     ENDIF
853                  END DO
854               ENDIF
855               DEALLOCATE( ll_write )
856               
857               slext%inum = sobsgroups(jgroup)%ssurfdata%nextra
858               ALLOCATE( ll_write(sobsgroups(jgroup)%ssurfdata%nextra) )
859               ll_write(:) = .TRUE.
860               IF ( (sobsgroups(jgroup)%nextvars > 0) .AND. &
861                  & (sobsgroups(jgroup)%ssurfdata%nextra > sobsgroups(jgroup)%nextvars) ) THEN
862                  DO jext = sobsgroups(jgroup)%nextvars + 1, sobsgroups(jgroup)%ssurfdata%nextra
863                     DO jext2 = 1, sobsgroups(jgroup)%nextvars
864                        IF ( TRIM(sobsgroups(jgroup)%ssurfdata%cextvars(jext )) == &
865                           & TRIM(sobsgroups(jgroup)%ssurfdata%cextvars(jext2)) ) THEN
866                           slext%inum = slext%inum - 1
867                           ll_write(jext) = .FALSE.
868                        ENDIF
869                     END DO
870                  END DO
871               ENDIF
872               IF ( slext%inum > 0 ) THEN
873                  ALLOCATE( slext%ipoint(slext%inum),   &
874                     &      slext%cdname(slext%inum),   &
875                     &      slext%cdlong(slext%inum,1), &
876                     &      slext%cdunit(slext%inum,1) )
877                  jext2 = 0
878                  DO jext = 1, sobsgroups(jgroup)%ssurfdata%nextra
879                     IF ( ll_write(jext) ) THEN
880                        jext2 = jext2 + 1
881                        slext%ipoint(jext2)   = jext
882                        slext%cdname(jext2)   = sobsgroups(jgroup)%ssurfdata%cextvars(jext)
883                        slext%cdlong(jext2,1) = sobsgroups(jgroup)%ssurfdata%cextlong(jext)
884                        slext%cdunit(jext2,1) = sobsgroups(jgroup)%ssurfdata%cextunit(jext)
885                     ENDIF
886                  END DO
887               ENDIF
888               DEALLOCATE( ll_write )
889
890               CALL obs_wri_surf( sobsgroups(jgroup)%ssurfdata, sobsgroups(jgroup)%cgroupname, sladd, slext )
891
892               IF ( sladd%inum > 0 ) THEN
893                  DEALLOCATE( sladd%ipoint, sladd%cdname, sladd%cdlong, sladd%cdunit )
894               ENDIF
895               IF ( slext%inum > 0 ) THEN
896                  DEALLOCATE( slext%ipoint, slext%cdname, slext%cdlong, slext%cdunit )
897               ENDIF
898
899            ENDIF
900
901         ENDIF
902
903      END DO
904
905      IF( ln_timing )   CALL timing_stop('dia_obs_wri')
906
907   END SUBROUTINE dia_obs_wri
908
909   SUBROUTINE calc_date( kstp, ddobs )
910      !!----------------------------------------------------------------------
911      !!                    ***  ROUTINE calc_date  ***
912      !!         
913      !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format
914      !!
915      !! ** Method  : Get date in double precision YYYYMMDD.HHMMSS format
916      !!
917      !! ** Action  : Get date in double precision YYYYMMDD.HHMMSS format
918      !!
919      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
920      !!
921      !! History :
922      !!        !  06-03  (K. Mogensen)  Original code
923      !!        !  06-05  (K. Mogensen)  Reformatted
924      !!        !  06-10  (A. Weaver) Cleaning
925      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
926      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
927      !!        !  2014-09  (D. Lea) New generic routine now deals with arbitrary initial time of day
928      !!----------------------------------------------------------------------
929      USE phycst, ONLY : &            ! Physical constants
930         & rday
931      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
932         & rdt
933
934      IMPLICIT NONE
935
936      !! * Arguments
937      REAL(KIND=dp), INTENT(OUT) :: ddobs                        ! Date in YYYYMMDD.HHMMSS
938      INTEGER :: kstp
939
940      !! * Local declarations
941      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
942      INTEGER :: imon
943      INTEGER :: iday
944      INTEGER :: ihou
945      INTEGER :: imin
946      INTEGER :: imday       ! Number of days in month.
947      REAL(wp) :: zdayfrc    ! Fraction of day
948
949      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
950
951      !!----------------------------------------------------------------------
952      !! Initial date initialization (year, month, day, hour, minute)
953      !!----------------------------------------------------------------------
954      iyea =   ndate0 / 10000
955      imon = ( ndate0 - iyea * 10000 ) / 100
956      iday =   ndate0 - iyea * 10000 - imon * 100
957      ihou =   nn_time0 / 100
958      imin = ( nn_time0 - ihou * 100 ) 
959
960      !!----------------------------------------------------------------------
961      !! Compute number of days + number of hours + min since initial time
962      !!----------------------------------------------------------------------
963      zdayfrc = kstp * rdt / rday
964      zdayfrc = zdayfrc - aint(zdayfrc)
965      imin = imin + int( zdayfrc * 24 * 60 ) 
966      DO WHILE (imin >= 60) 
967        imin=imin-60
968        ihou=ihou+1
969      END DO
970      DO WHILE (ihou >= 24)
971        ihou=ihou-24
972        iday=iday+1
973      END DO
974      iday = iday + kstp * rdt / rday 
975
976      !-----------------------------------------------------------------------
977      ! Convert number of days (iday) into a real date
978      !----------------------------------------------------------------------
979
980      CALL calc_month_len( iyea, imonth_len )
981
982      DO WHILE ( iday > imonth_len(imon) )
983         iday = iday - imonth_len(imon)
984         imon = imon + 1 
985         IF ( imon > 12 ) THEN
986            imon = 1
987            iyea = iyea + 1
988            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
989         ENDIF
990      END DO
991
992      !----------------------------------------------------------------------
993      ! Convert it into YYYYMMDD.HHMMSS format.
994      !----------------------------------------------------------------------
995      ddobs = iyea * 10000_dp + imon * 100_dp + &
996         &    iday + ihou * 0.01_dp + imin * 0.0001_dp
997
998   END SUBROUTINE calc_date
999
1000   SUBROUTINE ini_date( ddobsini )
1001      !!----------------------------------------------------------------------
1002      !!                    ***  ROUTINE ini_date  ***
1003      !!         
1004      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1005      !!
1006      !! ** Method  :
1007      !!
1008      !! ** Action  :
1009      !!
1010      !! History :
1011      !!        !  06-03  (K. Mogensen)  Original code
1012      !!        !  06-05  (K. Mogensen)  Reformatted
1013      !!        !  06-10  (A. Weaver) Cleaning
1014      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1015      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
1016      !!----------------------------------------------------------------------
1017
1018      IMPLICIT NONE
1019
1020      !! * Arguments
1021      REAL(KIND=dp), INTENT(OUT) :: ddobsini                   ! Initial date in YYYYMMDD.HHMMSS
1022
1023      CALL calc_date( nit000 - 1, ddobsini )
1024
1025   END SUBROUTINE ini_date
1026
1027   SUBROUTINE fin_date( ddobsfin )
1028      !!----------------------------------------------------------------------
1029      !!                    ***  ROUTINE fin_date  ***
1030      !!         
1031      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1032      !!
1033      !! ** Method  :
1034      !!
1035      !! ** Action  :
1036      !!
1037      !! History :
1038      !!        !  06-03  (K. Mogensen)  Original code
1039      !!        !  06-05  (K. Mogensen)  Reformatted
1040      !!        !  06-10  (A. Weaver) Cleaning
1041      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1042      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
1043      !!----------------------------------------------------------------------
1044
1045      IMPLICIT NONE
1046
1047      !! * Arguments
1048      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1049
1050      CALL calc_date( nitend, ddobsfin )
1051
1052   END SUBROUTINE fin_date
1053
1054END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.