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

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9211

Last change on this file since 9211 was 9211, checked in by dford, 6 years ago

Change to loop over obs_prof_opt one variable at a time.

File size: 82.3 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_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres
51   LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres
52   LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres
53   LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres
54   LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
55
56   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint
57   REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint
58   REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint
59   REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint
60   REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint
61   REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint
62   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint
63   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint
64   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint
65   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint
66
67   INTEGER :: nn_1dint         !: Vertical interpolation method
68   INTEGER :: nn_2dint_default !: Default horizontal interpolation method
69   INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default)
70   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default)
71   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default)
72   INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default)
73 
74   INTEGER, DIMENSION(imaxavtypes) :: &
75      & nn_profdavtypes      !: Profile data types representing a daily average
76   INTEGER :: nproftypes     !: Number of profile obs types
77   INTEGER :: nsurftypes     !: Number of surface obs types
78   INTEGER, DIMENSION(:), ALLOCATABLE :: &
79      & nvarsprof, &         !: Number of profile variables
80      & nvarssurf            !: Number of surface variables
81   INTEGER, DIMENSION(:), ALLOCATABLE :: &
82      & nextrprof, &         !: Number of profile extra variables
83      & nextrsurf            !: Number of surface extra variables
84   INTEGER, DIMENSION(:), ALLOCATABLE :: &
85      & n2dintsurf           !: Interpolation option for surface variables
86   REAL(wp), DIMENSION(:), ALLOCATABLE :: &
87      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables
88      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables
89   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
90      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres
91      & llnightav            !: Logical for calculating night-time averages
92
93   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
94      & surfdata, &          !: Initial surface data
95      & surfdataqc           !: Surface data after quality control
96   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
97      & profdata, &          !: Initial profile data
98      & profdataqc           !: Profile data after quality control
99
100   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
101      & cobstypesprof, &     !: Profile obs types
102      & cobstypessurf        !: Surface obs types
103
104   !!----------------------------------------------------------------------
105   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
106   !! $Id$
107   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
108   !!----------------------------------------------------------------------
109
110   !! * Substitutions
111#  include "domzgr_substitute.h90"
112CONTAINS
113
114   SUBROUTINE dia_obs_init
115      !!----------------------------------------------------------------------
116      !!                    ***  ROUTINE dia_obs_init  ***
117      !!         
118      !! ** Purpose : Initialize and read observations
119      !!
120      !! ** Method  : Read the namelist and call reading routines
121      !!
122      !! ** Action  : Read the namelist and call reading routines
123      !!
124      !! History :
125      !!        !  06-03  (K. Mogensen) Original code
126      !!        !  06-05  (A. Weaver) Reformatted
127      !!        !  06-10  (A. Weaver) Cleaning and add controls
128      !!        !  07-03  (K. Mogensen) General handling of profiles
129      !!        !  14-08  (J.While) Incorporated SST bias correction
130      !!        !  15-02  (M. Martin) Simplification of namelist and code
131      !!----------------------------------------------------------------------
132
133      IMPLICIT NONE
134
135      !! * Local declarations
136      INTEGER, PARAMETER :: &
137         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
138      INTEGER, DIMENSION(:), ALLOCATABLE :: &
139         & ifilesprof, &         ! Number of profile files
140         & ifilessurf            ! Number of surface files
141      INTEGER :: ios             ! Local integer output status for namelist read
142      INTEGER :: jtype           ! Counter for obs types
143      INTEGER :: jvar            ! Counter for variables
144      INTEGER :: jfile           ! Counter for files
145      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply
146      INTEGER :: n2dint_type     ! Local version of nn_2dint*
147
148      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
149         & cn_profbfiles,      & ! T/S profile input filenames
150         & cn_sstfbfiles,      & ! Sea surface temperature input filenames
151         & cn_slafbfiles,      & ! Sea level anomaly input filenames
152         & cn_sicfbfiles,      & ! Seaice concentration input filenames
153         & cn_velfbfiles,      & ! Velocity profile input filenames
154         & cn_sssfbfiles,      & ! Sea surface salinity input filenames
155         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames
156         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames
157         & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames
158         & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames
159         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames
160         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames
161         & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames
162         & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames
163         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames
164         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames
165         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames
166         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames
167         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames
168         & cn_pno3fbfiles,     & ! Profile nitrate input filenames
169         & cn_psi4fbfiles,     & ! Profile silicate input filenames
170         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames
171         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames
172         & cn_palkfbfiles,     & ! Profile alkalinity input filenames
173         & cn_pphfbfiles,      & ! Profile pH input filenames
174         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames
175         & cn_sstbiasfiles       ! SST bias input filenames
176
177      CHARACTER(LEN=128) :: &
178         & cn_altbiasfile        ! Altimeter bias input filename
179
180
181      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
182      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
183      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies
184      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
185      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
186      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs
187      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
188      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs
189      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs
190      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs
191      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs
192      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs
193      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs
194      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs
195      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs
196      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs
197      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs
198      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs
199      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs
200      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs
201      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs
202      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs
203      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs
204      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs
205      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs
206      LOGICAL :: ln_pph          ! Logical switch for profile pH obs
207      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs
208      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
209      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
210      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST
211      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
212      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
213      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
214
215      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
216      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
217
218      REAL(wp) :: rtype_avglamscl ! Local version of rn_*_avglamscl
219      REAL(wp) :: rtype_avgphiscl ! Local version of rn_*_avgphiscl
220
221      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
222         & clproffiles, &        ! Profile filenames
223         & clsurffiles           ! Surface filenames
224
225      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read
226      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
227      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)
228
229      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
230         & zglam                 ! Model longitudes for profile variables
231      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
232         & zgphi                 ! Model latitudes for profile variables
233      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
234         & zmask                 ! Model land/sea mask associated with variables
235
236
237      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
238         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
239         &            ln_slchltot, ln_slchldia, ln_slchlnon,          &
240         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &
241         &            ln_slchlpic, ln_schltot,                        &
242         &            ln_sspm,     ln_sfco2,    ln_spco2,             &
243         &            ln_plchltot, ln_pchltot,  ln_pno3,              &
244         &            ln_psi4,     ln_ppo4,     ln_pdic,              &
245         &            ln_palk,     ln_pph,      ln_po2,               &
246         &            ln_altbias, ln_sstbias, ln_nea,                 &
247         &            ln_grid_global, ln_grid_search_lookup,          &
248         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
249         &            ln_sstnight, ln_default_fp_indegs,              &
250         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
251         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
252         &            cn_profbfiles, cn_slafbfiles,                   &
253         &            cn_sstfbfiles, cn_sicfbfiles,                   &
254         &            cn_velfbfiles, cn_sssfbfiles,                   &
255         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &
256         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         &
257         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         &
258         &            cn_slchlpicfbfiles,                             &
259         &            cn_schltotfbfiles, cn_sspmfbfiles,              &
260         &            cn_sfco2fbfiles, cn_spco2fbfiles,               &
261         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          &
262         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, &
263         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  &
264         &            cn_po2fbfiles,                                  &
265         &            cn_sstbiasfiles, cn_altbiasfile,                &
266         &            cn_gridsearchfile, rn_gridsearchres,            &
267         &            rn_dobsini, rn_dobsend,                         &
268         &            rn_default_avglamscl, rn_default_avgphiscl,     &
269         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
270         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
271         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
272         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
273         &            nn_1dint, nn_2dint_default,                     &
274         &            nn_2dint_sla, nn_2dint_sst,                     &
275         &            nn_2dint_sss, nn_2dint_sic,                     &
276         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
277         &            nn_profdavtypes
278
279      !-----------------------------------------------------------------------
280      ! Read namelist parameters
281      !-----------------------------------------------------------------------
282
283      ! Some namelist arrays need initialising
284      cn_profbfiles(:)      = ''
285      cn_slafbfiles(:)      = ''
286      cn_sstfbfiles(:)      = ''
287      cn_sicfbfiles(:)      = ''
288      cn_velfbfiles(:)      = ''
289      cn_sssfbfiles(:)      = ''
290      cn_slchltotfbfiles(:) = ''
291      cn_slchldiafbfiles(:) = ''
292      cn_slchlnonfbfiles(:) = ''
293      cn_slchldinfbfiles(:) = ''
294      cn_slchlmicfbfiles(:) = ''
295      cn_slchlnanfbfiles(:) = ''
296      cn_slchlpicfbfiles(:) = ''
297      cn_schltotfbfiles(:)  = ''
298      cn_sspmfbfiles(:)     = ''
299      cn_sfco2fbfiles(:)    = ''
300      cn_spco2fbfiles(:)    = ''
301      cn_plchltotfbfiles(:) = ''
302      cn_pchltotfbfiles(:)  = ''
303      cn_pno3fbfiles(:)     = ''
304      cn_psi4fbfiles(:)     = ''
305      cn_ppo4fbfiles(:)     = ''
306      cn_pdicfbfiles(:)     = ''
307      cn_palkfbfiles(:)     = ''
308      cn_pphfbfiles(:)      = ''
309      cn_po2fbfiles(:)      = ''
310      cn_sstbiasfiles(:)    = ''
311      nn_profdavtypes(:)    = -1
312
313      CALL ini_date( rn_dobsini )
314      CALL fin_date( rn_dobsend )
315
316      ! Read namelist namobs : control observation diagnostics
317      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
318      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
319901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
320
321      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
322      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
323902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
324      IF(lwm) WRITE ( numond, namobs )
325
326      lk_diaobs = .FALSE.
327#if defined key_diaobs
328      IF ( ln_diaobs ) lk_diaobs = .TRUE.
329#endif
330
331      IF ( .NOT. lk_diaobs ) THEN
332         IF(lwp) WRITE(numout,cform_war)
333         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs'
334         RETURN
335      ENDIF
336
337      IF(lwp) THEN
338         WRITE(numout,*)
339         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
340         WRITE(numout,*) '~~~~~~~~~~~~'
341         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
342         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
343         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
344         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
345         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
346         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
347         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
348         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
349         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot
350         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia
351         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon
352         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin
353         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic
354         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan
355         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic
356         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot
357         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm
358         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2
359         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2
360         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot
361         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot
362         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3
363         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4
364         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4
365         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic
366         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk
367         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph
368         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2
369         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global
370         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
371         IF (ln_grid_search_lookup) &
372            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
373         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
374         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
375         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
376         WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default
377         WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla
378         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst
379         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss
380         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic
381         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl
382         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl
383         WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs
384         WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl
385         WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl
386         WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs
387         WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl
388         WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl
389         WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs
390         WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl
391         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl
392         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs
393         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
394         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
395         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
396         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
397         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
398         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
399         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias
400         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
401         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
402         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
403      ENDIF
404      !-----------------------------------------------------------------------
405      ! Set up list of observation types to be used
406      ! and the files associated with each type
407      !-----------------------------------------------------------------------
408
409      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          &
410         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     &
411         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) )
412      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     &
413         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, &
414         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  &
415         &                  ln_sspm,     ln_sfco2,    ln_spco2 /) )
416
417      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
418         IF(lwp) WRITE(numout,cform_war)
419         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
420            &                    ' are set to .FALSE. so turning off calls to dia_obs'
421         nwarn = nwarn + 1
422         lk_diaobs = .FALSE.
423         RETURN
424      ENDIF
425
426      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes
427      IF ( nproftypes > 0 ) THEN
428
429         ALLOCATE( cobstypesprof(nproftypes) )
430         ALLOCATE( ifilesprof(nproftypes) )
431         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
432
433         jtype = 0
434         IF (ln_t3d .OR. ln_s3d) THEN
435            jtype = jtype + 1
436            cobstypesprof(jtype) = 'prof'
437            clproffiles(jtype,:) = cn_profbfiles
438         ENDIF
439         IF (ln_vel3d) THEN
440            jtype = jtype + 1
441            cobstypesprof(jtype) =  'vel'
442            clproffiles(jtype,:) = cn_velfbfiles
443         ENDIF
444         IF (ln_plchltot) THEN
445            jtype = jtype + 1
446            cobstypesprof(jtype) = 'plchltot'
447            clproffiles(jtype,:) = cn_plchltotfbfiles
448         ENDIF
449         IF (ln_pchltot) THEN
450            jtype = jtype + 1
451            cobstypesprof(jtype) = 'pchltot'
452            clproffiles(jtype,:) = cn_pchltotfbfiles
453         ENDIF
454         IF (ln_pno3) THEN
455            jtype = jtype + 1
456            cobstypesprof(jtype) = 'pno3'
457            clproffiles(jtype,:) = cn_pno3fbfiles
458         ENDIF
459         IF (ln_psi4) THEN
460            jtype = jtype + 1
461            cobstypesprof(jtype) = 'psi4'
462            clproffiles(jtype,:) = cn_psi4fbfiles
463         ENDIF
464         IF (ln_ppo4) THEN
465            jtype = jtype + 1
466            cobstypesprof(jtype) = 'ppo4'
467            clproffiles(jtype,:) = cn_ppo4fbfiles
468         ENDIF
469         IF (ln_pdic) THEN
470            jtype = jtype + 1
471            cobstypesprof(jtype) = 'pdic'
472            clproffiles(jtype,:) = cn_pdicfbfiles
473         ENDIF
474         IF (ln_palk) THEN
475            jtype = jtype + 1
476            cobstypesprof(jtype) = 'palk'
477            clproffiles(jtype,:) = cn_palkfbfiles
478         ENDIF
479         IF (ln_pph) THEN
480            jtype = jtype + 1
481            cobstypesprof(jtype) = 'pph'
482            clproffiles(jtype,:) = cn_pphfbfiles
483         ENDIF
484         IF (ln_po2) THEN
485            jtype = jtype + 1
486            cobstypesprof(jtype) = 'po2'
487            clproffiles(jtype,:) = cn_po2fbfiles
488         ENDIF
489
490         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles )
491
492      ENDIF
493
494      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
495      IF ( nsurftypes > 0 ) THEN
496
497         ALLOCATE( cobstypessurf(nsurftypes) )
498         ALLOCATE( ifilessurf(nsurftypes) )
499         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
500         ALLOCATE(n2dintsurf(nsurftypes))
501         ALLOCATE(ravglamscl(nsurftypes))
502         ALLOCATE(ravgphiscl(nsurftypes))
503         ALLOCATE(lfpindegs(nsurftypes))
504         ALLOCATE(llnightav(nsurftypes))
505
506         jtype = 0
507         IF (ln_sla) THEN
508            jtype = jtype + 1
509            cobstypessurf(jtype) = 'sla'
510            clsurffiles(jtype,:) = cn_slafbfiles
511         ENDIF
512         IF (ln_sst) THEN
513            jtype = jtype + 1
514            cobstypessurf(jtype) = 'sst'
515            clsurffiles(jtype,:) = cn_sstfbfiles
516         ENDIF
517         IF (ln_sic) THEN
518            jtype = jtype + 1
519            cobstypessurf(jtype) = 'sic'
520            clsurffiles(jtype,:) = cn_sicfbfiles
521         ENDIF
522         IF (ln_sss) THEN
523            jtype = jtype + 1
524            cobstypessurf(jtype) = 'sss'
525            clsurffiles(jtype,:) = cn_sssfbfiles
526         ENDIF
527         IF (ln_slchltot) THEN
528            jtype = jtype + 1
529            cobstypessurf(jtype) = 'slchltot'
530            clsurffiles(jtype,:) = cn_slchltotfbfiles
531         ENDIF
532         IF (ln_slchldia) THEN
533            jtype = jtype + 1
534            cobstypessurf(jtype) = 'slchldia'
535            clsurffiles(jtype,:) = cn_slchldiafbfiles
536         ENDIF
537         IF (ln_slchlnon) THEN
538            jtype = jtype + 1
539            cobstypessurf(jtype) = 'slchlnon'
540            clsurffiles(jtype,:) = cn_slchlnonfbfiles
541         ENDIF
542         IF (ln_slchldin) THEN
543            jtype = jtype + 1
544            cobstypessurf(jtype) = 'slchldin'
545            clsurffiles(jtype,:) = cn_slchldinfbfiles
546         ENDIF
547         IF (ln_slchlmic) THEN
548            jtype = jtype + 1
549            cobstypessurf(jtype) = 'slchlmic'
550            clsurffiles(jtype,:) = cn_slchlmicfbfiles
551         ENDIF
552         IF (ln_slchlnan) THEN
553            jtype = jtype + 1
554            cobstypessurf(jtype) = 'slchlnan'
555            clsurffiles(jtype,:) = cn_slchlnanfbfiles
556         ENDIF
557         IF (ln_slchlpic) THEN
558            jtype = jtype + 1
559            cobstypessurf(jtype) = 'slchlpic'
560            clsurffiles(jtype,:) = cn_slchlpicfbfiles
561         ENDIF
562         IF (ln_schltot) THEN
563            jtype = jtype + 1
564            cobstypessurf(jtype) = 'schltot'
565            clsurffiles(jtype,:) = cn_schltotfbfiles
566         ENDIF
567         IF (ln_sspm) THEN
568            jtype = jtype + 1
569            cobstypessurf(jtype) = 'sspm'
570            clsurffiles(jtype,:) = cn_sspmfbfiles
571         ENDIF
572         IF (ln_sfco2) THEN
573            jtype = jtype + 1
574            cobstypessurf(jtype) = 'sfco2'
575            clsurffiles(jtype,:) = cn_sfco2fbfiles
576         ENDIF
577         IF (ln_spco2) THEN
578            jtype = jtype + 1
579            cobstypessurf(jtype) = 'spco2'
580            clsurffiles(jtype,:) = cn_spco2fbfiles
581         ENDIF
582
583         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles )
584
585         DO jtype = 1, nsurftypes
586
587            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
588               IF ( nn_2dint_sla == -1 ) THEN
589                  n2dint_type  = nn_2dint_default
590               ELSE
591                  n2dint_type  = nn_2dint_sla
592               ENDIF
593               rtype_avglamscl = rn_sla_avglamscl
594               rtype_avgphiscl = rn_sla_avgphiscl
595               ltype_fp_indegs = ln_sla_fp_indegs
596               ltype_night     = .FALSE.
597            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
598               IF ( nn_2dint_sst == -1 ) THEN
599                  n2dint_type  = nn_2dint_default
600               ELSE
601                  n2dint_type  = nn_2dint_sst
602               ENDIF
603               rtype_avglamscl = rn_sst_avglamscl
604               rtype_avgphiscl = rn_sst_avgphiscl
605               ltype_fp_indegs = ln_sst_fp_indegs
606               ltype_night     = ln_sstnight
607            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
608               IF ( nn_2dint_sic == -1 ) THEN
609                  n2dint_type  = nn_2dint_default
610               ELSE
611                  n2dint_type  = nn_2dint_sic
612               ENDIF
613               rtype_avglamscl = rn_sic_avglamscl
614               rtype_avgphiscl = rn_sic_avgphiscl
615               ltype_fp_indegs = ln_sic_fp_indegs
616               ltype_night     = .FALSE.
617            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
618               IF ( nn_2dint_sss == -1 ) THEN
619                  n2dint_type  = nn_2dint_default
620               ELSE
621                  n2dint_type  = nn_2dint_sss
622               ENDIF
623               rtype_avglamscl = rn_sss_avglamscl
624               rtype_avgphiscl = rn_sss_avgphiscl
625               ltype_fp_indegs = ln_sss_fp_indegs
626               ltype_night     = .FALSE.
627            ELSE
628               n2dint_type     = nn_2dint_default
629               rtype_avglamscl = rn_default_avglamscl
630               rtype_avgphiscl = rn_default_avgphiscl
631               ltype_fp_indegs = ln_default_fp_indegs
632               ltype_night     = .FALSE.
633            ENDIF
634           
635            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), &
636               &                    nn_2dint_default, n2dint_type,                 &
637               &                    rtype_avglamscl, rtype_avgphiscl,              &
638               &                    ltype_fp_indegs, ltype_night,                  &
639               &                    n2dintsurf, ravglamscl, ravgphiscl,            &
640               &                    lfpindegs, llnightav )
641
642         END DO
643
644      ENDIF
645
646      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
647
648
649      !-----------------------------------------------------------------------
650      ! Obs operator parameter checking and initialisations
651      !-----------------------------------------------------------------------
652
653      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
654         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
655         RETURN
656      ENDIF
657
658      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
659         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
660            &                    ' is not available')
661      ENDIF
662
663      IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN
664         CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', &
665            &                    ' is not available')
666      ENDIF
667
668      CALL obs_typ_init
669
670      CALL mppmap_init
671
672      CALL obs_grid_setup( )
673
674      !-----------------------------------------------------------------------
675      ! Depending on switches read the various observation types
676      !-----------------------------------------------------------------------
677
678      IF ( nproftypes > 0 ) THEN
679
680         ALLOCATE(profdata(nproftypes))
681         ALLOCATE(profdataqc(nproftypes))
682         ALLOCATE(nvarsprof(nproftypes))
683         ALLOCATE(nextrprof(nproftypes))
684
685         DO jtype = 1, nproftypes
686
687            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
688               nvarsprof(jtype) = 2
689               nextrprof(jtype) = 1
690               ALLOCATE(llvar(nvarsprof(jtype)))
691               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
692               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
693               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
694               llvar(1)       = ln_t3d
695               llvar(2)       = ln_s3d
696               zglam(:,:,1)   = glamt(:,:)
697               zglam(:,:,2)   = glamt(:,:)
698               zgphi(:,:,1)   = gphit(:,:)
699               zgphi(:,:,2)   = gphit(:,:)
700               zmask(:,:,:,1) = tmask(:,:,:)
701               zmask(:,:,:,2) = tmask(:,:,:)
702            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
703               nvarsprof(jtype) = 2
704               nextrprof(jtype) = 2
705               ALLOCATE(llvar(nvarsprof(jtype)))
706               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
707               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
708               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
709               llvar(1)       = ln_vel3d
710               llvar(2)       = ln_vel3d
711               zglam(:,:,1)   = glamu(:,:)
712               zglam(:,:,2)   = glamv(:,:)
713               zgphi(:,:,1)   = gphiu(:,:)
714               zgphi(:,:,2)   = gphiv(:,:)
715               zmask(:,:,:,1) = umask(:,:,:)
716               zmask(:,:,:,2) = vmask(:,:,:)
717            ELSE
718               nvarsprof(jtype) = 1
719               nextrprof(jtype) = 0
720               ALLOCATE(llvar(nvarsprof(jtype)))
721               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
722               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
723               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
724               llvar(1)       = .TRUE.
725               zglam(:,:,1)   = glamt(:,:)
726               zgphi(:,:,1)   = gphit(:,:)
727               zmask(:,:,:,1) = tmask(:,:,:)
728            ENDIF
729
730            !Read in profile or profile obs types
731            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
732               &               clproffiles(jtype,1:ifilesprof(jtype)), &
733               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
734               &               rn_dobsini, rn_dobsend, llvar, &
735               &               ln_ignmis, ln_s_at_t, .FALSE., &
736               &               kdailyavtypes = nn_profdavtypes )
737
738            DO jvar = 1, nvarsprof(jtype)
739               CALL obs_prof_staend( profdata(jtype), jvar )
740            END DO
741
742            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
743               &               llvar, &
744               &               jpi, jpj, jpk, &
745               &               zmask, zglam, zgphi,  &
746               &               ln_nea, ln_bound_reject, &
747               &               kdailyavtypes = nn_profdavtypes )
748           
749            DEALLOCATE( llvar )
750            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zglam )
751            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zgphi )
752            CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
753
754         END DO
755
756         DEALLOCATE( ifilesprof, clproffiles )
757
758      ENDIF
759
760      IF ( nsurftypes > 0 ) THEN
761
762         ALLOCATE(surfdata(nsurftypes))
763         ALLOCATE(surfdataqc(nsurftypes))
764         ALLOCATE(nvarssurf(nsurftypes))
765         ALLOCATE(nextrsurf(nsurftypes))
766
767         DO jtype = 1, nsurftypes
768
769            nvarssurf(jtype) = 1
770            nextrsurf(jtype) = 0
771            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
772
773            !Read in surface obs types
774            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
775               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
776               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
777               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
778
779            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
780
781            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
782               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
783               IF ( ln_altbias ) &
784                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
785            ENDIF
786
787            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
788               jnumsstbias = 0
789               DO jfile = 1, jpmaxnfiles
790                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
791                     &  jnumsstbias = jnumsstbias + 1
792               END DO
793               IF ( jnumsstbias == 0 ) THEN
794                  CALL ctl_stop("ln_sstbias set but no bias files to read in")   
795               ENDIF
796
797               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
798                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
799
800            ENDIF
801
802         END DO
803
804         DEALLOCATE( ifilessurf, clsurffiles )
805
806      ENDIF
807
808   END SUBROUTINE dia_obs_init
809
810   SUBROUTINE dia_obs( kstp )
811      !!----------------------------------------------------------------------
812      !!                    ***  ROUTINE dia_obs  ***
813      !!         
814      !! ** Purpose : Call the observation operators on each time step
815      !!
816      !! ** Method  : Call the observation operators on each time step to
817      !!              compute the model equivalent of the following data:
818      !!               - Profile data, currently T/S or U/V
819      !!               - Surface data, currently SST, SLA or sea-ice concentration.
820      !!
821      !! ** Action  :
822      !!
823      !! History :
824      !!        !  06-03  (K. Mogensen) Original code
825      !!        !  06-05  (K. Mogensen) Reformatted
826      !!        !  06-10  (A. Weaver) Cleaning
827      !!        !  07-03  (K. Mogensen) General handling of profiles
828      !!        !  07-04  (G. Smith) Generalized surface operators
829      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
830      !!        !  15-08  (M. Martin) Combined surface/profile routines.
831      !!----------------------------------------------------------------------
832      !! * Modules used
833      USE phycst, ONLY : &         ! Physical constants
834         & rday
835      USE oce, ONLY : &            ! Ocean dynamics and tracers variables
836         & tsn,       &
837         & un,        &
838         & vn,        &
839         & sshn
840#if defined  key_lim3
841      USE ice, ONLY : &            ! LIM3 Ice model variables
842         & frld
843#endif
844#if defined key_lim2
845      USE ice_2, ONLY : &          ! LIM2 Ice model variables
846         & frld
847#endif
848#if defined key_cice
849      USE sbc_oce, ONLY : fr_i     ! ice fraction
850#endif
851#if defined key_hadocc
852      USE trc, ONLY :  &           ! HadOCC variables
853         & trn, &
854         & HADOCC_CHL, &
855         & HADOCC_FCO2, &
856         & HADOCC_PCO2, &
857         & HADOCC_FILL_FLT
858      USE par_hadocc
859#elif defined key_medusa && defined key_foam_medusa
860      USE trc, ONLY :  &           ! MEDUSA variables
861         & trn
862      USE par_medusa
863#if defined key_roam
864      USE sms_medusa, ONLY: &
865         & f2_pco2w, &
866         & f2_fco2w, &
867         & f3_pH
868#endif
869#elif defined key_fabm
870      USE fabm
871      USE par_fabm
872#endif
873#if defined key_spm
874      USE par_spm, ONLY: &         ! ERSEM/SPM sediments
875         & jp_spm
876      USE trc, ONLY :  &
877         & trn
878#endif
879
880      IMPLICIT NONE
881
882      !! * Arguments
883      INTEGER, INTENT(IN) :: kstp  ! Current timestep
884      !! * Local declarations
885      INTEGER :: idaystp           ! Number of timesteps per day
886      INTEGER :: jtype             ! Data loop variable
887      INTEGER :: jvar              ! Variable number
888      INTEGER :: ji, jj            ! Loop counters
889      REAL(wp) :: tiny             ! small number
890      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
891         & zprofvar                ! Model values for variables in a prof ob
892      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
893         & zprofmask               ! Mask associated with zprofvar
894      REAL(wp), POINTER, DIMENSION(:,:) :: &
895         & zsurfvar, &             ! Model values equivalent to surface ob.
896         & zsurfmask               ! Mask associated with surface variable
897      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
898         & zglam,    &             ! Model longitudes for prof variables
899         & zgphi                   ! Model latitudes for prof variables
900      LOGICAL :: llog10            ! Perform log10 transform of variable
901
902
903      IF(lwp) THEN
904         WRITE(numout,*)
905         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
906         WRITE(numout,*) '~~~~~~~'
907         CALL FLUSH(numout)
908      ENDIF
909
910      idaystp = NINT( rday / rdt )
911
912      !-----------------------------------------------------------------------
913      ! Call the profile and surface observation operators
914      !-----------------------------------------------------------------------
915
916      IF ( nproftypes > 0 ) THEN
917
918         DO jtype = 1, nproftypes
919
920            ! Allocate local work arrays
921            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
922            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
923            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
924            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
925           
926            ! Defaults which might change
927            DO jvar = 1, profdataqc(jtype)%nvar
928               zprofmask(:,:,:,jvar) = tmask(:,:,:)
929               zglam(:,:,jvar)       = glamt(:,:)
930               zgphi(:,:,jvar)       = gphit(:,:)
931            END DO
932
933            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
934
935            CASE('prof')
936               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
937               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
938
939            CASE('vel')
940               zprofvar(:,:,:,1) = un(:,:,:)
941               zprofvar(:,:,:,2) = vn(:,:,:)
942               zprofmask(:,:,:,1) = umask(:,:,:)
943               zprofmask(:,:,:,2) = vmask(:,:,:)
944               zglam(:,:,1) = glamu(:,:)
945               zglam(:,:,2) = glamv(:,:)
946               zgphi(:,:,1) = gphiu(:,:)
947               zgphi(:,:,2) = gphiv(:,:)
948
949            CASE('plchltot')
950#if defined key_hadocc
951               ! Chlorophyll from HadOCC
952               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
953#elif defined key_medusa && defined key_foam_medusa
954               ! Add non-diatom and diatom chlorophyll from MEDUSA
955               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
956#elif defined key_fabm
957               ! Add all chlorophyll groups from ERSEM
958               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &
959                  &                trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)
960#else
961               CALL ctl_stop( ' Trying to run plchltot observation operator', &
962                  &           ' but no biogeochemical model appears to have been defined' )
963#endif
964               ! Take the log10 where we can, otherwise exclude
965               tiny = 1.0e-20
966               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )
967                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:))
968               ELSEWHERE
969                  zprofvar(:,:,:,:)  = obfillflt
970                  zprofmask(:,:,:,:) = 0
971               END WHERE
972
973            CASE('pchltot')
974#if defined key_hadocc
975               ! Chlorophyll from HadOCC
976               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
977#elif defined key_medusa && defined key_foam_medusa
978               ! Add non-diatom and diatom chlorophyll from MEDUSA
979               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
980#elif defined key_fabm
981               ! Add all chlorophyll groups from ERSEM
982               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_chl1) + trn(:,:,:,jp_fabm_chl2) + &
983                  &                trn(:,:,:,jp_fabm_chl3) + trn(:,:,:,jp_fabm_chl4)
984#else
985               CALL ctl_stop( ' Trying to run pchltot observation operator', &
986                  &           ' but no biogeochemical model appears to have been defined' )
987#endif
988
989            CASE('pno3')
990#if defined key_hadocc
991               ! Dissolved inorganic nitrogen from HadOCC
992               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)
993#elif defined key_medusa && defined key_foam_medusa
994               ! Dissolved inorganic nitrogen from MEDUSA
995               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
996#elif defined key_fabm
997               ! Nitrate from ERSEM
998               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n3n)
999#else
1000               CALL ctl_stop( ' Trying to run pno3 observation operator', &
1001                  &           ' but no biogeochemical model appears to have been defined' )
1002#endif
1003
1004            CASE('psi4')
1005#if defined key_hadocc
1006               CALL ctl_stop( ' Trying to run psi4 observation operator', &
1007                  &           ' but HadOCC does not simulate silicate' )
1008#elif defined key_medusa && defined key_foam_medusa
1009               ! Silicate from MEDUSA
1010               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
1011#elif defined key_fabm
1012               ! Silicate from ERSEM
1013               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n5s)
1014#else
1015               CALL ctl_stop( ' Trying to run psi4 observation operator', &
1016                  &           ' but no biogeochemical model appears to have been defined' )
1017#endif
1018
1019            CASE('ppo4')
1020#if defined key_hadocc
1021               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
1022                  &           ' but HadOCC does not simulate phosphate' )
1023#elif defined key_medusa && defined key_foam_medusa
1024               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
1025                  &           ' but MEDUSA does not simulate phosphate' )
1026#elif defined key_fabm
1027               ! Phosphate from ERSEM
1028               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_n1p)
1029#else
1030               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
1031                  &           ' but no biogeochemical model appears to have been defined' )
1032#endif
1033
1034            CASE('pdic')
1035#if defined key_hadocc
1036               ! Dissolved inorganic carbon from HadOCC
1037               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)
1038#elif defined key_medusa && defined key_foam_medusa
1039               ! Dissolved inorganic carbon from MEDUSA
1040               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
1041#elif defined key_fabm
1042               ! Dissolved inorganic carbon from ERSEM
1043               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3c)
1044#else
1045               CALL ctl_stop( ' Trying to run pdic observation operator', &
1046                  &           ' but no biogeochemical model appears to have been defined' )
1047#endif
1048
1049            CASE('palk')
1050#if defined key_hadocc
1051               ! Alkalinity from HadOCC
1052               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)
1053#elif defined key_medusa && defined key_foam_medusa
1054               ! Alkalinity from MEDUSA
1055               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
1056#elif defined key_fabm
1057               ! Alkalinity from ERSEM
1058               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3a)
1059#else
1060               CALL ctl_stop( ' Trying to run palk observation operator', &
1061                  &           ' but no biogeochemical model appears to have been defined' )
1062#endif
1063
1064            CASE('pph')
1065#if defined key_hadocc
1066               CALL ctl_stop( ' Trying to run pph observation operator', &
1067                  &           ' but HadOCC has no pH diagnostic defined' )
1068#elif defined key_medusa && defined key_foam_medusa
1069               ! pH from MEDUSA
1070               zprofvar(:,:,:,1) = f3_pH(:,:,:)
1071#elif defined key_fabm
1072               ! pH from ERSEM
1073               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o3ph)
1074#else
1075               CALL ctl_stop( ' Trying to run pph observation operator', &
1076                  &           ' but no biogeochemical model appears to have been defined' )
1077#endif
1078
1079            CASE('po2')
1080#if defined key_hadocc
1081               CALL ctl_stop( ' Trying to run po2 observation operator', &
1082                  &           ' but HadOCC does not simulate oxygen' )
1083#elif defined key_medusa && defined key_foam_medusa
1084               ! Oxygen from MEDUSA
1085               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
1086#elif defined key_fabm
1087               ! Oxygen from ERSEM
1088               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_o2o)
1089#else
1090               CALL ctl_stop( ' Trying to run po2 observation operator', &
1091                  &           ' but no biogeochemical model appears to have been defined' )
1092#endif
1093
1094            CASE DEFAULT
1095               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
1096
1097            END SELECT
1098
1099            DO jvar = 1, profdataqc(jtype)%nvar
1100               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
1101                  &               nit000, idaystp, jvar,                   &
1102                  &               zprofvar(:,:,:,jvar),                    &
1103                  &               fsdept(:,:,:), fsdepw(:,:,:),            & 
1104                  &               zprofmask(:,:,:,jvar),                   &
1105                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
1106                  &               nn_1dint, nn_2dint_default,              &
1107                  &               kdailyavtypes = nn_profdavtypes )
1108            END DO
1109
1110            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
1111            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
1112            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
1113            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
1114
1115         END DO
1116
1117      ENDIF
1118
1119      IF ( nsurftypes > 0 ) THEN
1120
1121         !Allocate local work arrays
1122         CALL wrk_alloc( jpi, jpj, zsurfvar )
1123         CALL wrk_alloc( jpi, jpj, zsurfmask )
1124
1125         DO jtype = 1, nsurftypes
1126
1127            !Defaults which might be changed
1128            zsurfmask(:,:) = tmask(:,:,1)
1129            llog10 = .FALSE.
1130
1131            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
1132            CASE('sst')
1133               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
1134            CASE('sla')
1135               zsurfvar(:,:) = sshn(:,:)
1136            CASE('sss')
1137               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
1138            CASE('sic')
1139               IF ( kstp == 0 ) THEN
1140                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
1141                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
1142                        &           'time-step but some obs are valid then.' )
1143                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
1144                        &           ' sea-ice obs will be missed'
1145                  ENDIF
1146                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
1147                     &                        surfdataqc(jtype)%nsstp(1)
1148                  CYCLE
1149               ELSE
1150#if defined key_cice
1151                  zsurfvar(:,:) = fr_i(:,:)
1152#elif defined key_lim2 || defined key_lim3
1153                  zsurfvar(:,:) = 1._wp - frld(:,:)
1154#else
1155               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
1156                  &           ' but no sea-ice model appears to have been defined' )
1157#endif
1158               ENDIF
1159
1160            CASE('slchltot')
1161#if defined key_hadocc
1162               ! Surface chlorophyll from HadOCC
1163               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
1164#elif defined key_medusa && defined key_foam_medusa
1165               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
1166               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
1167#elif defined key_fabm
1168               ! Add all surface chlorophyll groups from ERSEM
1169               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl2) + &
1170                  &            trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,jp_fabm_chl4)
1171#else
1172               CALL ctl_stop( ' Trying to run slchltot observation operator', &
1173                  &           ' but no biogeochemical model appears to have been defined' )
1174#endif
1175               llog10 = .TRUE.
1176
1177            CASE('slchldia')
1178#if defined key_hadocc
1179               CALL ctl_stop( ' Trying to run slchldia observation operator', &
1180                  &           ' but HadOCC does not explicitly simulate diatoms' )
1181#elif defined key_medusa && defined key_foam_medusa
1182               ! Diatom surface chlorophyll from MEDUSA
1183               zsurfvar(:,:) = trn(:,:,1,jpchd)
1184#elif defined key_fabm
1185               ! Diatom surface chlorophyll from ERSEM
1186               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1)
1187#else
1188               CALL ctl_stop( ' Trying to run slchldia observation operator', &
1189                  &           ' but no biogeochemical model appears to have been defined' )
1190#endif
1191               llog10 = .TRUE.
1192
1193            CASE('slchlnon')
1194#if defined key_hadocc
1195               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
1196                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
1197#elif defined key_medusa && defined key_foam_medusa
1198               ! Non-diatom surface chlorophyll from MEDUSA
1199               zsurfvar(:,:) = trn(:,:,1,jpchn)
1200#elif defined key_fabm
1201               ! Add all non-diatom surface chlorophyll groups from ERSEM
1202               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl2) + &
1203                  &            trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,jp_fabm_chl4)
1204#else
1205               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
1206                  &           ' but no biogeochemical model appears to have been defined' )
1207#endif
1208               llog10 = .TRUE.
1209
1210            CASE('slchldin')
1211#if defined key_hadocc
1212               CALL ctl_stop( ' Trying to run slchldin observation operator', &
1213                  &           ' but HadOCC does not explicitly simulate dinoflagellates' )
1214#elif defined key_medusa && defined key_foam_medusa
1215               CALL ctl_stop( ' Trying to run slchldin observation operator', &
1216                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' )
1217#elif defined key_fabm
1218               ! Dinoflagellate surface chlorophyll from ERSEM
1219               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl4)
1220#else
1221               CALL ctl_stop( ' Trying to run slchldin observation operator', &
1222                  &           ' but no biogeochemical model appears to have been defined' )
1223#endif
1224               llog10 = .TRUE.
1225
1226            CASE('slchlmic')
1227#if defined key_hadocc
1228               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
1229                  &           ' but HadOCC does not explicitly simulate microphytoplankton' )
1230#elif defined key_medusa && defined key_foam_medusa
1231               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
1232                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' )
1233#elif defined key_fabm
1234               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM
1235               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl4)
1236#else
1237               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
1238                  &           ' but no biogeochemical model appears to have been defined' )
1239#endif
1240               llog10 = .TRUE.
1241
1242            CASE('slchlnan')
1243#if defined key_hadocc
1244               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
1245                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' )
1246#elif defined key_medusa && defined key_foam_medusa
1247               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
1248                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' )
1249#elif defined key_fabm
1250               ! Nanophytoplankton surface chlorophyll from ERSEM
1251               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl2)
1252#else
1253               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
1254                  &           ' but no biogeochemical model appears to have been defined' )
1255#endif
1256               llog10 = .TRUE.
1257
1258            CASE('slchlpic')
1259#if defined key_hadocc
1260               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
1261                  &           ' but HadOCC does not explicitly simulate picophytoplankton' )
1262#elif defined key_medusa && defined key_foam_medusa
1263               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
1264                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' )
1265#elif defined key_fabm
1266               ! Picophytoplankton surface chlorophyll from ERSEM
1267               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl3)
1268#else
1269               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
1270                  &           ' but no biogeochemical model appears to have been defined' )
1271#endif
1272               llog10 = .TRUE.
1273
1274            CASE('schltot')
1275#if defined key_hadocc
1276               ! Surface chlorophyll from HadOCC
1277               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
1278#elif defined key_medusa && defined key_foam_medusa
1279               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
1280               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
1281#elif defined key_fabm
1282               ! Add all surface chlorophyll groups from ERSEM
1283               zsurfvar(:,:) = trn(:,:,1,jp_fabm_chl1) + trn(:,:,1,jp_fabm_chl2) + &
1284                  &            trn(:,:,1,jp_fabm_chl3) + trn(:,:,1,jp_fabm_chl4)
1285#else
1286               CALL ctl_stop( ' Trying to run schltot observation operator', &
1287                  &           ' but no biogeochemical model appears to have been defined' )
1288#endif
1289
1290            CASE('sspm')
1291#if defined key_spm
1292               zsurfvar(:,:) = 0.0
1293               DO jn = 1, jp_spm
1294                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
1295               END DO
1296#else
1297               CALL ctl_stop( ' Trying to run sspm observation operator', &
1298                  &           ' but no spm model appears to have been defined' )
1299#endif
1300
1301            CASE('sfco2')
1302#if defined key_hadocc
1303               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
1304               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
1305                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
1306                  zsurfvar(:,:) = obfillflt
1307                  zsurfmask(:,:) = 0
1308                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
1309                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
1310               ENDIF
1311#elif defined key_medusa && defined key_foam_medusa && defined key_roam
1312               zsurfvar(:,:) = f2_fco2w(:,:)
1313#elif defined key_fabm
1314               ! First, get pCO2 from FABM
1315               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
1316               zsurfvar(:,:) = pco2_3d(:,:,1)
1317               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
1318               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
1319               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
1320               ! and
1321               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
1322               ! Marine Chemistry, 2: 203-215.
1323               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
1324               ! not explicitly included - atmospheric pressure is not necessarily available so this is
1325               ! the best assumption.
1326               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
1327               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
1328               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
1329               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
1330               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + &
1331                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
1332                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
1333                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
1334                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
1335                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
1336#else
1337               CALL ctl_stop( ' Trying to run sfco2 observation operator', &
1338                  &           ' but no biogeochemical model appears to have been defined' )
1339#endif
1340
1341            CASE('spco2')
1342#if defined key_hadocc
1343               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
1344               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
1345                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
1346                  zsurfvar(:,:) = obfillflt
1347                  zsurfmask(:,:) = 0
1348                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
1349                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
1350               ENDIF
1351#elif defined key_medusa && defined key_foam_medusa && defined key_roam
1352               zsurfvar(:,:) = f2_pco2w(:,:)
1353#elif defined key_fabm
1354               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
1355               zsurfvar(:,:) = pco2_3d(:,:,1)
1356#else
1357               CALL ctl_stop( ' Trying to run spco2 observation operator', &
1358                  &           ' but no biogeochemical model appears to have been defined' )
1359#endif
1360
1361            CASE DEFAULT
1362
1363               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
1364
1365            END SELECT
1366           
1367            IF ( llog10 ) THEN
1368               ! Take the log10 where we can, otherwise exclude
1369               tiny = 1.0e-20
1370               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
1371                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
1372               ELSEWHERE
1373                  zsurfvar(:,:)  = obfillflt
1374                  zsurfmask(:,:) = 0
1375               END WHERE
1376            ENDIF
1377
1378            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
1379               &               nit000, idaystp, zsurfvar, zsurfmask,    &
1380               &               n2dintsurf(jtype), llnightav(jtype),     &
1381               &               ravglamscl(jtype), ravgphiscl(jtype),     &
1382               &               lfpindegs(jtype) )
1383
1384         END DO
1385
1386         CALL wrk_dealloc( jpi, jpj, zsurfvar )
1387         CALL wrk_dealloc( jpi, jpj, zsurfmask )
1388
1389      ENDIF
1390
1391   END SUBROUTINE dia_obs
1392
1393   SUBROUTINE dia_obs_wri
1394      !!----------------------------------------------------------------------
1395      !!                    ***  ROUTINE dia_obs_wri  ***
1396      !!         
1397      !! ** Purpose : Call observation diagnostic output routines
1398      !!
1399      !! ** Method  : Call observation diagnostic output routines
1400      !!
1401      !! ** Action  :
1402      !!
1403      !! History :
1404      !!        !  06-03  (K. Mogensen) Original code
1405      !!        !  06-05  (K. Mogensen) Reformatted
1406      !!        !  06-10  (A. Weaver) Cleaning
1407      !!        !  07-03  (K. Mogensen) General handling of profiles
1408      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
1409      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
1410      !!----------------------------------------------------------------------
1411      !! * Modules used
1412      USE obs_rot_vel          ! Rotation of velocities
1413
1414      IMPLICIT NONE
1415
1416      !! * Local declarations
1417      INTEGER :: jtype                    ! Data set loop variable
1418      INTEGER :: jo, jvar, jk
1419      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
1420         & zu, &
1421         & zv
1422
1423      !-----------------------------------------------------------------------
1424      ! Depending on switches call various observation output routines
1425      !-----------------------------------------------------------------------
1426
1427      IF ( nproftypes > 0 ) THEN
1428
1429         DO jtype = 1, nproftypes
1430
1431            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
1432
1433               ! For velocity data, rotate the model velocities to N/S, E/W
1434               ! using the compressed data structure.
1435               ALLOCATE( &
1436                  & zu(profdataqc(jtype)%nvprot(1)), &
1437                  & zv(profdataqc(jtype)%nvprot(2))  &
1438                  & )
1439
1440               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
1441
1442               DO jo = 1, profdataqc(jtype)%nprof
1443                  DO jvar = 1, 2
1444                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
1445
1446                        IF ( jvar == 1 ) THEN
1447                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
1448                        ELSE
1449                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
1450                        ENDIF
1451
1452                     END DO
1453                  END DO
1454               END DO
1455
1456               DEALLOCATE( zu )
1457               DEALLOCATE( zv )
1458
1459            END IF
1460
1461            CALL obs_prof_decompress( profdataqc(jtype), &
1462               &                      profdata(jtype), .TRUE., numout )
1463
1464            CALL obs_wri_prof( profdata(jtype) )
1465
1466         END DO
1467
1468      ENDIF
1469
1470      IF ( nsurftypes > 0 ) THEN
1471
1472         DO jtype = 1, nsurftypes
1473
1474            CALL obs_surf_decompress( surfdataqc(jtype), &
1475               &                      surfdata(jtype), .TRUE., numout )
1476
1477            CALL obs_wri_surf( surfdata(jtype) )
1478
1479         END DO
1480
1481      ENDIF
1482
1483   END SUBROUTINE dia_obs_wri
1484
1485   SUBROUTINE dia_obs_dealloc
1486      IMPLICIT NONE
1487      !!----------------------------------------------------------------------
1488      !!                    *** ROUTINE dia_obs_dealloc ***
1489      !!
1490      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1491      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1492      !!
1493      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1494      !!
1495      !!  ** Action :
1496      !!
1497      !!----------------------------------------------------------------------
1498      ! obs_grid deallocation
1499      CALL obs_grid_deallocate
1500
1501      ! diaobs deallocation
1502      IF ( nproftypes > 0 ) &
1503         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
1504
1505      IF ( nsurftypes > 0 ) &
1506         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
1507         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
1508
1509   END SUBROUTINE dia_obs_dealloc
1510
1511   SUBROUTINE ini_date( ddobsini )
1512      !!----------------------------------------------------------------------
1513      !!                    ***  ROUTINE ini_date  ***
1514      !!
1515      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1516      !!
1517      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
1518      !!
1519      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
1520      !!
1521      !! History :
1522      !!        !  06-03  (K. Mogensen)  Original code
1523      !!        !  06-05  (K. Mogensen)  Reformatted
1524      !!        !  06-10  (A. Weaver) Cleaning
1525      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1526      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1527      !!----------------------------------------------------------------------
1528      USE phycst, ONLY : &            ! Physical constants
1529         & rday
1530      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1531         & rdt
1532
1533      IMPLICIT NONE
1534
1535      !! * Arguments
1536      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
1537
1538      !! * Local declarations
1539      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1540      INTEGER :: imon
1541      INTEGER :: iday
1542      INTEGER :: ihou
1543      INTEGER :: imin
1544      INTEGER :: imday       ! Number of days in month.
1545      INTEGER, DIMENSION(12) :: &
1546         &       imonth_len  ! Length in days of the months of the current year
1547      REAL(wp) :: zdayfrc    ! Fraction of day
1548
1549      !----------------------------------------------------------------------
1550      ! Initial date initialization (year, month, day, hour, minute)
1551      ! (This assumes that the initial date is for 00z))
1552      !----------------------------------------------------------------------
1553      iyea =   ndate0 / 10000
1554      imon = ( ndate0 - iyea * 10000 ) / 100
1555      iday =   ndate0 - iyea * 10000 - imon * 100
1556      ihou = 0
1557      imin = 0
1558
1559      !----------------------------------------------------------------------
1560      ! Compute number of days + number of hours + min since initial time
1561      !----------------------------------------------------------------------
1562      iday = iday + ( nit000 -1 ) * rdt / rday
1563      zdayfrc = ( nit000 -1 ) * rdt / rday
1564      zdayfrc = zdayfrc - aint(zdayfrc)
1565      ihou = int( zdayfrc * 24 )
1566      imin = int( (zdayfrc * 24 - ihou) * 60 )
1567
1568      !-----------------------------------------------------------------------
1569      ! Convert number of days (iday) into a real date
1570      !----------------------------------------------------------------------
1571
1572      CALL calc_month_len( iyea, imonth_len )
1573
1574      DO WHILE ( iday > imonth_len(imon) )
1575         iday = iday - imonth_len(imon)
1576         imon = imon + 1 
1577         IF ( imon > 12 ) THEN
1578            imon = 1
1579            iyea = iyea + 1
1580            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1581         ENDIF
1582      END DO
1583
1584      !----------------------------------------------------------------------
1585      ! Convert it into YYYYMMDD.HHMMSS format.
1586      !----------------------------------------------------------------------
1587      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1588         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1589
1590
1591   END SUBROUTINE ini_date
1592
1593   SUBROUTINE fin_date( ddobsfin )
1594      !!----------------------------------------------------------------------
1595      !!                    ***  ROUTINE fin_date  ***
1596      !!
1597      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1598      !!
1599      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
1600      !!
1601      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
1602      !!
1603      !! History :
1604      !!        !  06-03  (K. Mogensen)  Original code
1605      !!        !  06-05  (K. Mogensen)  Reformatted
1606      !!        !  06-10  (A. Weaver) Cleaning
1607      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1608      !!----------------------------------------------------------------------
1609      USE phycst, ONLY : &            ! Physical constants
1610         & rday
1611      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1612         & rdt
1613
1614      IMPLICIT NONE
1615
1616      !! * Arguments
1617      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1618
1619      !! * Local declarations
1620      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1621      INTEGER :: imon
1622      INTEGER :: iday
1623      INTEGER :: ihou
1624      INTEGER :: imin
1625      INTEGER :: imday       ! Number of days in month.
1626      INTEGER, DIMENSION(12) :: &
1627         &       imonth_len  ! Length in days of the months of the current year
1628      REAL(wp) :: zdayfrc    ! Fraction of day
1629
1630      !-----------------------------------------------------------------------
1631      ! Initial date initialization (year, month, day, hour, minute)
1632      ! (This assumes that the initial date is for 00z)
1633      !-----------------------------------------------------------------------
1634      iyea =   ndate0 / 10000
1635      imon = ( ndate0 - iyea * 10000 ) / 100
1636      iday =   ndate0 - iyea * 10000 - imon * 100
1637      ihou = 0
1638      imin = 0
1639
1640      !-----------------------------------------------------------------------
1641      ! Compute number of days + number of hours + min since initial time
1642      !-----------------------------------------------------------------------
1643      iday    = iday +  nitend  * rdt / rday
1644      zdayfrc =  nitend  * rdt / rday
1645      zdayfrc = zdayfrc - AINT( zdayfrc )
1646      ihou    = INT( zdayfrc * 24 )
1647      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1648
1649      !-----------------------------------------------------------------------
1650      ! Convert number of days (iday) into a real date
1651      !----------------------------------------------------------------------
1652
1653      CALL calc_month_len( iyea, imonth_len )
1654
1655      DO WHILE ( iday > imonth_len(imon) )
1656         iday = iday - imonth_len(imon)
1657         imon = imon + 1 
1658         IF ( imon > 12 ) THEN
1659            imon = 1
1660            iyea = iyea + 1
1661            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1662         ENDIF
1663      END DO
1664
1665      !-----------------------------------------------------------------------
1666      ! Convert it into YYYYMMDD.HHMMSS format
1667      !-----------------------------------------------------------------------
1668      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1669         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1670
1671    END SUBROUTINE fin_date
1672
1673    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
1674
1675       INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
1676       INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
1677       INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
1678          &                   ifiles      ! Out number of files for each type
1679       CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
1680          &                   cobstypes   ! List of obs types
1681       CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
1682          &                   cfiles      ! List of files for all types
1683
1684       !Local variables
1685       INTEGER :: jfile
1686       INTEGER :: jtype
1687
1688       DO jtype = 1, ntypes
1689
1690          ifiles(jtype) = 0
1691          DO jfile = 1, jpmaxnfiles
1692             IF ( trim(cfiles(jtype,jfile)) /= '' ) &
1693                       ifiles(jtype) = ifiles(jtype) + 1
1694          END DO
1695
1696          IF ( ifiles(jtype) == 0 ) THEN
1697               CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   &
1698                  &           ' set to true but no files available to read' )
1699          ENDIF
1700
1701          IF(lwp) THEN   
1702             WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
1703             DO jfile = 1, ifiles(jtype)
1704                WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
1705             END DO
1706          ENDIF
1707
1708       END DO
1709
1710    END SUBROUTINE obs_settypefiles
1711
1712    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
1713               &                  n2dint_default, n2dint_type,        &
1714               &                  ravglamscl_type, ravgphiscl_type,   &
1715               &                  lfp_indegs_type, lavnight_type,     &
1716               &                  n2dint, ravglamscl, ravgphiscl,     &
1717               &                  lfpindegs, lavnight )
1718
1719       INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
1720       INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
1721       INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
1722       INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
1723       REAL(wp), INTENT(IN) :: &
1724          &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
1725          &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
1726       LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
1727       LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
1728       CHARACTER(len=8), INTENT(IN) :: ctypein 
1729
1730       INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1731          &                    n2dint 
1732       REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1733          &                    ravglamscl, ravgphiscl
1734       LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1735          &                    lfpindegs, lavnight
1736
1737       lavnight(jtype) = lavnight_type
1738
1739       IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN
1740          n2dint(jtype) = n2dint_type
1741       ELSE IF ( n2dint_type == -1 ) THEN
1742          n2dint(jtype) = n2dint_default
1743       ELSE
1744          CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
1745            &                    ' is not available')
1746       ENDIF
1747
1748       ! For averaging observation footprints set options for size of footprint
1749       IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
1750          IF ( ravglamscl_type > 0._wp ) THEN
1751             ravglamscl(jtype) = ravglamscl_type
1752          ELSE
1753             CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1754                            'scale (ravglamscl) for observation type '//TRIM(ctypein) )     
1755          ENDIF
1756
1757          IF ( ravgphiscl_type > 0._wp ) THEN
1758             ravgphiscl(jtype) = ravgphiscl_type
1759          ELSE
1760             CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1761                            'scale (ravgphiscl) for observation type '//TRIM(ctypein) )     
1762          ENDIF
1763
1764          lfpindegs(jtype) = lfp_indegs_type 
1765
1766       ENDIF
1767
1768       ! Write out info
1769       IF(lwp) THEN
1770          IF ( n2dint(jtype) <= 4 ) THEN
1771             WRITE(numout,*) '             '//TRIM(ctypein)// &
1772                &            ' model counterparts will be interpolated horizontally'
1773          ELSE IF ( n2dint(jtype) <= 6 ) THEN
1774             WRITE(numout,*) '             '//TRIM(ctypein)// &
1775                &            ' model counterparts will be averaged horizontally'
1776             WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
1777             WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
1778             IF ( lfpindegs(jtype) ) THEN
1779                 WRITE(numout,*) '             '//'    (in degrees)'
1780             ELSE
1781                 WRITE(numout,*) '             '//'    (in metres)'
1782             ENDIF
1783          ENDIF
1784       ENDIF
1785
1786    END SUBROUTINE obs_setinterpopts
1787
1788END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.