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/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6044

Last change on this file since 6044 was 6044, checked in by timgraham, 9 years ago

Merged branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER into branch

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