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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 12140

Last change on this file since 12140 was 12140, checked in by kingr, 4 years ago
Added option to remove tides from SLA bkg by taking average over 24h50m.
File size: 95.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: