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

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

Some bug fixes.
File size: 50.0 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               CALL obs_pre_prof( sobsgroups(jgroup)%sprofdata,     &
246                  &               sobsgroups(jgroup)%sprofdataqc,   &
247                  &               llvar,                            &
248                  &               jpi, jpj, jpk,                    &
249                  &               sobsgroups(jgroup)%rmask,         &
250                  &               sobsgroups(jgroup)%rglam,         &
251                  &               sobsgroups(jgroup)%rgphi,         &
252                  &               sobsgroups(jgroup)%lnea,          &
253                  &               sobsgroups(jgroup)%lbound_reject, &
254                  &               kdailyavtypes = sobsgroups(jgroup)%nprofdavtypes )
255               !
256               DEALLOCATE( llvar )
257               !
258            ELSEIF (sobsgroups(jgroup)%lsurf) THEN
259               !
260               ! Read in surface obs types
261               !
262               CALL obs_rea_surf( sobsgroups(jgroup)%ssurfdata,         &
263                  &               sobsgroups(jgroup)%nobsfiles,         &
264                  &               sobsgroups(jgroup)%cobsfiles,         &
265                  &               sobsgroups(jgroup)%nobstypes,         &
266                  &               sobsgroups(jgroup)%naddvars,          &
267                  &               sobsgroups(jgroup)%nextvars,          &
268                  &               nitend-nit000+2,                      &
269                  &               rn_dobsini,                           &
270                  &               rn_dobsend,                           &
271                  &               sobsgroups(jgroup)%rtime_mean_period, &
272                  &               sobsgroups(jgroup)%ltime_mean_bkg,    &
273                  &               sobsgroups(jgroup)%lignmis,           &
274                  &               .FALSE.,                              &
275                  &               sobsgroups(jgroup)%lnight,            &
276                  &               sobsgroups(jgroup)%cobstypes )
277               !
278               IF( sobsgroups(jgroup)%lsla ) THEN
279                  sobsgroups(jgroup)%ssurfdata%cextvars(sobsgroups(jgroup)%next_mdt) = 'MDT'
280                  sobsgroups(jgroup)%ssurfdata%cextlong(sobsgroups(jgroup)%next_mdt) = 'Mean dynamic topography'
281                  sobsgroups(jgroup)%ssurfdata%cextunit(sobsgroups(jgroup)%next_mdt) = 'Metres'
282                  sobsgroups(jgroup)%ssurfdata%caddvars(sobsgroups(jgroup)%nadd_ssh) = 'SSH'
283                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
284                     sobsgroups(jgroup)%ssurfdata%caddlong(sobsgroups(jgroup)%nadd_ssh,jvar) = 'Model Sea surface height'
285                     sobsgroups(jgroup)%ssurfdata%caddunit(sobsgroups(jgroup)%nadd_ssh,jvar) = 'Metres'
286                  END DO
287               ENDIF
288               !
289               IF( sobsgroups(jgroup)%lfbd ) THEN
290                  sobsgroups(jgroup)%ssurfdata%cextvars(sobsgroups(jgroup)%next_snow) = 'SNOW'
291                  sobsgroups(jgroup)%ssurfdata%cextlong(sobsgroups(jgroup)%next_snow) = 'Snow thickness'
292                  sobsgroups(jgroup)%ssurfdata%cextunit(sobsgroups(jgroup)%next_snow) = 'Metres'
293                  sobsgroups(jgroup)%ssurfdata%caddvars(sobsgroups(jgroup)%nadd_fbd)  = 'FBD'
294                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
295                     sobsgroups(jgroup)%ssurfdata%caddlong(sobsgroups(jgroup)%nadd_fbd,jvar) = 'Freeboard'
296                     sobsgroups(jgroup)%ssurfdata%caddunit(sobsgroups(jgroup)%nadd_fbd,jvar) = 'Metres'
297                  END DO
298               ENDIF
299               !
300               IF( sobsgroups(jgroup)%loutput_clim ) THEN
301                  sobsgroups(jgroup)%ssurfdata%caddvars(sobsgroups(jgroup)%nadd_clm)  = 'CLM'
302                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
303                     sobsgroups(jgroup)%ssurfdata%radd(:,:,jvar) = fbrmdi
304                     sobsgroups(jgroup)%ssurfdata%caddlong(sobsgroups(jgroup)%nadd_clm,jvar) = 'Climatology'
305                     sobsgroups(jgroup)%ssurfdata%caddunit(sobsgroups(jgroup)%nadd_clm,jvar) = sobsgroups(jgroup)%ssurfdata%cunit(jvar)
306                  END DO
307               ENDIF
308               !
309               CALL obs_pre_surf( sobsgroups(jgroup)%ssurfdata,      &
310                  &               sobsgroups(jgroup)%ssurfdataqc,    &
311                  &               jpi, jpj,                          &
312                  &               sobsgroups(jgroup)%rmask(:,:,1,:), &
313                  &               sobsgroups(jgroup)%rglam,          &
314                  &               sobsgroups(jgroup)%rgphi,          &
315                  &               sobsgroups(jgroup)%lnea,           &
316                  &               sobsgroups(jgroup)%lbound_reject )
317               !
318               IF( sobsgroups(jgroup)%lsla ) THEN
319                  CALL obs_rea_mdt( sobsgroups(jgroup)%ssurfdataqc, &
320                     &              sobsgroups(jgroup)%n2dint,      &
321                     &              sobsgroups(jgroup)%next_mdt,    &
322                     &              sobsgroups(jgroup)%nmsshc,      &
323                     &              sobsgroups(jgroup)%rmdtcorr,    &
324                     &              sobsgroups(jgroup)%rmdtcutoff )
325                  IF( sobsgroups(jgroup)%laltbias ) THEN
326                     !CALL obs_rea_altbias( sobsgroups(jgroup)%ssurfdataqc, &
327                     !   &                  sobsgroups(jgroup)%n2dint,      &
328                     !   &                  sobsgroups(jgroup)%caltbiasfile )
329                     CALL obs_app_bias( sobsgroups(jgroup)%ssurfdataqc,   &
330                        &               sobsgroups(jgroup)%next_mdt,      & 
331                        &               sobsgroups(jgroup)%n2dint,        & 
332                        &               1,                                &
333                        &               sobsgroups(jgroup)%caltbiasfile,  &
334                        &               'altbias',                        &
335                        &               ld_extvar=.TRUE. ) 
336                  ENDIF
337               ENDIF
338               !
339#if defined key_si3
340               IF( sobsgroups(jgroup)%lfbd ) THEN
341                  CALL obs_rea_snowdepth( sobsgroups(jgroup)%ssurfdataqc, &
342                     &                    sobsgroups(jgroup)%n2dint,      &
343                     &                    sobsgroups(jgroup)%next_snow,   &
344                     &                    hm_s(:,:) )
345               ENDIF
346#elif defined key_cice
347               IF( sobsgroups(jgroup)%lfbd ) THEN
348                  CALL obs_rea_snowdepth( sobsgroups(jgroup)%ssurfdataqc, &
349                     &                    sobsgroups(jgroup)%n2dint,      &
350                     &                    sobsgroups(jgroup)%next_snow,   &
351                     &                    thick_s(:,:) )
352               ENDIF
353#endif
354               !
355               IF( sobsgroups(jgroup)%lobsbias ) THEN
356                  CALL obs_app_bias( sobsgroups(jgroup)%ssurfdataqc,   &
357                     &               sobsgroups(jgroup)%nbiasvar,      & 
358                     &               sobsgroups(jgroup)%n2dint,        & 
359                     &               sobsgroups(jgroup)%nobsbiasfiles, &
360                     &               sobsgroups(jgroup)%cobsbiasfiles, &
361                     &               sobsgroups(jgroup)%cbiasvarname ) 
362               ENDIF
363               !
364            ENDIF
365         ENDIF
366      END DO
367      !
368   END SUBROUTINE dia_obs_init
369
370
371   SUBROUTINE dia_obs( kstp )
372      !!----------------------------------------------------------------------
373      !!                    ***  ROUTINE dia_obs  ***
374      !!         
375      !! ** Purpose : Call the observation operators on each time step
376      !!
377      !! ** Method  : Call the observation operators on each time step to
378      !!              compute the model equivalent of the following data:
379      !!               - Profile data, currently T/S or U/V
380      !!               - Surface data, currently SST, SLA or sea-ice concentration.
381      !!
382      !! ** Action  :
383      !!----------------------------------------------------------------------
384      USE dom_oce, ONLY : gdept_n, gdept_1d   ! Ocean space and time domain variables
385      USE phycst , ONLY : rday                ! Physical constants
386      USE oce    , ONLY : tsn, un, vn, sshn   ! Ocean dynamics and tracers variables
387      USE phycst , ONLY : rday                ! Physical constants
388#if defined key_si3
389      USE ice    , ONLY : at_i, hm_i          ! SI3 Ice model variables
390#elif defined key_cice
391      USE sbc_oce, ONLY : fr_i, thick_i       ! CICE Ice model variables
392#endif
393      USE tradmp,  ONLY : tclim, sclim        ! T&S climatology
394      USE obs_fbm, ONLY : fbrmdi              ! Real missing data indicator
395
396      IMPLICIT NONE
397
398      !! * Arguments
399      INTEGER, INTENT(IN) :: kstp  ! Current timestep
400      !! * Local declarations
401      INTEGER :: idaystp           ! Number of timesteps per day
402      INTEGER :: imeanstp          ! Number of timesteps for time averaging
403      INTEGER :: jtype             ! Data loop variable
404      INTEGER :: jvar              ! Variable number
405      INTEGER :: jgroup
406      INTEGER :: ji, jj, jobs      ! Loop counters
407      LOGICAL :: lstp0             ! Flag special treatment on zeroth time step
408      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
409         & zprofvar, &             ! Model values for variables in a prof ob
410         & zprofclim               ! Climatology values for variables in a prof ob
411      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
412         & zsurfvar, &             ! Model values for variables in a surf ob
413         & zsurfclim               ! Climatology values for variables in a surf ob
414
415      !-----------------------------------------------------------------------
416
417      IF( ln_timing )   CALL timing_start('dia_obs')
418
419      IF(lwp) THEN
420         WRITE(numout,*)
421         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
422         WRITE(numout,*) '~~~~~~~'
423      ENDIF
424
425      idaystp = NINT( rday / rdt )
426
427      !-----------------------------------------------------------------------
428      ! Call the profile and surface observation operators
429      !-----------------------------------------------------------------------
430
431      ALLOCATE( zprofvar(jpi,jpj,jpk),  &
432         &      zprofclim(jpi,jpj,jpk), &
433         &      zsurfvar(jpi,jpj),      &
434         &      zsurfclim(jpi,jpj) )
435
436      DO jgroup = 1, nn_obsgroups
437         IF (sobsgroups(jgroup)%lenabled) THEN
438
439            IF (sobsgroups(jgroup)%lprof) THEN
440
441               zprofclim(:,:,:) = fbrmdi
442
443               DO jvar = 1, sobsgroups(jgroup)%nobstypes
444
445                  SELECT CASE ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) )
446                  CASE('POTM')
447                     zprofvar(:,:,:) = tsn(:,:,:,jp_tem)
448                     IF (sobsgroups(jgroup)%loutput_clim) THEN
449                        zprofclim(:,:,:) = tclim(:,:,:)
450                     ENDIF
451                  CASE('PSAL')
452                     zprofvar(:,:,:) = tsn(:,:,:,jp_sal)
453                     IF (sobsgroups(jgroup)%loutput_clim) THEN
454                        zprofclim(:,:,:) = sclim(:,:,:)
455                     ENDIF
456                  CASE('UVEL')
457                     zprofvar(:,:,:) = un(:,:,:)
458                  CASE('VVEL')
459                     zprofvar(:,:,:) = vn(:,:,:)
460                  CASE DEFAULT
461                     CALL ctl_stop( 'Unknown profile observation type '//TRIM(sobsgroups(jgroup)%cobstypes(jvar))//' in dia_obs' )
462                  END SELECT
463
464                  CALL obs_prof_opt( sobsgroups(jgroup)%sprofdataqc,       &
465                     &               kstp, jpi, jpj, jpk,                  &
466                     &               nit000, idaystp, jvar,                &
467                     &               zprofvar,                             &
468                     &               sobsgroups(jgroup)%loutput_clim,      &
469                     &               sobsgroups(jgroup)%nadd_clm,          &
470                     &               zprofclim,                            &
471                     &               gdept_n,                              &
472                     &               gdepw_n,                              & 
473                     &               sobsgroups(jgroup)%rmask(:,:,:,jvar), &
474                     &               sobsgroups(jgroup)%rglam(:,:,jvar),   &
475                     &               sobsgroups(jgroup)%rgphi(:,:,jvar),   &
476                     &               sobsgroups(jgroup)%n1dint,            &
477                     &               sobsgroups(jgroup)%n2dint,            &
478                     &               kdailyavtypes = sobsgroups(jgroup)%nprofdavtypes )
479
480               END DO
481
482            ELSEIF (sobsgroups(jgroup)%lsurf) THEN
483
484               zsurfclim(:,:) = fbrmdi
485
486               DO jvar = 1, sobsgroups(jgroup)%nobstypes
487
488                  lstp0 = .FALSE.
489                  SELECT CASE ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) )
490                  CASE('SST')
491                     zsurfvar(:,:) = tsn(:,:,1,jp_tem)
492                     IF (sobsgroups(jgroup)%loutput_clim) THEN
493                        zsurfclim(:,:) = tclim(:,:,1)
494                     ENDIF
495                  CASE('SLA')
496                     zsurfvar(:,:) = sshn(:,:)
497                  CASE('SSS')
498                     zsurfvar(:,:) = tsn(:,:,1,jp_sal)
499                     IF (sobsgroups(jgroup)%loutput_clim) THEN
500                        zsurfclim(:,:) = sclim(:,:,1)
501                     ENDIF
502                  CASE('UVEL')
503                     zsurfvar(:,:) = un(:,:,1)
504                  CASE('VVEL')
505                     zsurfvar(:,:) = vn(:,:,1)
506                  CASE('ICECONC')
507                     IF ( kstp == 0 ) THEN
508                        lstp0 = .TRUE.
509                     ELSE
510#if defined key_si3
511                        zsurfvar(:,:) = at_i(:,:)
512#elif defined key_cice
513                        zsurfvar(:,:) = fr_i(:,:)
514#else
515                        CALL ctl_stop( ' Trying to run sea-ice observation operator', &
516                           &           ' but no sea-ice model appears to have been defined' )
517#endif
518                     ENDIF
519                  CASE('SIT','FBD')
520                     IF ( kstp == 0 ) THEN
521                        lstp0 = .TRUE.
522                     ELSE
523#if defined key_si3
524                        zsurfvar(:,:) = hm_i(:,:)
525#elif defined key_cice
526                        zsurfvar(:,:) = thick_i(:,:)
527#else
528                        CALL ctl_stop( ' Trying to run sea-ice observation operator', &
529                           &           ' but no sea-ice model appears to have been defined' )
530#endif
531                     ENDIF
532                  END SELECT
533
534                  IF ( lstp0 ) THEN
535                     IF ( sobsgroups(jgroup)%ssurfdataqc%nsstpmpp(1) > 0 ) THEN
536                        DO jobs = sobsgroups(jgroup)%ssurfdataqc%nsurfup + 1, &
537                           &      sobsgroups(jgroup)%ssurfdataqc%nsurfup + sobsgroups(jgroup)%ssurfdataqc%nsstp(1)
538                           sobsgroups(jgroup)%ssurfdata%nqc(jobs) = IBSET(sobsgroups(jgroup)%ssurfdata%nqc(jobs),13)
539                        END DO
540                        IF ( lwp ) THEN
541                           CALL ctl_warn( TRIM(sobsgroups(jgroup)%cobstypes(jvar))// &
542                              &           ' not initialised on zeroth '           // &
543                              &           'time-step but some obs are valid then.' )
544                           WRITE(numout,*)sobsgroups(jgroup)%ssurfdataqc%nsstpmpp(1), &
545                              &           TRIM(sobsgroups(jgroup)%cobstypes(jvar)),   &
546                              &           'observations will be flagged as bad'
547                        ENDIF
548                     ENDIF
549                     IF ( jvar == sobsgroups(jgroup)%ssurfdataqc%nvar ) THEN
550                        sobsgroups(jgroup)%ssurfdataqc%nsurfup = sobsgroups(jgroup)%ssurfdataqc%nsurfup + &
551                           &                                     sobsgroups(jgroup)%ssurfdataqc%nsstp(1)
552                     ENDIF
553                  ELSE
554                     IF ( sobsgroups(jgroup)%ltime_mean_bkg ) THEN
555                        ! Number of time-steps in meaning period
556                        imeanstp = NINT( ( sobsgroups(jgroup)%rtime_mean_period * 60.0 * 60.0 ) / rdt )
557                     ENDIF
558                     CALL obs_surf_opt( sobsgroups(jgroup)%ssurfdataqc,       &
559                        &               kstp, jpi, jpj,                       &
560                        &               nit000, idaystp,                      &
561                        &               jvar, zsurfvar,                       &
562                        &               sobsgroups(jgroup)%loutput_clim,      &
563                        &               sobsgroups(jgroup)%nadd_clm,          &
564                        &               zsurfclim,                            &
565                        &               sobsgroups(jgroup)%rmask(:,:,1,jvar), &
566                        &               sobsgroups(jgroup)%n2dint,            &
567                        &               sobsgroups(jgroup)%lnight,            &
568                        &               sobsgroups(jgroup)%ravglamscl,        &
569                        &               sobsgroups(jgroup)%ravgphiscl,        &
570                        &               sobsgroups(jgroup)%lfp_indegs,        &
571                        &               sobsgroups(jgroup)%ltime_mean_bkg,    &
572                        &               imeanstp,                             &
573                        &               kssh=sobsgroups(jgroup)%nadd_ssh,     &
574                        &               kmdt=sobsgroups(jgroup)%next_mdt,     &
575                        &               kfbd=sobsgroups(jgroup)%nadd_fbd,     &
576                        &               ksnow=sobsgroups(jgroup)%next_snow )
577                  ENDIF
578
579               END DO
580
581            ENDIF
582
583         ENDIF
584      END DO
585
586      DEALLOCATE( zprofvar, zprofclim, &
587         &        zsurfvar, zsurfclim )
588
589      IF( ln_timing )   CALL timing_stop('dia_obs')
590
591   END SUBROUTINE dia_obs
592
593   SUBROUTINE dia_obs_wri
594      !!----------------------------------------------------------------------
595      !!                    ***  ROUTINE dia_obs_wri  ***
596      !!         
597      !! ** Purpose : Call observation diagnostic output routines
598      !!
599      !! ** Method  : Call observation diagnostic output routines
600      !!
601      !! ** Action  :
602      !!
603      !! History :
604      !!        !  06-03  (K. Mogensen) Original code
605      !!        !  06-05  (K. Mogensen) Reformatted
606      !!        !  06-10  (A. Weaver) Cleaning
607      !!        !  07-03  (K. Mogensen) General handling of profiles
608      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
609      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
610      !!----------------------------------------------------------------------
611      !! * Modules used
612      USE obs_rot_vel          ! Rotation of velocities
613
614      IMPLICIT NONE
615
616      !! * Local declarations
617      INTEGER :: jgroup                   ! Data set loop variable
618      INTEGER :: jo, jvar, jk, jadd, jext, jadd2, jext2
619      INTEGER :: iuvar, ivvar
620      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
621         & zu, &
622         & zv
623      LOGICAL, DIMENSION(:), ALLOCATABLE :: ll_write
624      TYPE(obswriinfo) :: sladd, slext
625
626      IF( ln_timing )   CALL timing_start('dia_obs_wri')
627
628      !-----------------------------------------------------------------------
629      ! Depending on switches call various observation output routines
630      !-----------------------------------------------------------------------
631
632      DO jgroup = 1, nn_obsgroups
633         IF (sobsgroups(jgroup)%lenabled) THEN
634
635            IF (sobsgroups(jgroup)%lprof) THEN
636
637               IF (sobsgroups(jgroup)%lvel) THEN
638                  iuvar = 0
639                  ivvar = 0
640                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
641                     IF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_uvel ) THEN
642                        iuvar = jvar
643                     ELSEIF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_vvel ) THEN
644                        ivvar = jvar
645                     ENDIF
646                  END DO
647                  IF ( (iuvar > 0) .AND. (ivvar > 0) ) THEN
648
649                     ! For velocity data, rotate the model velocities to N/S, E/W
650                     ! using the compressed data structure.
651                     ALLOCATE( &
652                        & zu(sobsgroups(jgroup)%sprofdataqc%nvprot(iuvar)), &
653                        & zv(sobsgroups(jgroup)%sprofdataqc%nvprot(ivvar))  &
654                        & )
655
656                     CALL obs_rotvel_pro( sobsgroups(jgroup)%sprofdataqc, sobsgroups(jgroup)%n2dint, &
657                        &                 iuvar, ivvar, zu, zv )
658
659                     DO jo = 1, sobsgroups(jgroup)%sprofdataqc%nprof
660                        DO jk = sobsgroups(jgroup)%sprofdataqc%npvsta(jo,iuvar), sobsgroups(jgroup)%sprofdataqc%npvend(jo,iuvar)
661                           sobsgroups(jgroup)%sprofdataqc%var(iuvar)%vmod(jk) = zu(jk)
662                        END DO
663                        DO jk = sobsgroups(jgroup)%sprofdataqc%npvsta(jo,ivvar), sobsgroups(jgroup)%sprofdataqc%npvend(jo,ivvar)
664                           sobsgroups(jgroup)%sprofdataqc%var(ivvar)%vmod(jk) = zv(jk)
665                        END DO
666                     END DO
667
668                     DEALLOCATE( zu )
669                     DEALLOCATE( zv )
670
671                  ELSE
672                     CALL ctl_stop( 'Could not identify velocity observation variables to rotate' )
673                  END IF
674               END IF
675
676               CALL obs_prof_decompress( sobsgroups(jgroup)%sprofdataqc, &
677                  &                      sobsgroups(jgroup)%sprofdata, .TRUE., numout )
678
679               ! Put additional and extra variable information into obswriinfo structure
680               ! used by obs_write.
681               ! add/ext variables generated by the OBS code (1...sobsgroups(jgroup)%naddvars)
682               ! may duplicate ones read in (%naddvars+1...sobsgroups(jgroup)%sprofdata%nadd)
683               ! Check for this, and if so only write out the version generated by the OBS code
684               sladd%inum = sobsgroups(jgroup)%sprofdata%nadd
685               ALLOCATE( ll_write(sobsgroups(jgroup)%sprofdata%nadd) )
686               ll_write(:) = .TRUE.
687               IF ( (sobsgroups(jgroup)%naddvars > 0) .AND. &
688                  & (sobsgroups(jgroup)%sprofdata%nadd > sobsgroups(jgroup)%naddvars) ) THEN
689                  DO jadd = sobsgroups(jgroup)%naddvars + 1, sobsgroups(jgroup)%sprofdata%nadd
690                     DO jadd2 = 1, sobsgroups(jgroup)%naddvars
691                        IF ( TRIM(sobsgroups(jgroup)%sprofdata%caddvars(jadd )) == &
692                           & TRIM(sobsgroups(jgroup)%sprofdata%caddvars(jadd2)) ) THEN
693                           sladd%inum = sladd%inum - 1
694                           ll_write(jadd) = .FALSE.
695                        ENDIF
696                     END DO
697                  END DO
698               ENDIF
699               IF ( sladd%inum > 0 ) THEN
700                  ALLOCATE( sladd%ipoint(sladd%inum),                                   &
701                     &      sladd%cdname(sladd%inum),                                   &
702                     &      sladd%cdlong(sladd%inum,sobsgroups(jgroup)%sprofdata%nvar), &
703                     &      sladd%cdunit(sladd%inum,sobsgroups(jgroup)%sprofdata%nvar) )
704                  jadd2 = 0
705                  DO jadd = 1, sobsgroups(jgroup)%sprofdata%nadd
706                     IF ( ll_write(jadd) ) THEN
707                        jadd2 = jadd2 + 1
708                        sladd%ipoint(jadd2) = jadd
709                        sladd%cdname(jadd2) = sobsgroups(jgroup)%sprofdata%caddvars(jadd)
710                        DO jvar = 1, sobsgroups(jgroup)%sprofdata%nvar
711                           sladd%cdlong(jadd2,jvar) = sobsgroups(jgroup)%sprofdata%caddlong(jadd,jvar)
712                           sladd%cdunit(jadd2,jvar) = sobsgroups(jgroup)%sprofdata%caddunit(jadd,jvar)
713                        END DO
714                     ENDIF
715                  END DO
716               ENDIF
717               DEALLOCATE( ll_write )
718               
719               slext%inum = sobsgroups(jgroup)%sprofdata%next
720               ALLOCATE( ll_write(sobsgroups(jgroup)%sprofdata%next) )
721               ll_write(:) = .TRUE.
722               IF ( (sobsgroups(jgroup)%nextvars > 0) .AND. &
723                  & (sobsgroups(jgroup)%sprofdata%next > sobsgroups(jgroup)%nextvars) ) THEN
724                  DO jext = sobsgroups(jgroup)%nextvars + 1, sobsgroups(jgroup)%sprofdata%next
725                     DO jext2 = 1, sobsgroups(jgroup)%nextvars
726                        IF ( TRIM(sobsgroups(jgroup)%sprofdata%cextvars(jext )) == &
727                           & TRIM(sobsgroups(jgroup)%sprofdata%cextvars(jext2)) ) THEN
728                           slext%inum = slext%inum - 1
729                           ll_write(jext) = .FALSE.
730                        ENDIF
731                     END DO
732                  END DO
733               ENDIF
734               IF ( slext%inum > 0 ) THEN
735                  ALLOCATE( slext%ipoint(slext%inum),   &
736                     &      slext%cdname(slext%inum),   &
737                     &      slext%cdlong(slext%inum,1), &
738                     &      slext%cdunit(slext%inum,1) )
739                  jext2 = 0
740                  DO jext = 1, sobsgroups(jgroup)%sprofdata%next
741                     IF ( ll_write(jext) ) THEN
742                        jext2 = jext2 + 1
743                        slext%ipoint(jext2)   = jext
744                        slext%cdname(jext2)   = sobsgroups(jgroup)%sprofdata%cextvars(jext)
745                        slext%cdlong(jext2,1) = sobsgroups(jgroup)%sprofdata%cextlong(jext)
746                        slext%cdunit(jext2,1) = sobsgroups(jgroup)%sprofdata%cextunit(jext)
747                     ENDIF
748                  END DO
749               ENDIF
750               DEALLOCATE( ll_write )
751
752               CALL obs_wri_prof( sobsgroups(jgroup)%sprofdata, sobsgroups(jgroup)%cgroupname, sladd, slext )
753
754               IF ( sladd%inum > 0 ) THEN
755                  DEALLOCATE( sladd%ipoint, sladd%cdname, sladd%cdlong, sladd%cdunit )
756               ENDIF
757               IF ( slext%inum > 0 ) THEN
758                  DEALLOCATE( slext%ipoint, slext%cdname, slext%cdlong, slext%cdunit )
759               ENDIF
760
761            ELSEIF (sobsgroups(jgroup)%lsurf) THEN
762
763               IF (sobsgroups(jgroup)%lvel) THEN
764                  iuvar = 0
765                  ivvar = 0
766                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
767                     IF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_uvel ) THEN
768                        iuvar = jvar
769                     ELSEIF ( TRIM(sobsgroups(jgroup)%cobstypes(jvar)) == cobsname_vvel ) THEN
770                        ivvar = jvar
771                     ENDIF
772                  END DO
773                  IF ( (iuvar > 0) .AND. (ivvar > 0) ) THEN
774
775                     ! For velocity data, rotate the model velocities to N/S, E/W
776                     ! using the compressed data structure.
777                     ALLOCATE( &
778                        & zu(sobsgroups(jgroup)%ssurfdataqc%nsurf), &
779                        & zv(sobsgroups(jgroup)%ssurfdataqc%nsurf)  &
780                        & )
781
782                     CALL obs_rotvel_surf( sobsgroups(jgroup)%ssurfdataqc, sobsgroups(jgroup)%n2dint, &
783                        &                  iuvar, ivvar, zu, zv )
784
785                     DO jo = 1, sobsgroups(jgroup)%ssurfdataqc%nsurf
786                        sobsgroups(jgroup)%ssurfdataqc%rmod(jo,iuvar) = zu(jo)
787                        sobsgroups(jgroup)%ssurfdataqc%rmod(jo,ivvar) = zv(jo)
788                     END DO
789
790                     DEALLOCATE( zu )
791                     DEALLOCATE( zv )
792
793                  ELSE
794                     CALL ctl_stop( 'Could not identify velocity observation variables to rotate' )
795                  END IF
796               END IF
797
798               CALL obs_surf_decompress( sobsgroups(jgroup)%ssurfdataqc, &
799                  &                      sobsgroups(jgroup)%ssurfdata, .TRUE., numout )
800
801               IF (sobsgroups(jgroup)%lfbd) THEN
802                  ! Input observations were freeboard, but we're outputting ice thickness
803                  DO jvar = 1, sobsgroups(jgroup)%nobstypes
804                     IF ( sobsgroups(jgroup)%cobstypes(jvar) == cobsname_fbd ) THEN
805                        sobsgroups(jgroup)%ssurfdata%cvars(jvar) = 'SIT'
806                        sobsgroups(jgroup)%ssurfdata%clong(jvar) = 'Sea ice thickness'
807                        sobsgroups(jgroup)%ssurfdata%cunit(jvar) = 'm'
808                        EXIT
809                     ENDIF
810                  END DO
811               ENDIF
812
813               ! Put additional and extra variable information into obswriinfo structure
814               ! used by obs_write.
815               ! add/ext variables generated by the OBS code (1...sobsgroups(jgroup)%naddvars)
816               ! may duplicate ones read in (%naddvars+1...sobsgroups(jgroup)%ssurfdata%nadd)
817               ! Check for this, and if so only write out the version generated by the OBS code
818               sladd%inum = sobsgroups(jgroup)%ssurfdata%nadd
819               ALLOCATE( ll_write(sobsgroups(jgroup)%ssurfdata%nadd) )
820               ll_write(:) = .TRUE.
821               IF ( (sobsgroups(jgroup)%naddvars > 0) .AND. &
822                  & (sobsgroups(jgroup)%ssurfdata%nadd > sobsgroups(jgroup)%naddvars) ) THEN
823                  DO jadd = sobsgroups(jgroup)%naddvars + 1, sobsgroups(jgroup)%ssurfdata%nadd
824                     DO jadd2 = 1, sobsgroups(jgroup)%naddvars
825                        IF ( TRIM(sobsgroups(jgroup)%ssurfdata%caddvars(jadd )) == &
826                           & TRIM(sobsgroups(jgroup)%ssurfdata%caddvars(jadd2)) ) THEN
827                           sladd%inum = sladd%inum - 1
828                           ll_write(jadd) = .FALSE.
829                        ENDIF
830                     END DO
831                  END DO
832               ENDIF
833               IF ( sladd%inum > 0 ) THEN
834                  ALLOCATE( sladd%ipoint(sladd%inum),                                   &
835                     &      sladd%cdname(sladd%inum),                                   &
836                     &      sladd%cdlong(sladd%inum,sobsgroups(jgroup)%ssurfdata%nvar), &
837                     &      sladd%cdunit(sladd%inum,sobsgroups(jgroup)%ssurfdata%nvar) )
838                  jadd2 = 0
839                  DO jadd = 1, sobsgroups(jgroup)%ssurfdata%nadd
840                     IF ( ll_write(jadd) ) THEN
841                        jadd2 = jadd2 + 1
842                        sladd%ipoint(jadd2) = jadd
843                        sladd%cdname(jadd2) = sobsgroups(jgroup)%ssurfdata%caddvars(jadd)
844                        DO jvar = 1, sobsgroups(jgroup)%ssurfdata%nvar
845                           sladd%cdlong(jadd2,jvar) = sobsgroups(jgroup)%ssurfdata%caddlong(jadd,jvar)
846                           sladd%cdunit(jadd2,jvar) = sobsgroups(jgroup)%ssurfdata%caddunit(jadd,jvar)
847                        END DO
848                     ENDIF
849                  END DO
850               ENDIF
851               DEALLOCATE( ll_write )
852               
853               slext%inum = sobsgroups(jgroup)%ssurfdata%nextra
854               ALLOCATE( ll_write(sobsgroups(jgroup)%ssurfdata%nextra) )
855               ll_write(:) = .TRUE.
856               IF ( (sobsgroups(jgroup)%nextvars > 0) .AND. &
857                  & (sobsgroups(jgroup)%ssurfdata%nextra > sobsgroups(jgroup)%nextvars) ) THEN
858                  DO jext = sobsgroups(jgroup)%nextvars + 1, sobsgroups(jgroup)%ssurfdata%nextra
859                     DO jext2 = 1, sobsgroups(jgroup)%nextvars
860                        IF ( TRIM(sobsgroups(jgroup)%ssurfdata%cextvars(jext )) == &
861                           & TRIM(sobsgroups(jgroup)%ssurfdata%cextvars(jext2)) ) THEN
862                           slext%inum = slext%inum - 1
863                           ll_write(jext) = .FALSE.
864                        ENDIF
865                     END DO
866                  END DO
867               ENDIF
868               IF ( slext%inum > 0 ) THEN
869                  ALLOCATE( slext%ipoint(slext%inum),   &
870                     &      slext%cdname(slext%inum),   &
871                     &      slext%cdlong(slext%inum,1), &
872                     &      slext%cdunit(slext%inum,1) )
873                  jext2 = 0
874                  DO jext = 1, sobsgroups(jgroup)%ssurfdata%nextra
875                     IF ( ll_write(jext) ) THEN
876                        jext2 = jext2 + 1
877                        slext%ipoint(jext2)   = jext
878                        slext%cdname(jext2)   = sobsgroups(jgroup)%ssurfdata%cextvars(jext)
879                        slext%cdlong(jext2,1) = sobsgroups(jgroup)%ssurfdata%cextlong(jext)
880                        slext%cdunit(jext2,1) = sobsgroups(jgroup)%ssurfdata%cextunit(jext)
881                     ENDIF
882                  END DO
883               ENDIF
884               DEALLOCATE( ll_write )
885
886               CALL obs_wri_surf( sobsgroups(jgroup)%ssurfdata, sobsgroups(jgroup)%cgroupname, sladd, slext )
887
888               IF ( sladd%inum > 0 ) THEN
889                  DEALLOCATE( sladd%ipoint, sladd%cdname, sladd%cdlong, sladd%cdunit )
890               ENDIF
891               IF ( slext%inum > 0 ) THEN
892                  DEALLOCATE( slext%ipoint, slext%cdname, slext%cdlong, slext%cdunit )
893               ENDIF
894
895            ENDIF
896
897         ENDIF
898
899      END DO
900
901      IF( ln_timing )   CALL timing_stop('dia_obs_wri')
902
903   END SUBROUTINE dia_obs_wri
904
905   SUBROUTINE calc_date( kstp, ddobs )
906      !!----------------------------------------------------------------------
907      !!                    ***  ROUTINE calc_date  ***
908      !!         
909      !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format
910      !!
911      !! ** Method  : Get date in double precision YYYYMMDD.HHMMSS format
912      !!
913      !! ** Action  : Get date in double precision YYYYMMDD.HHMMSS format
914      !!
915      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
916      !!
917      !! History :
918      !!        !  06-03  (K. Mogensen)  Original code
919      !!        !  06-05  (K. Mogensen)  Reformatted
920      !!        !  06-10  (A. Weaver) Cleaning
921      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
922      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
923      !!        !  2014-09  (D. Lea) New generic routine now deals with arbitrary initial time of day
924      !!----------------------------------------------------------------------
925      USE phycst, ONLY : &            ! Physical constants
926         & rday
927      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
928         & rdt
929
930      IMPLICIT NONE
931
932      !! * Arguments
933      REAL(KIND=dp), INTENT(OUT) :: ddobs                        ! Date in YYYYMMDD.HHMMSS
934      INTEGER :: kstp
935
936      !! * Local declarations
937      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
938      INTEGER :: imon
939      INTEGER :: iday
940      INTEGER :: ihou
941      INTEGER :: imin
942      INTEGER :: imday       ! Number of days in month.
943      REAL(wp) :: zdayfrc    ! Fraction of day
944
945      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
946
947      !!----------------------------------------------------------------------
948      !! Initial date initialization (year, month, day, hour, minute)
949      !!----------------------------------------------------------------------
950      iyea =   ndate0 / 10000
951      imon = ( ndate0 - iyea * 10000 ) / 100
952      iday =   ndate0 - iyea * 10000 - imon * 100
953      ihou =   nn_time0 / 100
954      imin = ( nn_time0 - ihou * 100 ) 
955
956      !!----------------------------------------------------------------------
957      !! Compute number of days + number of hours + min since initial time
958      !!----------------------------------------------------------------------
959      zdayfrc = kstp * rdt / rday
960      zdayfrc = zdayfrc - aint(zdayfrc)
961      imin = imin + int( zdayfrc * 24 * 60 ) 
962      DO WHILE (imin >= 60) 
963        imin=imin-60
964        ihou=ihou+1
965      END DO
966      DO WHILE (ihou >= 24)
967        ihou=ihou-24
968        iday=iday+1
969      END DO
970      iday = iday + kstp * rdt / rday 
971
972      !-----------------------------------------------------------------------
973      ! Convert number of days (iday) into a real date
974      !----------------------------------------------------------------------
975
976      CALL calc_month_len( iyea, imonth_len )
977
978      DO WHILE ( iday > imonth_len(imon) )
979         iday = iday - imonth_len(imon)
980         imon = imon + 1 
981         IF ( imon > 12 ) THEN
982            imon = 1
983            iyea = iyea + 1
984            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
985         ENDIF
986      END DO
987
988      !----------------------------------------------------------------------
989      ! Convert it into YYYYMMDD.HHMMSS format.
990      !----------------------------------------------------------------------
991      ddobs = iyea * 10000_dp + imon * 100_dp + &
992         &    iday + ihou * 0.01_dp + imin * 0.0001_dp
993
994   END SUBROUTINE calc_date
995
996   SUBROUTINE ini_date( ddobsini )
997      !!----------------------------------------------------------------------
998      !!                    ***  ROUTINE ini_date  ***
999      !!         
1000      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1001      !!
1002      !! ** Method  :
1003      !!
1004      !! ** Action  :
1005      !!
1006      !! History :
1007      !!        !  06-03  (K. Mogensen)  Original code
1008      !!        !  06-05  (K. Mogensen)  Reformatted
1009      !!        !  06-10  (A. Weaver) Cleaning
1010      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1011      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
1012      !!----------------------------------------------------------------------
1013
1014      IMPLICIT NONE
1015
1016      !! * Arguments
1017      REAL(KIND=dp), INTENT(OUT) :: ddobsini                   ! Initial date in YYYYMMDD.HHMMSS
1018
1019      CALL calc_date( nit000 - 1, ddobsini )
1020
1021   END SUBROUTINE ini_date
1022
1023   SUBROUTINE fin_date( ddobsfin )
1024      !!----------------------------------------------------------------------
1025      !!                    ***  ROUTINE fin_date  ***
1026      !!         
1027      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1028      !!
1029      !! ** Method  :
1030      !!
1031      !! ** Action  :
1032      !!
1033      !! History :
1034      !!        !  06-03  (K. Mogensen)  Original code
1035      !!        !  06-05  (K. Mogensen)  Reformatted
1036      !!        !  06-10  (A. Weaver) Cleaning
1037      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1038      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
1039      !!----------------------------------------------------------------------
1040
1041      IMPLICIT NONE
1042
1043      !! * Arguments
1044      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1045
1046      CALL calc_date( nitend, ddobsfin )
1047
1048   END SUBROUTINE fin_date
1049
1050END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.