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 branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864

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