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

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

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