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/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8733

Last change on this file since 8733 was 8733, checked in by dancopsey, 6 years ago

Remove svn keywords.

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