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 @ 9229

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

Add options for phytoplankton carbon observations.

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