source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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