New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaobs.F90 in branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

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

Last change on this file since 7911 was 7911, checked in by timgraham, 7 years ago

Minor bug fixes with pointers in bdy and also remove wrk_alloc calls in test cases

  • Property svn:keywords set to Id
File size: 39.2 KB
RevLine 
[2128]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
[6140]12   !!   calc_date    : Compute the date of timestep in YYYYMMDD.HHMMSS format
[2128]13   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS
14   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
15   !!----------------------------------------------------------------------
[6140]16   !! * Modules used
[2128]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
[6140]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
[2128]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
[2715]34   USE lib_mpp                  ! For ctl_warn/stop
[2128]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
[4245]42      &   dia_obs_wri,  &  ! Write model equivalent to observations
[6140]43      &   dia_obs_dealloc, &  ! Deallocate dia_obs data
44      &   calc_date           ! Compute the date of a timestep
[2128]45
46   !! * Module variables
[6140]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
[2128]52   INTEGER, DIMENSION(imaxavtypes) :: &
[6140]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
[2128]69
[6140]70   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
71      & cobstypesprof, &     !: Profile obs types
72      & cobstypessurf        !: Surface obs types
[2128]73
[2287]74   !!----------------------------------------------------------------------
75   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
76   !! $Id$
77   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
78   !!----------------------------------------------------------------------
79
[2128]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
[6140]97      !!        !  14-08  (J.While) Incorporated SST bias correction 
98      !!        !  15-02  (M. Martin) Simplification of namelist and code
[2128]99      !!----------------------------------------------------------------------
100
101      IMPLICIT NONE
102
103      !! * Local declarations
[6140]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
[7910]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
[6140]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,             &
[2128]155         &            ln_grid_search_lookup,                          &
[6140]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
[2128]164
[6140]165      INTEGER :: jnumsstbias
[2128]166
167      !-----------------------------------------------------------------------
168      ! Read namelist parameters
169      !-----------------------------------------------------------------------
[6140]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
[2128]181
[6140]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
[4147]187      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
188901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
[2128]189
[6140]190      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
[4147]191      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
192902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
[4624]193      IF(lwm) WRITE ( numond, namobs )
[4147]194
[6140]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'
[2128]198         RETURN
199      ENDIF
200     
201      !-----------------------------------------------------------------------
[6140]202      ! Set up list of observation types to be used
203      ! and the files associated with each type
[2128]204      !-----------------------------------------------------------------------
205
[6140]206      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
207      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) )
[2128]208
[6140]209      IF (ln_sstbias) THEN
210         lmask(:) = .FALSE. 
211         WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 
212         jnumsstbias = COUNT(lmask) 
213         lmask(:) = .FALSE. 
214      ENDIF     
[2128]215
[6140]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
[2128]225
[6140]226      IF ( nproftypes > 0 ) THEN
[2128]227
[6140]228         ALLOCATE( cobstypesprof(nproftypes) )
229         ALLOCATE( ifilesprof(nproftypes) )
230         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
[2128]231
[6140]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
[2128]241            END DO
242         ENDIF
[6140]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
[2128]253
[6140]254      ENDIF
[2128]255
[6140]256      IF ( nsurftypes > 0 ) THEN
[2128]257
[6140]258         ALLOCATE( cobstypessurf(nsurftypes) )
259         ALLOCATE( ifilessurf(nsurftypes) )
260         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
[2128]261
[6140]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
[2128]271            END DO
272         ENDIF
[6140]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
[2128]281            END DO
282         ENDIF
[6140]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
[2128]295
296      ENDIF
297
[6140]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
[2128]328
[6140]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
[2128]336         ENDIF
337
[6140]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
[2128]346         ENDIF
[6140]347         WRITE(numout,*) '~~~~~~~~~~~~'
[2128]348
[6140]349      ENDIF
[2128]350
[6140]351      !-----------------------------------------------------------------------
352      ! Obs operator parameter checking and initialisations
353      !-----------------------------------------------------------------------
[2128]354
[6140]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
[2128]359
[6140]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' )
[2128]362      ENDIF
363
[6140]364      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
365         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
366            &                    ' is not available')
[2128]367      ENDIF
368
[6140]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
[2128]373
[6140]374      CALL obs_typ_init
375      IF(ln_grid_global) THEN
376         CALL mppmap_init
377      ENDIF
[2128]378
[6140]379      CALL obs_grid_setup( )
[2128]380
[6140]381      !-----------------------------------------------------------------------
382      ! Depending on switches read the various observation types
383      !-----------------------------------------------------------------------
[2128]384
[6140]385      IF ( nproftypes > 0 ) THEN
[2128]386
[6140]387         ALLOCATE(profdata(nproftypes))
388         ALLOCATE(profdataqc(nproftypes))
389         ALLOCATE(nvarsprof(nproftypes))
390         ALLOCATE(nextrprof(nproftypes))
[2128]391
[6140]392         DO jtype = 1, nproftypes
[2128]393
[6140]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
[2128]417
[6140]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 )
[2128]425
[6140]426            DO jvar = 1, nvarsprof(jtype)
427               CALL obs_prof_staend( profdata(jtype), jvar )
428            END DO
[3651]429
[6140]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 )
[2128]435
[6140]436         END DO
[3651]437
[6140]438         DEALLOCATE( ifilesprof, clproffiles )
[2128]439
440      ENDIF
441
[6140]442      IF ( nsurftypes > 0 ) THEN
[2128]443
[6140]444         ALLOCATE(surfdata(nsurftypes))
445         ALLOCATE(surfdataqc(nsurftypes))
446         ALLOCATE(nvarssurf(nsurftypes))
447         ALLOCATE(nextrsurf(nsurftypes))
[2128]448
[6140]449         DO jtype = 1, nsurftypes
[2128]450
[6140]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
[2128]456
[6140]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 )
[2128]462         
463         
[6140]464            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea )
[2128]465
[6140]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
[2128]479
[6140]480         DEALLOCATE( ifilessurf, clsurffiles )
[2128]481
[6140]482      ENDIF
[2128]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
[6140]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.
[2128]497      !!
[6140]498      !! ** Action  :
[2128]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
[6140]507      !!        !  14-08  (J. While) observation operator for profiles in
508      !!                             generalised vertical coordinates
509      !!        !  15-08  (M. Martin) Combined surface/profile routines.
[2128]510      !!----------------------------------------------------------------------
511      !! * Modules used
512      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
[6140]513         & gdept_n,       &     
514         & gdept_1d     
[2128]515      USE phycst, ONLY : &              ! Physical constants
516         & rday                         
517      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
[3294]518         & tsn,  &             
[6140]519         & un, vn, &
520         & sshn 
521      USE phycst, ONLY : &         ! Physical constants
522         & rday
[3294]523#if defined  key_lim3
[6140]524      USE ice, ONLY : &            ! LIM3 Ice model variables
[2128]525         & frld
526#endif
[3294]527#if defined key_lim2
[6140]528      USE ice_2, ONLY : &          ! LIM2 Ice model variables
[3294]529         & frld
[2715]530#endif
[2128]531      IMPLICIT NONE
532
533      !! * Arguments
[6140]534      INTEGER, INTENT(IN) :: kstp  ! Current timestep
[2128]535      !! * Local declarations
[6140]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
[7910]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
[7911]549#if ! defined key_lim2 && ! defined key_lim3
[7910]550      REAL(wp), DIMENSION(jpi,jpj) :: frld
[3294]551#endif
[6140]552      LOGICAL :: llnightav        ! Logical for calculating night-time average
553
554      !Allocate local work arrays
[3294]555#if ! defined key_lim2 && ! defined key_lim3
[2715]556#endif
557
[2128]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      !-----------------------------------------------------------------------
[3294]569#if ! defined key_lim2 && ! defined key_lim3
[2128]570      frld(:,:) = 0.0_wp
571#endif
572      !-----------------------------------------------------------------------
[6140]573      ! Call the profile and surface observation operators
[2128]574      !-----------------------------------------------------------------------
575
[6140]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 ) 
[2128]617            ELSE
[6140]618               CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// &
619                         'yet working for velocity data (turn off velocity observations')
[2128]620            ENDIF
[6140]621
[2128]622         END DO
623
624      ENDIF
625
[6140]626      IF ( nsurftypes > 0 ) THEN
[2128]627
[6140]628         DO jtype = 1, nsurftypes
[2128]629
[6140]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.
[3294]637#if defined key_lim2 || defined key_lim3
[6140]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.
[2128]654#endif
[6140]655            END SELECT
[2128]656
[6140]657            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
658               &               nit000, idaystp, zsurfvar, tmask(:,:,1), &
659               &               nn_2dint, llnightav )
660
[2128]661         END DO
[6140]662
[2128]663      ENDIF
664
[3294]665#if ! defined key_lim2 && ! defined key_lim3
[2715]666#endif
667
[2128]668   END SUBROUTINE dia_obs
[6140]669
670   SUBROUTINE dia_obs_wri
[2128]671      !!----------------------------------------------------------------------
672      !!                    ***  ROUTINE dia_obs_wri  ***
673      !!         
674      !! ** Purpose : Call observation diagnostic output routines
675      !!
676      !! ** Method  : Call observation diagnostic output routines
677      !!
[6140]678      !! ** Action  :
[2128]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
[6140]686      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
[2128]687      !!----------------------------------------------------------------------
[6140]688      !! * Modules used
689      USE obs_rot_vel          ! Rotation of velocities
690
[2128]691      IMPLICIT NONE
692
693      !! * Local declarations
[6140]694      INTEGER :: jtype                    ! Data set loop variable
695      INTEGER :: jo, jvar, jk
696      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
697         & zu, &
698         & zv
[2128]699
700      !-----------------------------------------------------------------------
701      ! Depending on switches call various observation output routines
702      !-----------------------------------------------------------------------
703
[6140]704      IF ( nproftypes > 0 ) THEN
[2128]705
[6140]706         DO jtype = 1, nproftypes
[2128]707
[6140]708            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
[2128]709
[6140]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                  & )
[2128]716
[6140]717               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
[2128]718
[6140]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)
[2128]722
[6140]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
[2128]728
[6140]729                     END DO
730                  END DO
731               END DO
[2128]732
[6140]733               DEALLOCATE( zu )
734               DEALLOCATE( zv )
[2128]735
[6140]736            END IF
[2128]737
[6140]738            CALL obs_prof_decompress( profdataqc(jtype), &
739               &                      profdata(jtype), .TRUE., numout )
[2128]740
[6140]741            CALL obs_wri_prof( profdata(jtype) )
[2128]742
743         END DO
744
745      ENDIF
746
[6140]747      IF ( nsurftypes > 0 ) THEN
[2128]748
[6140]749         DO jtype = 1, nsurftypes
[2128]750
[6140]751            CALL obs_surf_decompress( surfdataqc(jtype), &
752               &                      surfdata(jtype), .TRUE., numout )
[2128]753
[6140]754            CALL obs_wri_surf( surfdata(jtype) )
[2128]755
756         END DO
757
758      ENDIF
759
760   END SUBROUTINE dia_obs_wri
761
[4245]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      !!----------------------------------------------------------------------
[6140]775      ! obs_grid deallocation
[4245]776      CALL obs_grid_deallocate
777
[6140]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
[4245]785   END SUBROUTINE dia_obs_dealloc
786
[6140]787   SUBROUTINE calc_date( kstp, ddobs )
[2128]788      !!----------------------------------------------------------------------
[6140]789      !!                    ***  ROUTINE calc_date  ***
[2128]790      !!         
[6140]791      !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format
[2128]792      !!
[6140]793      !! ** Method  : Get date in double precision YYYYMMDD.HHMMSS format
[2128]794      !!
[6140]795      !! ** Action  : Get date in double precision YYYYMMDD.HHMMSS format
[2128]796      !!
[6140]797      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
798      !!
[2128]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
[6140]805      !!        !  2014-09  (D. Lea) New generic routine now deals with arbitrary initial time of day
[2128]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
[6140]815      REAL(KIND=dp), INTENT(OUT) :: ddobs                        ! Date in YYYYMMDD.HHMMSS
816      INTEGER :: kstp
[2128]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
[6140]824      INTEGER :: imday       ! Number of days in month.
825      REAL(wp) :: zdayfrc    ! Fraction of day
[2128]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
[6140]835      ihou =   nn_time0 / 100
836      imin = ( nn_time0 - ihou * 100 ) 
[2128]837
838      !!----------------------------------------------------------------------
839      !! Compute number of days + number of hours + min since initial time
840      !!----------------------------------------------------------------------
[6140]841      zdayfrc = kstp * rdt / rday
[2128]842      zdayfrc = zdayfrc - aint(zdayfrc)
[6140]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 
[2128]853
[6140]854      !-----------------------------------------------------------------------
855      ! Convert number of days (iday) into a real date
856      !----------------------------------------------------------------------
[2128]857
858      CALL calc_month_len( iyea, imonth_len )
[6140]859
[2128]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
[6140]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 )
[2128]879      !!----------------------------------------------------------------------
[6140]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
[2128]894      !!----------------------------------------------------------------------
895
[6140]896      IMPLICIT NONE
[2128]897
[6140]898      !! * Arguments
899      REAL(KIND=dp), INTENT(OUT) :: ddobsini                   ! Initial date in YYYYMMDD.HHMMSS
900
901      CALL calc_date( nit000 - 1, ddobsini )
902
[2128]903   END SUBROUTINE ini_date
904
905   SUBROUTINE fin_date( ddobsfin )
906      !!----------------------------------------------------------------------
907      !!                    ***  ROUTINE fin_date  ***
908      !!         
[6140]909      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
[2128]910      !!
[6140]911      !! ** Method  :
[2128]912      !!
[6140]913      !! ** Action  :
[2128]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
[6140]920      !!        !  2014-09  (D. Lea) Change to call generic routine calc_date
[2128]921      !!----------------------------------------------------------------------
922
923      IMPLICIT NONE
924
925      !! * Arguments
[6140]926      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
[2128]927
[6140]928      CALL calc_date( nitend, ddobsfin )
[2128]929
[6140]930   END SUBROUTINE fin_date
931   
[2128]932END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.