source: branches/UKMO/dev_r5518_obs_oper_update_medusa/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8455

Last change on this file since 8455 was 8455, checked in by dford, 3 years ago

Update getting chlorophyll, pCO2 and fCO2 from MEDUSA.

File size: 58.5 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   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS
13   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE wrk_nemo                 ! Memory Allocation
17   USE par_kind                 ! Precision variables
18   USE in_out_manager           ! I/O manager
19   USE par_oce
20   USE dom_oce                  ! Ocean space and time domain variables
21   USE obs_read_prof            ! Reading and allocation of profile obs
22   USE obs_read_surf            ! Reading and allocation of surface obs
23   USE obs_readmdt              ! Reading and allocation of MDT for SLA.
24   USE obs_prep                 ! Preparation of obs. (grid search etc).
25   USE obs_oper                 ! Observation operators
26   USE obs_write                ! Writing of observation related diagnostics
27   USE obs_grid                 ! Grid searching
28   USE obs_read_altbias         ! Bias treatment for altimeter
29   USE obs_sstbias              ! Bias correction routine for SST
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
45   !! * Module variables
46   LOGICAL, PUBLIC :: &
47      &       lk_diaobs = .TRUE.  !: Include this for backwards compatibility at NEMO 3.6.
48   LOGICAL :: ln_diaobs           !: Logical switch for the obs operator
49   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs
50   LOGICAL :: ln_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres
51   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres
52   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres
53   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
54
55   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
56   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
57   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
58   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
59   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
60   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
61   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
62   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
63
64   INTEGER :: nn_1dint       !: Vertical interpolation method
65   INTEGER :: nn_2dint       !: Default horizontal interpolation method
66   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method
67   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method
68   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method
69   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method
70 
71   INTEGER, DIMENSION(imaxavtypes) :: &
72      & nn_profdavtypes      !: Profile data types representing a daily average
73   INTEGER :: nproftypes     !: Number of profile obs types
74   INTEGER :: nsurftypes     !: Number of surface obs types
75   INTEGER, DIMENSION(:), ALLOCATABLE :: &
76      & nvarsprof, &         !: Number of profile variables
77      & nvarssurf            !: Number of surface variables
78   INTEGER, DIMENSION(:), ALLOCATABLE :: &
79      & nextrprof, &         !: Number of profile extra variables
80      & nextrsurf            !: Number of surface extra variables
81   INTEGER, DIMENSION(:), ALLOCATABLE :: &
82      & n2dintsurf           !: Interpolation option for surface variables
83   REAL(wp), DIMENSION(:), ALLOCATABLE :: &
84      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables
85      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables
86   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
87      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres
88      & llnightav            !: Logical for calculating night-time averages
89
90   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
91      & surfdata, &          !: Initial surface data
92      & surfdataqc           !: Surface data after quality control
93   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
94      & profdata, &          !: Initial profile data
95      & profdataqc           !: Profile data after quality control
96
97   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
98      & cobstypesprof, &     !: Profile obs types
99      & cobstypessurf        !: Surface obs types
100
101   !!----------------------------------------------------------------------
102   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
103   !! $Id$
104   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
105   !!----------------------------------------------------------------------
106
107   !! * Substitutions
108#  include "domzgr_substitute.h90"
109CONTAINS
110
111   SUBROUTINE dia_obs_init
112      !!----------------------------------------------------------------------
113      !!                    ***  ROUTINE dia_obs_init  ***
114      !!         
115      !! ** Purpose : Initialize and read observations
116      !!
117      !! ** Method  : Read the namelist and call reading routines
118      !!
119      !! ** Action  : Read the namelist and call reading routines
120      !!
121      !! History :
122      !!        !  06-03  (K. Mogensen) Original code
123      !!        !  06-05  (A. Weaver) Reformatted
124      !!        !  06-10  (A. Weaver) Cleaning and add controls
125      !!        !  07-03  (K. Mogensen) General handling of profiles
126      !!        !  14-08  (J.While) Incorporated SST bias correction
127      !!        !  15-02  (M. Martin) Simplification of namelist and code
128      !!----------------------------------------------------------------------
129
130      IMPLICIT NONE
131
132      !! * Local declarations
133      INTEGER, PARAMETER :: &
134         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
135      INTEGER, DIMENSION(:), ALLOCATABLE :: &
136         & ifilesprof, &         ! Number of profile files
137         & ifilessurf            ! Number of surface files
138      INTEGER :: ios             ! Local integer output status for namelist read
139      INTEGER :: jtype           ! Counter for obs types
140      INTEGER :: jvar            ! Counter for variables
141      INTEGER :: jfile           ! Counter for files
142      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply
143
144      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
145         & cn_profbfiles,    &   ! T/S profile input filenames
146         & cn_sstfbfiles,    &   ! Sea surface temperature input filenames
147         & cn_slafbfiles,    &   ! Sea level anomaly input filenames
148         & cn_sicfbfiles,    &   ! Seaice concentration input filenames
149         & cn_velfbfiles,    &   ! Velocity profile input filenames
150         & cn_sssfbfiles,    &   ! Sea surface salinity input filenames
151         & cn_logchlfbfiles, &   ! Log(Chl) input filenames
152         & cn_spmfbfiles,    &   ! Sediment input filenames
153         & cn_fco2fbfiles,   &   ! fco2 input filenames
154         & cn_pco2fbfiles,   &   ! pco2 input filenames
155         & cn_sstbiasfiles       ! SST bias input filenames
156
157      CHARACTER(LEN=128) :: &
158         & cn_altbiasfile        ! Altimeter bias input filename
159
160
161      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
162      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
163      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies
164      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
165      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
166      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs
167      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
168      LOGICAL :: ln_logchl       ! Logical switch for log(Chl) obs
169      LOGICAL :: ln_spm          ! Logical switch for sediment obs
170      LOGICAL :: ln_fco2         ! Logical switch for fco2 obs
171      LOGICAL :: ln_pco2         ! Logical switch for pco2 obs
172      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
173      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
174      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST
175      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
176      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
177      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
178
179      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
180      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
181
182      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
183         & clproffiles, &        ! Profile filenames
184         & clsurffiles           ! Surface filenames
185
186      LOGICAL :: llvar1          ! Logical for profile variable 1
187      LOGICAL :: llvar2          ! Logical for profile variable 1
188
189      REAL(wp), POINTER, DIMENSION(:,:) :: &
190         & zglam1, &             ! Model longitudes for profile variable 1
191         & zglam2                ! Model longitudes for profile variable 2
192      REAL(wp), POINTER, DIMENSION(:,:) :: &
193         & zgphi1, &             ! Model latitudes for profile variable 1
194         & zgphi2                ! Model latitudes for profile variable 2
195      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
196         & zmask1, &             ! Model land/sea mask associated with variable 1
197         & zmask2                ! Model land/sea mask associated with variable 2
198
199
200      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
201         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
202         &            ln_logchl, ln_spm, ln_fco2, ln_pco2,            &
203         &            ln_altbias, ln_sstbias, ln_nea,                 &
204         &            ln_grid_global, ln_grid_search_lookup,          &
205         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
206         &            ln_sstnight,                                    &
207         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
208         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
209         &            cn_profbfiles, cn_slafbfiles,                   &
210         &            cn_sstfbfiles, cn_sicfbfiles,                   &
211         &            cn_velfbfiles, cn_sssfbfiles,                   &
212         &            cn_logchlfbfiles, cn_spmfbfiles,                &
213         &            cn_fco2fbfiles, cn_pco2fbfiles,                 &
214         &            cn_sstbiasfiles, cn_altbiasfile,                &
215         &            cn_gridsearchfile, rn_gridsearchres,            &
216         &            rn_dobsini, rn_dobsend,                         &
217         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
218         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
219         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
220         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
221         &            nn_1dint, nn_2dint,                             &
222         &            nn_2dint_sla, nn_2dint_sst,                     &
223         &            nn_2dint_sss, nn_2dint_sic,                     &
224         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
225         &            nn_profdavtypes
226
227      CALL wrk_alloc( jpi, jpj, zglam1 )
228      CALL wrk_alloc( jpi, jpj, zglam2 )
229      CALL wrk_alloc( jpi, jpj, zgphi1 )
230      CALL wrk_alloc( jpi, jpj, zgphi2 )
231      CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
232      CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
233
234      !-----------------------------------------------------------------------
235      ! Read namelist parameters
236      !-----------------------------------------------------------------------
237
238      ! Some namelist arrays need initialising
239      cn_profbfiles(:)    = ''
240      cn_slafbfiles(:)    = ''
241      cn_sstfbfiles(:)    = ''
242      cn_sicfbfiles(:)    = ''
243      cn_velfbfiles(:)    = ''
244      cn_sssfbfiles(:)    = ''
245      cn_logchlfbfiles(:) = ''
246      cn_spmfbfiles(:)    = ''
247      cn_fco2fbfiles(:)   = ''
248      cn_pco2fbfiles(:)   = ''
249      cn_sstbiasfiles(:)  = ''
250      nn_profdavtypes(:)  = -1
251
252      CALL ini_date( rn_dobsini )
253      CALL fin_date( rn_dobsend )
254
255      ! Read namelist namobs : control observation diagnostics
256      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
257      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
258901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
259
260      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
261      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
262902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
263      IF(lwm) WRITE ( numond, namobs )
264
265      lk_diaobs = .FALSE.
266#if defined key_diaobs
267      IF ( ln_diaobs ) lk_diaobs = .TRUE.
268#endif
269
270      IF ( .NOT. lk_diaobs ) THEN
271         IF(lwp) WRITE(numout,cform_war)
272         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs'
273         RETURN
274      ENDIF
275
276      IF(lwp) THEN
277         WRITE(numout,*)
278         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
279         WRITE(numout,*) '~~~~~~~~~~~~'
280         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
281         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
282         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
283         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
284         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
285         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
286         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
287         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
288         WRITE(numout,*) '             Logical switch for log(Chl) observations              ln_logchl = ', ln_logchl
289         WRITE(numout,*) '             Logical switch for SPM observations                      ln_spm = ', ln_spm
290         WRITE(numout,*) '             Logical switch for FCO2 observations                    ln_fco2 = ', ln_fco2
291         WRITE(numout,*) '             Logical switch for PCO2 observations                    ln_pco2 = ', ln_pco2
292         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global
293         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
294         IF (ln_grid_search_lookup) &
295            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
296         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
297         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
298         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
299         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint
300         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
301         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
302         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
303         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
304         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
305         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
306         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias
307         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
308         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
309         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
310      ENDIF
311      !-----------------------------------------------------------------------
312      ! Set up list of observation types to be used
313      ! and the files associated with each type
314      !-----------------------------------------------------------------------
315
316      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
317      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, &
318         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) )
319
320      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
321         IF(lwp) WRITE(numout,cform_war)
322         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
323            &                    ' are set to .FALSE. so turning off calls to dia_obs'
324         nwarn = nwarn + 1
325         lk_diaobs = .FALSE.
326         RETURN
327      ENDIF
328
329      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes
330      IF ( nproftypes > 0 ) THEN
331
332         ALLOCATE( cobstypesprof(nproftypes) )
333         ALLOCATE( ifilesprof(nproftypes) )
334         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
335
336         jtype = 0
337         IF (ln_t3d .OR. ln_s3d) THEN
338            jtype = jtype + 1
339            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', &
340               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
341         ENDIF
342         IF (ln_vel3d) THEN
343            jtype = jtype + 1
344            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', &
345               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
346         ENDIF
347
348      ENDIF
349
350      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
351      IF ( nsurftypes > 0 ) THEN
352
353         ALLOCATE( cobstypessurf(nsurftypes) )
354         ALLOCATE( ifilessurf(nsurftypes) )
355         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
356         ALLOCATE(n2dintsurf(nsurftypes))
357         ALLOCATE(ravglamscl(nsurftypes))
358         ALLOCATE(ravgphiscl(nsurftypes))
359         ALLOCATE(lfpindegs(nsurftypes))
360         ALLOCATE(llnightav(nsurftypes))
361
362         jtype = 0
363         IF (ln_sla) THEN
364            jtype = jtype + 1
365            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', &
366               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
367            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      &
368               &                  nn_2dint, nn_2dint_sla,             &
369               &                  rn_sla_avglamscl, rn_sla_avgphiscl, &
370               &                  ln_sla_fp_indegs, .FALSE.,          &
371               &                  n2dintsurf, ravglamscl, ravgphiscl, &
372               &                  lfpindegs, llnightav )
373         ENDIF
374         IF (ln_sst) THEN
375            jtype = jtype + 1
376            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', &
377               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
378            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      &
379               &                  nn_2dint, nn_2dint_sst,             &
380               &                  rn_sst_avglamscl, rn_sst_avgphiscl, &
381               &                  ln_sst_fp_indegs, ln_sstnight,      &
382               &                  n2dintsurf, ravglamscl, ravgphiscl, &
383               &                  lfpindegs, llnightav )
384         ENDIF
385#if defined key_lim2 || defined key_lim3 || defined key_cice
386         IF (ln_sic) THEN
387            jtype = jtype + 1
388            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', &
389               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
390            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      &
391               &                  nn_2dint, nn_2dint_sic,             &
392               &                  rn_sic_avglamscl, rn_sic_avgphiscl, &
393               &                  ln_sic_fp_indegs, .FALSE.,          &
394               &                  n2dintsurf, ravglamscl, ravgphiscl, &
395               &                  lfpindegs, llnightav )
396         ENDIF
397#endif
398         IF (ln_sss) THEN
399            jtype = jtype + 1
400            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', &
401               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
402            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      &
403               &                  nn_2dint, nn_2dint_sss,             &
404               &                  rn_sss_avglamscl, rn_sss_avgphiscl, &
405               &                  ln_sss_fp_indegs, .FALSE.,          &
406               &                  n2dintsurf, ravglamscl, ravgphiscl, &
407               &                  lfpindegs, llnightav )
408         ENDIF
409
410         IF (ln_logchl) THEN
411            jtype = jtype + 1
412            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', &
413               &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles )
414            CALL obs_setinterpopts( nsurftypes, jtype, 'logchl',         &
415               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
416               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
417               &                    lfpindegs, llnightav )
418         ENDIF
419
420         IF (ln_spm) THEN
421            jtype = jtype + 1
422            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', &
423               &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles )
424            CALL obs_setinterpopts( nsurftypes, jtype, 'spm   ',         &
425               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
426               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
427               &                    lfpindegs, llnightav )
428         ENDIF
429
430         IF (ln_fco2) THEN
431            jtype = jtype + 1
432            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', &
433               &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
434            CALL obs_setinterpopts( nsurftypes, jtype, 'fco2  ',         &
435               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
436               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
437               &                    lfpindegs, llnightav )
438         ENDIF
439
440         IF (ln_pco2) THEN
441            jtype = jtype + 1
442            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', &
443               &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
444            CALL obs_setinterpopts( nsurftypes, jtype, 'pco2  ',         &
445               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
446               &                    n2dintsurf, ravglamscl, ravgphiscl,    &
447               &                    lfpindegs, llnightav )
448         ENDIF
449
450      ENDIF
451
452      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
453
454
455      !-----------------------------------------------------------------------
456      ! Obs operator parameter checking and initialisations
457      !-----------------------------------------------------------------------
458
459      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
460         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
461         RETURN
462      ENDIF
463
464      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
465         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
466            &                    ' is not available')
467      ENDIF
468
469      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN
470         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
471            &                    ' is not available')
472      ENDIF
473
474      CALL obs_typ_init
475
476      CALL mppmap_init
477
478      CALL obs_grid_setup( )
479
480      !-----------------------------------------------------------------------
481      ! Depending on switches read the various observation types
482      !-----------------------------------------------------------------------
483
484      IF ( nproftypes > 0 ) THEN
485
486         ALLOCATE(profdata(nproftypes))
487         ALLOCATE(profdataqc(nproftypes))
488         ALLOCATE(nvarsprof(nproftypes))
489         ALLOCATE(nextrprof(nproftypes))
490
491         DO jtype = 1, nproftypes
492
493            nvarsprof(jtype) = 2
494            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
495               nextrprof(jtype) = 1
496               llvar1 = ln_t3d
497               llvar2 = ln_s3d
498               zglam1 = glamt
499               zgphi1 = gphit
500               zmask1 = tmask
501               zglam2 = glamt
502               zgphi2 = gphit
503               zmask2 = tmask
504            ENDIF
505            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
506               nextrprof(jtype) = 2
507               llvar1 = ln_vel3d
508               llvar2 = ln_vel3d
509               zglam1 = glamu
510               zgphi1 = gphiu
511               zmask1 = umask
512               zglam2 = glamv
513               zgphi2 = gphiv
514               zmask2 = vmask
515            ENDIF
516
517            !Read in profile or profile obs types
518            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
519               &               clproffiles(jtype,1:ifilesprof(jtype)), &
520               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
521               &               rn_dobsini, rn_dobsend, llvar1, llvar2, &
522               &               ln_ignmis, ln_s_at_t, .FALSE., &
523               &               kdailyavtypes = nn_profdavtypes )
524
525            DO jvar = 1, nvarsprof(jtype)
526               CALL obs_prof_staend( profdata(jtype), jvar )
527            END DO
528
529            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
530               &               llvar1, llvar2, &
531               &               jpi, jpj, jpk, &
532               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  &
533               &               ln_nea, ln_bound_reject, &
534               &               kdailyavtypes = nn_profdavtypes )
535
536         END DO
537
538         DEALLOCATE( ifilesprof, clproffiles )
539
540      ENDIF
541
542      IF ( nsurftypes > 0 ) THEN
543
544         ALLOCATE(surfdata(nsurftypes))
545         ALLOCATE(surfdataqc(nsurftypes))
546         ALLOCATE(nvarssurf(nsurftypes))
547         ALLOCATE(nextrsurf(nsurftypes))
548
549         DO jtype = 1, nsurftypes
550
551            nvarssurf(jtype) = 1
552            nextrsurf(jtype) = 0
553            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
554
555            !Read in surface obs types
556            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
557               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
558               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
559               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
560
561            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
562
563            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
564               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
565               IF ( ln_altbias ) &
566                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
567            ENDIF
568
569            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
570               jnumsstbias = 0
571               DO jfile = 1, jpmaxnfiles
572                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
573                     &  jnumsstbias = jnumsstbias + 1
574               END DO
575               IF ( jnumsstbias == 0 ) THEN
576                  CALL ctl_stop("ln_sstbias set but no bias files to read in")   
577               ENDIF
578
579               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
580                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
581
582            ENDIF
583
584         END DO
585
586         DEALLOCATE( ifilessurf, clsurffiles )
587
588      ENDIF
589
590      CALL wrk_dealloc( jpi, jpj, zglam1 )
591      CALL wrk_dealloc( jpi, jpj, zglam2 )
592      CALL wrk_dealloc( jpi, jpj, zgphi1 )
593      CALL wrk_dealloc( jpi, jpj, zgphi2 )
594      CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
595      CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
596
597   END SUBROUTINE dia_obs_init
598
599   SUBROUTINE dia_obs( kstp )
600      !!----------------------------------------------------------------------
601      !!                    ***  ROUTINE dia_obs  ***
602      !!         
603      !! ** Purpose : Call the observation operators on each time step
604      !!
605      !! ** Method  : Call the observation operators on each time step to
606      !!              compute the model equivalent of the following data:
607      !!               - Profile data, currently T/S or U/V
608      !!               - Surface data, currently SST, SLA or sea-ice concentration.
609      !!
610      !! ** Action  :
611      !!
612      !! History :
613      !!        !  06-03  (K. Mogensen) Original code
614      !!        !  06-05  (K. Mogensen) Reformatted
615      !!        !  06-10  (A. Weaver) Cleaning
616      !!        !  07-03  (K. Mogensen) General handling of profiles
617      !!        !  07-04  (G. Smith) Generalized surface operators
618      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
619      !!        !  15-08  (M. Martin) Combined surface/profile routines.
620      !!----------------------------------------------------------------------
621      !! * Modules used
622      USE phycst, ONLY : &         ! Physical constants
623         & rday
624      USE oce, ONLY : &            ! Ocean dynamics and tracers variables
625         & tsn,       &
626         & un,        &
627         & vn,        &
628         & sshn
629#if defined  key_lim3
630      USE ice, ONLY : &            ! LIM3 Ice model variables
631         & frld
632#endif
633#if defined key_lim2
634      USE ice_2, ONLY : &          ! LIM2 Ice model variables
635         & frld
636#endif
637#if defined key_cice
638      USE sbc_oce, ONLY : fr_i     ! ice fraction
639#endif
640#if defined key_hadocc
641      USE trc, ONLY :  &           ! HadOCC chlorophyll, fCO2 and pCO2
642         & HADOCC_CHL, &
643         & HADOCC_FCO2, &
644         & HADOCC_PCO2, &
645         & HADOCC_FILL_FLT
646#elif defined key_medusa && defined key_foam_medusa
647      USE trc, ONLY :  &           ! MEDUSA chlorophyll, fCO2 and pCO2
648         & trn
649      USE par_medusa, ONLY: &
650         & jpchn, &
651         & jpchd
652#if defined key_roam
653      USE sms_medusa, ONLY: &
654         & f2_pco2w, &
655         & f2_fco2w
656#endif
657#elif defined key_fabm
658      USE fabm
659      USE par_fabm
660#endif
661#if defined key_spm
662      USE par_spm, ONLY: &         ! ERSEM/SPM sediments
663         & jp_spm
664      USE trc, ONLY :  &
665         & trn
666#endif
667
668      IMPLICIT NONE
669
670      !! * Arguments
671      INTEGER, INTENT(IN) :: kstp  ! Current timestep
672      !! * Local declarations
673      INTEGER :: idaystp           ! Number of timesteps per day
674      INTEGER :: jtype             ! Data loop variable
675      INTEGER :: jvar              ! Variable number
676      INTEGER :: ji, jj            ! Loop counters
677      REAL(wp) :: tiny                  ! small number
678      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
679         & zprofvar1, &            ! Model values for 1st variable in a prof ob
680         & zprofvar2               ! Model values for 2nd variable in a prof ob
681      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
682         & zprofmask1, &           ! Mask associated with zprofvar1
683         & zprofmask2              ! Mask associated with zprofvar2
684      REAL(wp), POINTER, DIMENSION(:,:) :: &
685         & zsurfvar, &             ! Model values equivalent to surface ob.
686         & zsurfmask               ! Mask associated with surface variable
687      REAL(wp), POINTER, DIMENSION(:,:) :: &
688         & zglam1,    &            ! Model longitudes for prof variable 1
689         & zglam2,    &            ! Model longitudes for prof variable 2
690         & zgphi1,    &            ! Model latitudes for prof variable 1
691         & zgphi2                  ! Model latitudes for prof variable 2
692
693
694      !Allocate local work arrays
695      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
696      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
697      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
698      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
699      CALL wrk_alloc( jpi, jpj, zsurfvar )
700      CALL wrk_alloc( jpi, jpj, zsurfmask )
701      CALL wrk_alloc( jpi, jpj, zglam1 )
702      CALL wrk_alloc( jpi, jpj, zglam2 )
703      CALL wrk_alloc( jpi, jpj, zgphi1 )
704      CALL wrk_alloc( jpi, jpj, zgphi2 )
705
706      IF(lwp) THEN
707         WRITE(numout,*)
708         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
709         WRITE(numout,*) '~~~~~~~'
710         CALL FLUSH(numout)
711      ENDIF
712
713      idaystp = NINT( rday / rdt )
714
715      !-----------------------------------------------------------------------
716      ! Call the profile and surface observation operators
717      !-----------------------------------------------------------------------
718
719      IF ( nproftypes > 0 ) THEN
720
721         DO jtype = 1, nproftypes
722
723            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
724            CASE('prof')
725               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
726               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
727               zprofmask1(:,:,:) = tmask(:,:,:)
728               zprofmask2(:,:,:) = tmask(:,:,:)
729               zglam1(:,:) = glamt(:,:)
730               zglam2(:,:) = glamt(:,:)
731               zgphi1(:,:) = gphit(:,:)
732               zgphi2(:,:) = gphit(:,:)
733            CASE('vel')
734               zprofvar1(:,:,:) = un(:,:,:)
735               zprofvar2(:,:,:) = vn(:,:,:)
736               zprofmask1(:,:,:) = umask(:,:,:)
737               zprofmask2(:,:,:) = vmask(:,:,:)
738               zglam1(:,:) = glamu(:,:)
739               zglam2(:,:) = glamv(:,:)
740               zgphi1(:,:) = gphiu(:,:)
741               zgphi2(:,:) = gphiv(:,:)
742            CASE DEFAULT
743               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
744            END SELECT
745
746            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
747               &               nit000, idaystp,                         &
748               &               zprofvar1, zprofvar2,                    &
749               &               fsdept(:,:,:), fsdepw(:,:,:),            & 
750               &               zprofmask1, zprofmask2,                  &
751               &               zglam1, zglam2, zgphi1, zgphi2,          &
752               &               nn_1dint, nn_2dint,                      &
753               &               kdailyavtypes = nn_profdavtypes )
754
755         END DO
756
757      ENDIF
758
759      IF ( nsurftypes > 0 ) THEN
760
761         DO jtype = 1, nsurftypes
762
763            !Defaults which might be changed
764            zsurfmask(:,:) = tmask(:,:,1)
765
766            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
767            CASE('sst')
768               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
769            CASE('sla')
770               zsurfvar(:,:) = sshn(:,:)
771            CASE('sss')
772               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
773            CASE('sic')
774               IF ( kstp == 0 ) THEN
775                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
776                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
777                        &           'time-step but some obs are valid then.' )
778                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
779                        &           ' sea-ice obs will be missed'
780                  ENDIF
781                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
782                     &                        surfdataqc(jtype)%nsstp(1)
783                  CYCLE
784               ELSE
785#if defined key_cice
786                  zsurfvar(:,:) = fr_i(:,:)
787#elif defined key_lim2 || defined key_lim3
788                  zsurfvar(:,:) = 1._wp - frld(:,:)
789#else
790               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
791                  &           ' but no sea-ice model appears to have been defined' )
792#endif
793               ENDIF
794
795            CASE('logchl')
796#if defined key_hadocc
797               zsurfvar(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
798#elif defined key_medusa && defined key_foam_medusa
799               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
800#elif defined key_fabm
801               chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
802               zsurfvar(:,:) = chl_3d(:,:,1)
803#else
804               CALL ctl_stop( ' Trying to run logchl observation operator', &
805                  &           ' but no biogeochemical model appears to have been defined' )
806#endif
807               zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
808               ! Take the log10 where we can, otherwise exclude
809               tiny = 1.0e-20
810               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
811                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
812               ELSEWHERE
813                  zsurfvar(:,:)  = obfillflt
814                  zsurfmask(:,:) = 0
815               END WHERE
816            CASE('spm')
817#if defined key_spm
818               zsurfvar(:,:) = 0.0
819               DO jn = 1, jp_spm
820                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
821               END DO
822#else
823               CALL ctl_stop( ' Trying to run spm observation operator', &
824                  &           ' but no spm model appears to have been defined' )
825#endif
826            CASE('fco2')
827#if defined key_hadocc
828               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
829               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
830                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
831                  zsurfvar(:,:) = obfillflt
832                  zsurfmask(:,:) = 0
833                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
834                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
835               ENDIF
836#elif defined key_medusa && defined key_foam_medusa && defined key_roam
837               zsurfvar(:,:) = f2_fco2w(:,:)
838#elif defined key_fabm
839               ! First, get pCO2 from FABM
840               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
841               zsurfvar(:,:) = pco2_3d(:,:,1)
842               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
843               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
844               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
845               ! and
846               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
847               ! Marine Chemistry, 2: 203-215.
848               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
849               ! not explicitly included - atmospheric pressure is not necessarily available so this is
850               ! the best assumption.
851               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
852               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
853               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
854               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
855               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + &
856                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
857                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
858                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
859                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
860                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
861#else
862               CALL ctl_stop( ' Trying to run fco2 observation operator', &
863                  &           ' but no biogeochemical model appears to have been defined' )
864#endif
865            CASE('pco2')
866#if defined key_hadocc
867               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
868               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
869                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
870                  zsurfvar(:,:) = obfillflt
871                  zsurfmask(:,:) = 0
872                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
873                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
874               ENDIF
875#elif defined key_medusa && defined key_foam_medusa && defined key_roam
876               zsurfvar(:,:) = f2_pco2w(:,:)
877#elif defined key_fabm
878               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
879               zsurfvar(:,:) = pco2_3d(:,:,1)
880#else
881               CALL ctl_stop( ' Trying to run pCO2 observation operator', &
882                  &           ' but no biogeochemical model appears to have been defined' )
883#endif
884
885            CASE DEFAULT
886
887               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
888
889            END SELECT
890
891            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
892               &               nit000, idaystp, zsurfvar, zsurfmask,    &
893               &               n2dintsurf(jtype), llnightav(jtype),     &
894               &               ravglamscl(jtype), ravgphiscl(jtype),     &
895               &               lfpindegs(jtype) )
896
897         END DO
898
899      ENDIF
900
901      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
902      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
903      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
904      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
905      CALL wrk_dealloc( jpi, jpj, zsurfvar )
906      CALL wrk_dealloc( jpi, jpj, zsurfmask )
907      CALL wrk_dealloc( jpi, jpj, zglam1 )
908      CALL wrk_dealloc( jpi, jpj, zglam2 )
909      CALL wrk_dealloc( jpi, jpj, zgphi1 )
910      CALL wrk_dealloc( jpi, jpj, zgphi2 )
911
912   END SUBROUTINE dia_obs
913
914   SUBROUTINE dia_obs_wri
915      !!----------------------------------------------------------------------
916      !!                    ***  ROUTINE dia_obs_wri  ***
917      !!         
918      !! ** Purpose : Call observation diagnostic output routines
919      !!
920      !! ** Method  : Call observation diagnostic output routines
921      !!
922      !! ** Action  :
923      !!
924      !! History :
925      !!        !  06-03  (K. Mogensen) Original code
926      !!        !  06-05  (K. Mogensen) Reformatted
927      !!        !  06-10  (A. Weaver) Cleaning
928      !!        !  07-03  (K. Mogensen) General handling of profiles
929      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
930      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
931      !!----------------------------------------------------------------------
932      !! * Modules used
933      USE obs_rot_vel          ! Rotation of velocities
934
935      IMPLICIT NONE
936
937      !! * Local declarations
938      INTEGER :: jtype                    ! Data set loop variable
939      INTEGER :: jo, jvar, jk
940      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
941         & zu, &
942         & zv
943
944      !-----------------------------------------------------------------------
945      ! Depending on switches call various observation output routines
946      !-----------------------------------------------------------------------
947
948      IF ( nproftypes > 0 ) THEN
949
950         DO jtype = 1, nproftypes
951
952            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
953
954               ! For velocity data, rotate the model velocities to N/S, E/W
955               ! using the compressed data structure.
956               ALLOCATE( &
957                  & zu(profdataqc(jtype)%nvprot(1)), &
958                  & zv(profdataqc(jtype)%nvprot(2))  &
959                  & )
960
961               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
962
963               DO jo = 1, profdataqc(jtype)%nprof
964                  DO jvar = 1, 2
965                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
966
967                        IF ( jvar == 1 ) THEN
968                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
969                        ELSE
970                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
971                        ENDIF
972
973                     END DO
974                  END DO
975               END DO
976
977               DEALLOCATE( zu )
978               DEALLOCATE( zv )
979
980            END IF
981
982            CALL obs_prof_decompress( profdataqc(jtype), &
983               &                      profdata(jtype), .TRUE., numout )
984
985            CALL obs_wri_prof( profdata(jtype) )
986
987         END DO
988
989      ENDIF
990
991      IF ( nsurftypes > 0 ) THEN
992
993         DO jtype = 1, nsurftypes
994
995            CALL obs_surf_decompress( surfdataqc(jtype), &
996               &                      surfdata(jtype), .TRUE., numout )
997
998            CALL obs_wri_surf( surfdata(jtype) )
999
1000         END DO
1001
1002      ENDIF
1003
1004   END SUBROUTINE dia_obs_wri
1005
1006   SUBROUTINE dia_obs_dealloc
1007      IMPLICIT NONE
1008      !!----------------------------------------------------------------------
1009      !!                    *** ROUTINE dia_obs_dealloc ***
1010      !!
1011      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1012      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1013      !!
1014      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1015      !!
1016      !!  ** Action :
1017      !!
1018      !!----------------------------------------------------------------------
1019      ! obs_grid deallocation
1020      CALL obs_grid_deallocate
1021
1022      ! diaobs deallocation
1023      IF ( nproftypes > 0 ) &
1024         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
1025
1026      IF ( nsurftypes > 0 ) &
1027         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
1028         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
1029
1030   END SUBROUTINE dia_obs_dealloc
1031
1032   SUBROUTINE ini_date( ddobsini )
1033      !!----------------------------------------------------------------------
1034      !!                    ***  ROUTINE ini_date  ***
1035      !!
1036      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1037      !!
1038      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
1039      !!
1040      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
1041      !!
1042      !! History :
1043      !!        !  06-03  (K. Mogensen)  Original code
1044      !!        !  06-05  (K. Mogensen)  Reformatted
1045      !!        !  06-10  (A. Weaver) Cleaning
1046      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1047      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1048      !!----------------------------------------------------------------------
1049      USE phycst, ONLY : &            ! Physical constants
1050         & rday
1051      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1052         & rdt
1053
1054      IMPLICIT NONE
1055
1056      !! * Arguments
1057      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
1058
1059      !! * Local declarations
1060      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1061      INTEGER :: imon
1062      INTEGER :: iday
1063      INTEGER :: ihou
1064      INTEGER :: imin
1065      INTEGER :: imday       ! Number of days in month.
1066      INTEGER, DIMENSION(12) :: &
1067         &       imonth_len  ! Length in days of the months of the current year
1068      REAL(wp) :: zdayfrc    ! Fraction of day
1069
1070      !----------------------------------------------------------------------
1071      ! Initial date initialization (year, month, day, hour, minute)
1072      ! (This assumes that the initial date is for 00z))
1073      !----------------------------------------------------------------------
1074      iyea =   ndate0 / 10000
1075      imon = ( ndate0 - iyea * 10000 ) / 100
1076      iday =   ndate0 - iyea * 10000 - imon * 100
1077      ihou = 0
1078      imin = 0
1079
1080      !----------------------------------------------------------------------
1081      ! Compute number of days + number of hours + min since initial time
1082      !----------------------------------------------------------------------
1083      iday = iday + ( nit000 -1 ) * rdt / rday
1084      zdayfrc = ( nit000 -1 ) * rdt / rday
1085      zdayfrc = zdayfrc - aint(zdayfrc)
1086      ihou = int( zdayfrc * 24 )
1087      imin = int( (zdayfrc * 24 - ihou) * 60 )
1088
1089      !-----------------------------------------------------------------------
1090      ! Convert number of days (iday) into a real date
1091      !----------------------------------------------------------------------
1092
1093      CALL calc_month_len( iyea, imonth_len )
1094
1095      DO WHILE ( iday > imonth_len(imon) )
1096         iday = iday - imonth_len(imon)
1097         imon = imon + 1 
1098         IF ( imon > 12 ) THEN
1099            imon = 1
1100            iyea = iyea + 1
1101            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1102         ENDIF
1103      END DO
1104
1105      !----------------------------------------------------------------------
1106      ! Convert it into YYYYMMDD.HHMMSS format.
1107      !----------------------------------------------------------------------
1108      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1109         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1110
1111
1112   END SUBROUTINE ini_date
1113
1114   SUBROUTINE fin_date( ddobsfin )
1115      !!----------------------------------------------------------------------
1116      !!                    ***  ROUTINE fin_date  ***
1117      !!
1118      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1119      !!
1120      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
1121      !!
1122      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
1123      !!
1124      !! History :
1125      !!        !  06-03  (K. Mogensen)  Original code
1126      !!        !  06-05  (K. Mogensen)  Reformatted
1127      !!        !  06-10  (A. Weaver) Cleaning
1128      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1129      !!----------------------------------------------------------------------
1130      USE phycst, ONLY : &            ! Physical constants
1131         & rday
1132      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1133         & rdt
1134
1135      IMPLICIT NONE
1136
1137      !! * Arguments
1138      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1139
1140      !! * Local declarations
1141      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1142      INTEGER :: imon
1143      INTEGER :: iday
1144      INTEGER :: ihou
1145      INTEGER :: imin
1146      INTEGER :: imday       ! Number of days in month.
1147      INTEGER, DIMENSION(12) :: &
1148         &       imonth_len  ! Length in days of the months of the current year
1149      REAL(wp) :: zdayfrc    ! Fraction of day
1150
1151      !-----------------------------------------------------------------------
1152      ! Initial date initialization (year, month, day, hour, minute)
1153      ! (This assumes that the initial date is for 00z)
1154      !-----------------------------------------------------------------------
1155      iyea =   ndate0 / 10000
1156      imon = ( ndate0 - iyea * 10000 ) / 100
1157      iday =   ndate0 - iyea * 10000 - imon * 100
1158      ihou = 0
1159      imin = 0
1160
1161      !-----------------------------------------------------------------------
1162      ! Compute number of days + number of hours + min since initial time
1163      !-----------------------------------------------------------------------
1164      iday    = iday +  nitend  * rdt / rday
1165      zdayfrc =  nitend  * rdt / rday
1166      zdayfrc = zdayfrc - AINT( zdayfrc )
1167      ihou    = INT( zdayfrc * 24 )
1168      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1169
1170      !-----------------------------------------------------------------------
1171      ! Convert number of days (iday) into a real date
1172      !----------------------------------------------------------------------
1173
1174      CALL calc_month_len( iyea, imonth_len )
1175
1176      DO WHILE ( iday > imonth_len(imon) )
1177         iday = iday - imonth_len(imon)
1178         imon = imon + 1 
1179         IF ( imon > 12 ) THEN
1180            imon = 1
1181            iyea = iyea + 1
1182            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1183         ENDIF
1184      END DO
1185
1186      !-----------------------------------------------------------------------
1187      ! Convert it into YYYYMMDD.HHMMSS format
1188      !-----------------------------------------------------------------------
1189      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1190         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1191
1192    END SUBROUTINE fin_date
1193
1194    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
1195       &                         cfilestype, ifiles, cobstypes, cfiles )
1196
1197    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
1198    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
1199    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs
1200    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1201       &                   ifiles      ! Out appended number of files for this type
1202
1203    CHARACTER(len=6), INTENT(IN) :: ctypein 
1204    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
1205       &                   cfilestype  ! In list of files for this obs type
1206    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
1207       &                   cobstypes   ! Out appended list of obs types
1208    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
1209       &                   cfiles      ! Out appended list of files for all types
1210
1211    !Local variables
1212    INTEGER :: jfile
1213
1214    cfiles(jtype,:) = cfilestype(:)
1215    cobstypes(jtype) = ctypein
1216    ifiles(jtype) = 0
1217    DO jfile = 1, jpmaxnfiles
1218       IF ( trim(cfiles(jtype,jfile)) /= '' ) &
1219                 ifiles(jtype) = ifiles(jtype) + 1
1220    END DO
1221
1222    IF ( ifiles(jtype) == 0 ) THEN
1223         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   &
1224            &           ' set to true but no files available to read' )
1225    ENDIF
1226
1227    IF(lwp) THEN   
1228       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
1229       DO jfile = 1, ifiles(jtype)
1230          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
1231       END DO
1232    ENDIF
1233
1234    END SUBROUTINE obs_settypefiles
1235
1236    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
1237               &                  n2dint_default, n2dint_type,        &
1238               &                  ravglamscl_type, ravgphiscl_type,   &
1239               &                  lfp_indegs_type, lavnight_type,     &
1240               &                  n2dint, ravglamscl, ravgphiscl,     &
1241               &                  lfpindegs, lavnight )
1242
1243    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
1244    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
1245    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
1246    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
1247    REAL(wp), INTENT(IN) :: &
1248       &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
1249       &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
1250    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
1251    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
1252    CHARACTER(len=6), INTENT(IN) :: ctypein 
1253
1254    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1255       &                    n2dint 
1256    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1257       &                    ravglamscl, ravgphiscl
1258    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1259       &                    lfpindegs, lavnight
1260
1261    lavnight(jtype) = lavnight_type
1262
1263    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
1264       n2dint(jtype) = n2dint_type
1265    ELSE
1266       n2dint(jtype) = n2dint_default
1267    ENDIF
1268
1269    ! For averaging observation footprints set options for size of footprint
1270    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
1271       IF ( ravglamscl_type > 0._wp ) THEN
1272          ravglamscl(jtype) = ravglamscl_type
1273       ELSE
1274          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1275                         'scale (ravglamscl) for observation type '//TRIM(ctypein) )     
1276       ENDIF
1277
1278       IF ( ravgphiscl_type > 0._wp ) THEN
1279          ravgphiscl(jtype) = ravgphiscl_type
1280       ELSE
1281          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1282                         'scale (ravgphiscl) for observation type '//TRIM(ctypein) )     
1283       ENDIF
1284
1285       lfpindegs(jtype) = lfp_indegs_type 
1286
1287    ENDIF
1288
1289    ! Write out info
1290    IF(lwp) THEN
1291       IF ( n2dint(jtype) <= 4 ) THEN
1292          WRITE(numout,*) '             '//TRIM(ctypein)// &
1293             &            ' model counterparts will be interpolated horizontally'
1294       ELSE IF ( n2dint(jtype) <= 6 ) THEN
1295          WRITE(numout,*) '             '//TRIM(ctypein)// &
1296             &            ' model counterparts will be averaged horizontally'
1297          WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
1298          WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
1299          IF ( lfpindegs(jtype) ) THEN
1300              WRITE(numout,*) '             '//'    (in degrees)'
1301          ELSE
1302              WRITE(numout,*) '             '//'    (in metres)'
1303          ENDIF
1304       ENDIF
1305    ENDIF
1306
1307    END SUBROUTINE obs_setinterpopts
1308
1309END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.