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

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

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 7862

Last change on this file since 7862 was 7862, checked in by mattmartin, 5 years ago

Change standard output writing.

File size: 59.6 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         & MEDUSA_CHL, &
649         & MEDUSA_FCO2, &
650         & MEDUSA_PCO2, &
651         & MEDUSA_FILL_FLT
652#elif defined key_fabm
653      USE fabm
654      USE par_fabm
655#endif
656#if defined key_spm
657      USE par_spm, ONLY: &         ! ERSEM/SPM sediments
658         & jp_spm
659      USE trc, ONLY :  &
660         & trn
661#endif
662
663      IMPLICIT NONE
664
665      !! * Arguments
666      INTEGER, INTENT(IN) :: kstp  ! Current timestep
667      !! * Local declarations
668      INTEGER :: idaystp           ! Number of timesteps per day
669      INTEGER :: jtype             ! Data loop variable
670      INTEGER :: jvar              ! Variable number
671      INTEGER :: ji, jj            ! Loop counters
672      REAL(wp) :: tiny                  ! small number
673      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
674         & zprofvar1, &            ! Model values for 1st variable in a prof ob
675         & zprofvar2               ! Model values for 2nd variable in a prof ob
676      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
677         & zprofmask1, &           ! Mask associated with zprofvar1
678         & zprofmask2              ! Mask associated with zprofvar2
679      REAL(wp), POINTER, DIMENSION(:,:) :: &
680         & zsurfvar, &             ! Model values equivalent to surface ob.
681         & zsurfmask               ! Mask associated with surface variable
682      REAL(wp), POINTER, DIMENSION(:,:) :: &
683         & zglam1,    &            ! Model longitudes for prof variable 1
684         & zglam2,    &            ! Model longitudes for prof variable 2
685         & zgphi1,    &            ! Model latitudes for prof variable 1
686         & zgphi2                  ! Model latitudes for prof variable 2
687
688
689      !Allocate local work arrays
690      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
691      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
692      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
693      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
694      CALL wrk_alloc( jpi, jpj, zsurfvar )
695      CALL wrk_alloc( jpi, jpj, zsurfmask )
696      CALL wrk_alloc( jpi, jpj, zglam1 )
697      CALL wrk_alloc( jpi, jpj, zglam2 )
698      CALL wrk_alloc( jpi, jpj, zgphi1 )
699      CALL wrk_alloc( jpi, jpj, zgphi2 )
700
701      IF(lwp) THEN
702         WRITE(numout,*)
703         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
704         WRITE(numout,*) '~~~~~~~'
705         CALL FLUSH(numout)
706      ENDIF
707
708      idaystp = NINT( rday / rdt )
709
710      !-----------------------------------------------------------------------
711      ! Call the profile and surface observation operators
712      !-----------------------------------------------------------------------
713
714      IF ( nproftypes > 0 ) THEN
715
716         DO jtype = 1, nproftypes
717
718            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
719            CASE('prof')
720               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
721               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
722               zprofmask1(:,:,:) = tmask(:,:,:)
723               zprofmask2(:,:,:) = tmask(:,:,:)
724               zglam1(:,:) = glamt(:,:)
725               zglam2(:,:) = glamt(:,:)
726               zgphi1(:,:) = gphit(:,:)
727               zgphi2(:,:) = gphit(:,:)
728            CASE('vel')
729               zprofvar1(:,:,:) = un(:,:,:)
730               zprofvar2(:,:,:) = vn(:,:,:)
731               zprofmask1(:,:,:) = umask(:,:,:)
732               zprofmask2(:,:,:) = vmask(:,:,:)
733               zglam1(:,:) = glamu(:,:)
734               zglam2(:,:) = glamv(:,:)
735               zgphi1(:,:) = gphiu(:,:)
736               zgphi2(:,:) = gphiv(:,:)
737            CASE DEFAULT
738               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
739            END SELECT
740
741            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
742               &               nit000, idaystp,                         &
743               &               zprofvar1, zprofvar2,                    &
744               &               fsdept(:,:,:), fsdepw(:,:,:),            & 
745               &               zprofmask1, zprofmask2,                  &
746               &               zglam1, zglam2, zgphi1, zgphi2,          &
747               &               nn_1dint, nn_2dint,                      &
748               &               kdailyavtypes = nn_profdavtypes )
749
750         END DO
751
752      ENDIF
753
754      IF ( nsurftypes > 0 ) THEN
755
756         DO jtype = 1, nsurftypes
757
758            !Defaults which might be changed
759            zsurfmask(:,:) = tmask(:,:,1)
760
761            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
762            CASE('sst')
763               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
764            CASE('sla')
765               zsurfvar(:,:) = sshn(:,:)
766            CASE('sss')
767               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
768            CASE('sic')
769               IF ( kstp == 0 ) THEN
770                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
771                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
772                        &           'time-step but some obs are valid then.' )
773                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
774                        &           ' sea-ice obs will be missed'
775                  ENDIF
776                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
777                     &                        surfdataqc(jtype)%nsstp(1)
778                  CYCLE
779               ELSE
780#if defined key_cice
781                  zsurfvar(:,:) = fr_i(:,:)
782#elif defined key_lim2 || defined key_lim3
783                  zsurfvar(:,:) = 1._wp - frld(:,:)
784#else
785               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
786                  &           ' but no sea-ice model appears to have been defined' )
787#endif
788               ENDIF
789
790            CASE('logchl')
791#if defined key_hadocc
792               zsurfvar(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
793#elif defined key_medusa && defined key_foam_medusa
794               zsurfvar(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
795#elif defined key_fabm
796               chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
797               zsurfvar(:,:) = chl_3d(:,:,1)
798#else
799               CALL ctl_stop( ' Trying to run logchl observation operator', &
800                  &           ' but no biogeochemical model appears to have been defined' )
801#endif
802               zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
803               ! Take the log10 where we can, otherwise exclude
804               tiny = 1.0e-20
805               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
806                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
807               ELSEWHERE
808                  zsurfvar(:,:)  = obfillflt
809                  zsurfmask(:,:) = 0
810               END WHERE
811            CASE('spm')
812#if defined key_spm
813               zsurfvar(:,:) = 0.0
814               DO jn = 1, jp_spm
815                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
816               END DO
817#else
818               CALL ctl_stop( ' Trying to run spm observation operator', &
819                  &           ' but no spm model appears to have been defined' )
820#endif
821            CASE('fco2')
822#if defined key_hadocc
823               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
824               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
825                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
826                  zsurfvar(:,:) = obfillflt
827                  zsurfmask(:,:) = 0
828                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
829                     &           ' on timestep ' // TRIM(STR(kstp)),                              &
830                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
831               ENDIF
832#elif defined key_medusa && defined key_foam_medusa
833               zsurfmask(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA
834               IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) .AND. &
835                  & ( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
836                  zsurfvar(:,:) = obfillflt
837                  zsurfmask(:,:) = 0
838                  CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
839                     &           ' on timestep ' // TRIM(STR(kstp)),                              &
840                     &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
841               ENDIF
842#elif defined key_fabm
843               ! First, get pCO2 from FABM
844               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
845               zsurfvar(:,:) = pco2_3d(:,:,1)
846               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
847               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
848               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
849               ! and
850               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
851               ! Marine Chemistry, 2: 203-215.
852               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
853               ! not explicitly included - atmospheric pressure is not necessarily available so this is
854               ! the best assumption.
855               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
856               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
857               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
858               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
859               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + &
860                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
861                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
862                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
863                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
864                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
865#else
866               CALL ctl_stop( ' Trying to run fco2 observation operator', &
867                  &           ' but no biogeochemical model appears to have been defined' )
868#endif
869            CASE('pco2')
870#if defined key_hadocc
871               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
872               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
873                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
874                  zsurfvar(:,:) = obfillflt
875                  zsurfmask(:,:) = 0
876                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
877                     &           ' on timestep ' // TRIM(STR(kstp)),                              &
878                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
879               ENDIF
880#elif defined key_medusa && defined key_foam_medusa
881               zsurfvar(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA
882               IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) .AND. &
883                  & ( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
884                  zsurfvar(:,:) = obfillflt
885                  zsurfmask(:,:) = 0
886                  CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
887                     &           ' on timestep ' // TRIM(STR(kstp)),                              &
888                     &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
889               ENDIF
890#elif defined key_fabm
891               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
892               zsurfvar(:,:) = pco2_3d(:,:,1)
893#else
894               CALL ctl_stop( ' Trying to run pCO2 observation operator', &
895                  &           ' but no biogeochemical model appears to have been defined' )
896#endif
897
898            CASE DEFAULT
899
900               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
901
902            END SELECT
903
904            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
905               &               nit000, idaystp, zsurfvar, zsurfmask,    &
906               &               n2dintsurf(jtype), llnightav(jtype),     &
907               &               ravglamscl(jtype), ravgphiscl(jtype),     &
908               &               lfpindegs(jtype) )
909
910         END DO
911
912      ENDIF
913
914      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
915      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
916      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
917      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
918      CALL wrk_dealloc( jpi, jpj, zsurfvar )
919      CALL wrk_dealloc( jpi, jpj, zsurfmask )
920      CALL wrk_dealloc( jpi, jpj, zglam1 )
921      CALL wrk_dealloc( jpi, jpj, zglam2 )
922      CALL wrk_dealloc( jpi, jpj, zgphi1 )
923      CALL wrk_dealloc( jpi, jpj, zgphi2 )
924
925   END SUBROUTINE dia_obs
926
927   SUBROUTINE dia_obs_wri
928      !!----------------------------------------------------------------------
929      !!                    ***  ROUTINE dia_obs_wri  ***
930      !!         
931      !! ** Purpose : Call observation diagnostic output routines
932      !!
933      !! ** Method  : Call observation diagnostic output routines
934      !!
935      !! ** Action  :
936      !!
937      !! History :
938      !!        !  06-03  (K. Mogensen) Original code
939      !!        !  06-05  (K. Mogensen) Reformatted
940      !!        !  06-10  (A. Weaver) Cleaning
941      !!        !  07-03  (K. Mogensen) General handling of profiles
942      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
943      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
944      !!----------------------------------------------------------------------
945      !! * Modules used
946      USE obs_rot_vel          ! Rotation of velocities
947
948      IMPLICIT NONE
949
950      !! * Local declarations
951      INTEGER :: jtype                    ! Data set loop variable
952      INTEGER :: jo, jvar, jk
953      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
954         & zu, &
955         & zv
956
957      !-----------------------------------------------------------------------
958      ! Depending on switches call various observation output routines
959      !-----------------------------------------------------------------------
960
961      IF ( nproftypes > 0 ) THEN
962
963         DO jtype = 1, nproftypes
964
965            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
966
967               ! For velocity data, rotate the model velocities to N/S, E/W
968               ! using the compressed data structure.
969               ALLOCATE( &
970                  & zu(profdataqc(jtype)%nvprot(1)), &
971                  & zv(profdataqc(jtype)%nvprot(2))  &
972                  & )
973
974               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
975
976               DO jo = 1, profdataqc(jtype)%nprof
977                  DO jvar = 1, 2
978                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
979
980                        IF ( jvar == 1 ) THEN
981                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
982                        ELSE
983                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
984                        ENDIF
985
986                     END DO
987                  END DO
988               END DO
989
990               DEALLOCATE( zu )
991               DEALLOCATE( zv )
992
993            END IF
994
995            CALL obs_prof_decompress( profdataqc(jtype), &
996               &                      profdata(jtype), .TRUE., numout )
997
998            CALL obs_wri_prof( profdata(jtype) )
999
1000         END DO
1001
1002      ENDIF
1003
1004      IF ( nsurftypes > 0 ) THEN
1005
1006         DO jtype = 1, nsurftypes
1007
1008            CALL obs_surf_decompress( surfdataqc(jtype), &
1009               &                      surfdata(jtype), .TRUE., numout )
1010
1011            CALL obs_wri_surf( surfdata(jtype) )
1012
1013         END DO
1014
1015      ENDIF
1016
1017   END SUBROUTINE dia_obs_wri
1018
1019   SUBROUTINE dia_obs_dealloc
1020      IMPLICIT NONE
1021      !!----------------------------------------------------------------------
1022      !!                    *** ROUTINE dia_obs_dealloc ***
1023      !!
1024      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1025      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1026      !!
1027      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1028      !!
1029      !!  ** Action :
1030      !!
1031      !!----------------------------------------------------------------------
1032      ! obs_grid deallocation
1033      CALL obs_grid_deallocate
1034
1035      ! diaobs deallocation
1036      IF ( nproftypes > 0 ) &
1037         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
1038
1039      IF ( nsurftypes > 0 ) &
1040         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
1041         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
1042
1043   END SUBROUTINE dia_obs_dealloc
1044
1045   SUBROUTINE ini_date( ddobsini )
1046      !!----------------------------------------------------------------------
1047      !!                    ***  ROUTINE ini_date  ***
1048      !!
1049      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1050      !!
1051      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
1052      !!
1053      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
1054      !!
1055      !! History :
1056      !!        !  06-03  (K. Mogensen)  Original code
1057      !!        !  06-05  (K. Mogensen)  Reformatted
1058      !!        !  06-10  (A. Weaver) Cleaning
1059      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1060      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1061      !!----------------------------------------------------------------------
1062      USE phycst, ONLY : &            ! Physical constants
1063         & rday
1064      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1065         & rdt
1066
1067      IMPLICIT NONE
1068
1069      !! * Arguments
1070      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
1071
1072      !! * Local declarations
1073      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1074      INTEGER :: imon
1075      INTEGER :: iday
1076      INTEGER :: ihou
1077      INTEGER :: imin
1078      INTEGER :: imday       ! Number of days in month.
1079      INTEGER, DIMENSION(12) :: &
1080         &       imonth_len  ! Length in days of the months of the current year
1081      REAL(wp) :: zdayfrc    ! Fraction of day
1082
1083      !----------------------------------------------------------------------
1084      ! Initial date initialization (year, month, day, hour, minute)
1085      ! (This assumes that the initial date is for 00z))
1086      !----------------------------------------------------------------------
1087      iyea =   ndate0 / 10000
1088      imon = ( ndate0 - iyea * 10000 ) / 100
1089      iday =   ndate0 - iyea * 10000 - imon * 100
1090      ihou = 0
1091      imin = 0
1092
1093      !----------------------------------------------------------------------
1094      ! Compute number of days + number of hours + min since initial time
1095      !----------------------------------------------------------------------
1096      iday = iday + ( nit000 -1 ) * rdt / rday
1097      zdayfrc = ( nit000 -1 ) * rdt / rday
1098      zdayfrc = zdayfrc - aint(zdayfrc)
1099      ihou = int( zdayfrc * 24 )
1100      imin = int( (zdayfrc * 24 - ihou) * 60 )
1101
1102      !-----------------------------------------------------------------------
1103      ! Convert number of days (iday) into a real date
1104      !----------------------------------------------------------------------
1105
1106      CALL calc_month_len( iyea, imonth_len )
1107
1108      DO WHILE ( iday > imonth_len(imon) )
1109         iday = iday - imonth_len(imon)
1110         imon = imon + 1 
1111         IF ( imon > 12 ) THEN
1112            imon = 1
1113            iyea = iyea + 1
1114            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1115         ENDIF
1116      END DO
1117
1118      !----------------------------------------------------------------------
1119      ! Convert it into YYYYMMDD.HHMMSS format.
1120      !----------------------------------------------------------------------
1121      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1122         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1123
1124
1125   END SUBROUTINE ini_date
1126
1127   SUBROUTINE fin_date( ddobsfin )
1128      !!----------------------------------------------------------------------
1129      !!                    ***  ROUTINE fin_date  ***
1130      !!
1131      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1132      !!
1133      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
1134      !!
1135      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
1136      !!
1137      !! History :
1138      !!        !  06-03  (K. Mogensen)  Original code
1139      !!        !  06-05  (K. Mogensen)  Reformatted
1140      !!        !  06-10  (A. Weaver) Cleaning
1141      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1142      !!----------------------------------------------------------------------
1143      USE phycst, ONLY : &            ! Physical constants
1144         & rday
1145      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1146         & rdt
1147
1148      IMPLICIT NONE
1149
1150      !! * Arguments
1151      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1152
1153      !! * Local declarations
1154      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1155      INTEGER :: imon
1156      INTEGER :: iday
1157      INTEGER :: ihou
1158      INTEGER :: imin
1159      INTEGER :: imday       ! Number of days in month.
1160      INTEGER, DIMENSION(12) :: &
1161         &       imonth_len  ! Length in days of the months of the current year
1162      REAL(wp) :: zdayfrc    ! Fraction of day
1163
1164      !-----------------------------------------------------------------------
1165      ! Initial date initialization (year, month, day, hour, minute)
1166      ! (This assumes that the initial date is for 00z)
1167      !-----------------------------------------------------------------------
1168      iyea =   ndate0 / 10000
1169      imon = ( ndate0 - iyea * 10000 ) / 100
1170      iday =   ndate0 - iyea * 10000 - imon * 100
1171      ihou = 0
1172      imin = 0
1173
1174      !-----------------------------------------------------------------------
1175      ! Compute number of days + number of hours + min since initial time
1176      !-----------------------------------------------------------------------
1177      iday    = iday +  nitend  * rdt / rday
1178      zdayfrc =  nitend  * rdt / rday
1179      zdayfrc = zdayfrc - AINT( zdayfrc )
1180      ihou    = INT( zdayfrc * 24 )
1181      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1182
1183      !-----------------------------------------------------------------------
1184      ! Convert number of days (iday) into a real date
1185      !----------------------------------------------------------------------
1186
1187      CALL calc_month_len( iyea, imonth_len )
1188
1189      DO WHILE ( iday > imonth_len(imon) )
1190         iday = iday - imonth_len(imon)
1191         imon = imon + 1 
1192         IF ( imon > 12 ) THEN
1193            imon = 1
1194            iyea = iyea + 1
1195            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1196         ENDIF
1197      END DO
1198
1199      !-----------------------------------------------------------------------
1200      ! Convert it into YYYYMMDD.HHMMSS format
1201      !-----------------------------------------------------------------------
1202      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1203         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1204
1205    END SUBROUTINE fin_date
1206
1207    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
1208       &                         cfilestype, ifiles, cobstypes, cfiles )
1209
1210    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
1211    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
1212    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs
1213    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1214       &                   ifiles      ! Out appended number of files for this type
1215
1216    CHARACTER(len=6), INTENT(IN) :: ctypein 
1217    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
1218       &                   cfilestype  ! In list of files for this obs type
1219    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
1220       &                   cobstypes   ! Out appended list of obs types
1221    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
1222       &                   cfiles      ! Out appended list of files for all types
1223
1224    !Local variables
1225    INTEGER :: jfile
1226
1227    cfiles(jtype,:) = cfilestype(:)
1228    cobstypes(jtype) = ctypein
1229    ifiles(jtype) = 0
1230    DO jfile = 1, jpmaxnfiles
1231       IF ( trim(cfiles(jtype,jfile)) /= '' ) &
1232                 ifiles(jtype) = ifiles(jtype) + 1
1233    END DO
1234
1235    IF ( ifiles(jtype) == 0 ) THEN
1236         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   &
1237            &           ' set to true but no files available to read' )
1238    ENDIF
1239
1240    IF(lwp) THEN   
1241       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
1242       DO jfile = 1, ifiles(jtype)
1243          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
1244       END DO
1245    ENDIF
1246
1247    END SUBROUTINE obs_settypefiles
1248
1249    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
1250               &                  n2dint_default, n2dint_type,        &
1251               &                  ravglamscl_type, ravgphiscl_type,   &
1252               &                  lfp_indegs_type, lavnight_type,     &
1253               &                  n2dint, ravglamscl, ravgphiscl,     &
1254               &                  lfpindegs, lavnight )
1255
1256    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
1257    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
1258    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
1259    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
1260    REAL(wp), INTENT(IN) :: &
1261       &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
1262       &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
1263    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
1264    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
1265    CHARACTER(len=6), INTENT(IN) :: ctypein 
1266
1267    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1268       &                    n2dint 
1269    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1270       &                    ravglamscl, ravgphiscl
1271    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1272       &                    lfpindegs, lavnight
1273
1274    lavnight(jtype) = lavnight_type
1275
1276    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
1277       n2dint(jtype) = n2dint_type
1278    ELSE
1279       n2dint(jtype) = n2dint_default
1280    ENDIF
1281
1282    ! For averaging observation footprints set options for size of footprint
1283    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
1284       IF ( ravglamscl_type > 0._wp ) THEN
1285          ravglamscl(jtype) = ravglamscl_type
1286       ELSE
1287          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1288                         'scale (ravglamscl) for observation type '//TRIM(ctypein) )     
1289       ENDIF
1290
1291       IF ( ravgphiscl_type > 0._wp ) THEN
1292          ravgphiscl(jtype) = ravgphiscl_type
1293       ELSE
1294          CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1295                         'scale (ravgphiscl) for observation type '//TRIM(ctypein) )     
1296       ENDIF
1297
1298       lfpindegs(jtype) = lfp_indegs_type 
1299
1300    ENDIF
1301
1302    ! Write out info
1303    IF(lwp) THEN
1304       IF ( n2dint(jtype) <= 4 ) THEN
1305          WRITE(numout,*) '             '//TRIM(ctypein)// &
1306             &            ' model counterparts will be interpolated horizontally'
1307       ELSE IF ( n2dint(jtype) <= 6 ) THEN
1308          WRITE(numout,*) '             '//TRIM(ctypein)// &
1309             &            ' model counterparts will be averaged horizontally'
1310          WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
1311          WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
1312          IF ( lfpindegs(jtype) ) THEN
1313              WRITE(numout,*) '             '//'    (in degrees)'
1314          ELSE
1315              WRITE(numout,*) '             '//'    (in metres)'
1316          ENDIF
1317       ENDIF
1318    ENDIF
1319
1320    END SUBROUTINE obs_setinterpopts
1321
1322END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.