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

Last change on this file since 10729 was 10729, checked in by dford, 5 years ago
Minor bug fixes to using observation operator with FABM-ERSEM. See https://code.metoffice.gov.uk/trac/utils/ticket/174.
File size: 87.6 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
[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 :: &
[9306]	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
[9306]	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
[9306]	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
[9306]	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]	112	CONTAINS
	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
[9306]	146	INTEGER :: n2dint_type ! Local version of nn_2dint*
[7992]	147
	148	CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
[9306]	149	& cn_profbfiles, & ! T/S profile input filenames
	150	& cn_sstfbfiles, & ! Sea surface temperature input filenames
	151	& cn_slafbfiles, & ! Sea level anomaly input filenames
	152	& cn_sicfbfiles, & ! Seaice concentration input filenames
	153	& cn_velfbfiles, & ! Velocity profile input filenames
	154	& cn_sssfbfiles, & ! Sea surface salinity input filenames
	155	& cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames
	156	& cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames
	157	& cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames
	158	& cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames
	159	& cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames
	160	& cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames
	161	& cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames
	162	& cn_schltotfbfiles, & ! Surface total chlorophyll input filenames
	163	& cn_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
	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
	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
[7992]	178	& cn_sstbiasfiles ! SST bias input filenames
	179
	180	CHARACTER(LEN=128) :: &
	181	& cn_altbiasfile ! Altimeter bias input filename
	182
	183
	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
[9306]	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
	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
	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
	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
[7992]	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
	220
	221	REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS
	222	REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS
	223
[9306]	224	REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl
	225	REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl
	226
[7992]	227	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
	228	& clproffiles, & ! Profile filenames
	229	& clsurffiles ! Surface filenames
	230
[9306]	231	LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read
	232	LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
	233	LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables)
[7992]	234
	235	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
[9306]	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, &
[9306]	245	& ln_slchltot, ln_slchldia, ln_slchlnon, &
	246	& ln_slchldin, ln_slchlmic, ln_slchlnan, &
	247	& ln_slchlpic, ln_schltot, &
	248	& ln_slphytot, ln_slphydia, ln_slphynon, &
	249	& ln_sspm, ln_sfco2, ln_spco2, &
	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, &
[9306]	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, &
[9306]	262	& cn_slchltotfbfiles, cn_slchldiafbfiles, &
	263	& cn_slchlnonfbfiles, cn_slchldinfbfiles, &
	264	& cn_slchlmicfbfiles, cn_slchlnanfbfiles, &
	265	& cn_slchlpicfbfiles, cn_schltotfbfiles, &
	266	& cn_slphytotfbfiles, cn_slphydiafbfiles, &
	267	& cn_slphynonfbfiles, cn_sspmfbfiles, &
	268	& cn_sfco2fbfiles, cn_spco2fbfiles, &
	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, &
[9306]	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, &
[9306]	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
[9306]	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(:) = ''
	306	cn_slphytotfbfiles(:) = ''
	307	cn_slphydiafbfiles(:) = ''
	308	cn_slphynonfbfiles(:) = ''
	309	cn_sspmfbfiles(:) = ''
	310	cn_sfco2fbfiles(:) = ''
	311	cn_spco2fbfiles(:) = ''
	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(:) = ''
	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)
	330	901 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 )
	334	902 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
[9306]	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
	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
	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
	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
[9306]	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
[9306]	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 /) )
	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, &
	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
[9306]	451	cobstypesprof(jtype) = 'prof'
	452	clproffiles(jtype,:) = cn_profbfiles
[2128]	453	ENDIF
[7992]	454	IF (ln_vel3d) THEN
	455	jtype = jtype + 1
[9306]	456	cobstypesprof(jtype) = 'vel'
	457	clproffiles(jtype,:) = cn_velfbfiles
[2128]	458	ENDIF
[9306]	459	IF (ln_plchltot) THEN
	460	jtype = jtype + 1
	461	cobstypesprof(jtype) = 'plchltot'
	462	clproffiles(jtype,:) = cn_plchltotfbfiles
	463	ENDIF
	464	IF (ln_pchltot) THEN
	465	jtype = jtype + 1
	466	cobstypesprof(jtype) = 'pchltot'
	467	clproffiles(jtype,:) = cn_pchltotfbfiles
	468	ENDIF
	469	IF (ln_pno3) THEN
	470	jtype = jtype + 1
	471	cobstypesprof(jtype) = 'pno3'
	472	clproffiles(jtype,:) = cn_pno3fbfiles
	473	ENDIF
	474	IF (ln_psi4) THEN
	475	jtype = jtype + 1
	476	cobstypesprof(jtype) = 'psi4'
	477	clproffiles(jtype,:) = cn_psi4fbfiles
	478	ENDIF
	479	IF (ln_ppo4) THEN
	480	jtype = jtype + 1
	481	cobstypesprof(jtype) = 'ppo4'
	482	clproffiles(jtype,:) = cn_ppo4fbfiles
	483	ENDIF
	484	IF (ln_pdic) THEN
	485	jtype = jtype + 1
	486	cobstypesprof(jtype) = 'pdic'
	487	clproffiles(jtype,:) = cn_pdicfbfiles
	488	ENDIF
	489	IF (ln_palk) THEN
	490	jtype = jtype + 1
	491	cobstypesprof(jtype) = 'palk'
	492	clproffiles(jtype,:) = cn_palkfbfiles
	493	ENDIF
	494	IF (ln_pph) THEN
	495	jtype = jtype + 1
	496	cobstypesprof(jtype) = 'pph'
	497	clproffiles(jtype,:) = cn_pphfbfiles
	498	ENDIF
	499	IF (ln_po2) THEN
	500	jtype = jtype + 1
	501	cobstypesprof(jtype) = 'po2'
	502	clproffiles(jtype,:) = cn_po2fbfiles
	503	ENDIF
[2128]	504
[9306]	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
[9306]	524	cobstypessurf(jtype) = 'sla'
	525	clsurffiles(jtype,:) = cn_slafbfiles
[2128]	526	ENDIF
[7992]	527	IF (ln_sst) THEN
	528	jtype = jtype + 1
[9306]	529	cobstypessurf(jtype) = 'sst'
	530	clsurffiles(jtype,:) = cn_sstfbfiles
[2128]	531	ENDIF
[7992]	532	IF (ln_sic) THEN
	533	jtype = jtype + 1
[9306]	534	cobstypessurf(jtype) = 'sic'
	535	clsurffiles(jtype,:) = cn_sicfbfiles
[2128]	536	ENDIF
[7992]	537	IF (ln_sss) THEN
	538	jtype = jtype + 1
[9306]	539	cobstypessurf(jtype) = 'sss'
	540	clsurffiles(jtype,:) = cn_sssfbfiles
[2128]	541	ENDIF
[9306]	542	IF (ln_slchltot) THEN
[7992]	543	jtype = jtype + 1
[9306]	544	cobstypessurf(jtype) = 'slchltot'
	545	clsurffiles(jtype,:) = cn_slchltotfbfiles
[7992]	546	ENDIF
[9306]	547	IF (ln_slchldia) THEN
[7992]	548	jtype = jtype + 1
[9306]	549	cobstypessurf(jtype) = 'slchldia'
	550	clsurffiles(jtype,:) = cn_slchldiafbfiles
[7992]	551	ENDIF
[9306]	552	IF (ln_slchlnon) THEN
[7992]	553	jtype = jtype + 1
[9306]	554	cobstypessurf(jtype) = 'slchlnon'
	555	clsurffiles(jtype,:) = cn_slchlnonfbfiles
[2128]	556	ENDIF
[9306]	557	IF (ln_slchldin) THEN
[7992]	558	jtype = jtype + 1
[9306]	559	cobstypessurf(jtype) = 'slchldin'
	560	clsurffiles(jtype,:) = cn_slchldinfbfiles
[2128]	561	ENDIF
[9306]	562	IF (ln_slchlmic) THEN
	563	jtype = jtype + 1
	564	cobstypessurf(jtype) = 'slchlmic'
	565	clsurffiles(jtype,:) = cn_slchlmicfbfiles
	566	ENDIF
	567	IF (ln_slchlnan) THEN
	568	jtype = jtype + 1
	569	cobstypessurf(jtype) = 'slchlnan'
	570	clsurffiles(jtype,:) = cn_slchlnanfbfiles
	571	ENDIF
	572	IF (ln_slchlpic) THEN
	573	jtype = jtype + 1
	574	cobstypessurf(jtype) = 'slchlpic'
	575	clsurffiles(jtype,:) = cn_slchlpicfbfiles
	576	ENDIF
	577	IF (ln_schltot) THEN
	578	jtype = jtype + 1
	579	cobstypessurf(jtype) = 'schltot'
	580	clsurffiles(jtype,:) = cn_schltotfbfiles
	581	ENDIF
	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
	597	IF (ln_sspm) THEN
	598	jtype = jtype + 1
	599	cobstypessurf(jtype) = 'sspm'
	600	clsurffiles(jtype,:) = cn_sspmfbfiles
	601	ENDIF
	602	IF (ln_sfco2) THEN
	603	jtype = jtype + 1
	604	cobstypessurf(jtype) = 'sfco2'
	605	clsurffiles(jtype,:) = cn_sfco2fbfiles
	606	ENDIF
	607	IF (ln_spco2) THEN
	608	jtype = jtype + 1
	609	cobstypessurf(jtype) = 'spco2'
	610	clsurffiles(jtype,:) = cn_spco2fbfiles
	611	ENDIF
[2128]	612
[9306]	613	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles )
	614
	615	DO jtype = 1, nsurftypes
	616
	617	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
	618	IF ( nn_2dint_sla == -1 ) THEN
	619	n2dint_type = nn_2dint_default
	620	ELSE
	621	n2dint_type = nn_2dint_sla
	622	ENDIF
	623	ztype_avglamscl = rn_sla_avglamscl
	624	ztype_avgphiscl = rn_sla_avgphiscl
	625	ltype_fp_indegs = ln_sla_fp_indegs
	626	ltype_night = .FALSE.
	627	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
	628	IF ( nn_2dint_sst == -1 ) THEN
	629	n2dint_type = nn_2dint_default
	630	ELSE
	631	n2dint_type = nn_2dint_sst
	632	ENDIF
	633	ztype_avglamscl = rn_sst_avglamscl
	634	ztype_avgphiscl = rn_sst_avgphiscl
	635	ltype_fp_indegs = ln_sst_fp_indegs
	636	ltype_night = ln_sstnight
	637	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
	638	IF ( nn_2dint_sic == -1 ) THEN
	639	n2dint_type = nn_2dint_default
	640	ELSE
	641	n2dint_type = nn_2dint_sic
	642	ENDIF
	643	ztype_avglamscl = rn_sic_avglamscl
	644	ztype_avgphiscl = rn_sic_avgphiscl
	645	ltype_fp_indegs = ln_sic_fp_indegs
	646	ltype_night = .FALSE.
	647	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
	648	IF ( nn_2dint_sss == -1 ) THEN
	649	n2dint_type = nn_2dint_default
	650	ELSE
	651	n2dint_type = nn_2dint_sss
	652	ENDIF
	653	ztype_avglamscl = rn_sss_avglamscl
	654	ztype_avgphiscl = rn_sss_avgphiscl
	655	ltype_fp_indegs = ln_sss_fp_indegs
	656	ltype_night = .FALSE.
	657	ELSE
	658	n2dint_type = nn_2dint_default
	659	ztype_avglamscl = rn_default_avglamscl
	660	ztype_avgphiscl = rn_default_avgphiscl
	661	ltype_fp_indegs = ln_default_fp_indegs
	662	ltype_night = .FALSE.
	663	ENDIF
	664
	665	CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), &
	666	& nn_2dint_default, n2dint_type, &
	667	& ztype_avglamscl, ztype_avgphiscl, &
	668	& ltype_fp_indegs, ltype_night, &
	669	& n2dintsurf, ravglamscl, ravgphiscl, &
	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
[9306]	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
[9306]	718	nvarsprof(jtype) = 2
[7992]	719	nextrprof(jtype) = 1
[9306]	720	ALLOCATE(llvar(nvarsprof(jtype)))
	721	CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam )
	722	CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi )
	723	CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
	724	llvar(1) = 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(:,:,:)
	732	ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
	733	nvarsprof(jtype) = 2
[7992]	734	nextrprof(jtype) = 2
[9306]	735	ALLOCATE(llvar(nvarsprof(jtype)))
	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(:,:,:)
	747	ELSE
	748	nvarsprof(jtype) = 1
	749	nextrprof(jtype) = 0
	750	ALLOCATE(llvar(nvarsprof(jtype)))
	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, &
[9306]	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), &
[9306]	773	& llvar, &
[7992]	774	& jpi, jpj, jpk, &
[9306]	775	& zmask, zglam, zgphi, &
[7992]	776	& ln_nea, ln_bound_reject, &
	777	& kdailyavtypes = nn_profdavtypes )
[9306]	778
	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
[10388]	864	#if defined key_fabm
	865	& rt0, &
	866	#endif
[7992]	867	& rday
	868	USE oce, ONLY : & ! Ocean dynamics and tracers variables
	869	& tsn, &
	870	& un, &
	871	& vn, &
[2128]	872	& sshn
[3294]	873	#if defined key_lim3
[7992]	874	USE ice, ONLY : & ! LIM3 Ice model variables
[2128]	875	& frld
	876	#endif
[3294]	877	#if defined key_lim2
[7992]	878	USE ice_2, ONLY : & ! LIM2 Ice model variables
[3294]	879	& frld
[2715]	880	#endif
[7992]	881	#if defined key_cice
	882	USE sbc_oce, ONLY : fr_i ! ice fraction
	883	#endif
[10388]	884	#if defined key_top
	885	USE trc, ONLY : & ! Biogeochemical state variables
	886	& trn
	887	#endif
[7992]	888	#if defined key_hadocc
[10388]	889	USE par_hadocc ! HadOCC parameters
	890	USE trc, ONLY : &
[7992]	891	& HADOCC_CHL, &
	892	& HADOCC_FCO2, &
	893	& HADOCC_PCO2, &
	894	& HADOCC_FILL_FLT
[9306]	895	USE had_bgc_const, ONLY: c2n_p
[10168]	896	#elif defined key_medusa
[10388]	897	USE par_medusa ! MEDUSA parameters
[9306]	898	USE sms_medusa, ONLY: &
	899	& xthetapn, &
	900	& xthetapd
[8653]	901	#if defined key_roam
	902	USE sms_medusa, ONLY: &
	903	& f2_pco2w, &
[9306]	904	& f2_fco2w, &
	905	& f3_pH
[8653]	906	#endif
[7992]	907	#elif defined key_fabm
[10388]	908	USE par_fabm ! FABM parameters
	909	USE fabm, ONLY: &
[10729]	910	& fabm_get_interior_diagnostic_data
[7992]	911	#endif
	912	#if defined key_spm
[10388]	913	USE par_spm, ONLY: & ! Sediment parameters
[7992]	914	& jp_spm
	915	#endif
	916
[2128]	917	IMPLICIT NONE
	918
	919	!! * Arguments
[7992]	920	INTEGER, INTENT(IN) :: kstp ! Current timestep
[2128]	921	!! * Local declarations
[7992]	922	INTEGER :: idaystp ! Number of timesteps per day
	923	INTEGER :: jtype ! Data loop variable
	924	INTEGER :: jvar ! Variable number
[9306]	925	INTEGER :: ji, jj, jk ! Loop counters
	926	REAL(wp) :: tiny ! small number
	927	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
	928	& zprofvar ! Model values for variables in a prof ob
	929	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
	930	& zprofmask ! Mask associated with zprofvar
[7992]	931	REAL(wp), POINTER, DIMENSION(:,:) :: &
	932	& zsurfvar, & ! Model values equivalent to surface ob.
	933	& zsurfmask ! Mask associated with surface variable
[9306]	934	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
	935	& zglam, & ! Model longitudes for prof variables
	936	& zgphi ! Model latitudes for prof variables
	937	LOGICAL :: llog10 ! Perform log10 transform of variable
[10388]	938	#if defined key_fabm
	939	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
	940	& pco2_3d ! 3D pCO2 from FABM
	941	#endif
[2715]	942
[2128]	943	IF(lwp) THEN
	944	WRITE(numout,*)
	945	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
	946	WRITE(numout,*) '~~~~~~~'
[7992]	947	CALL FLUSH(numout)
[2128]	948	ENDIF
	949
	950	idaystp = NINT( rday / rdt )
	951
	952	!-----------------------------------------------------------------------
[7992]	953	! Call the profile and surface observation operators
[2128]	954	!-----------------------------------------------------------------------
	955
[7992]	956	IF ( nproftypes > 0 ) THEN
[2128]	957
[7992]	958	DO jtype = 1, nproftypes
[2128]	959
[9306]	960	! Allocate local work arrays
	961	CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar )
	962	CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
	963	CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam )
	964	CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi )
	965
	966	! Defaults which might change
	967	DO jvar = 1, profdataqc(jtype)%nvar
	968	zprofmask(:,:,:,jvar) = tmask(:,:,:)
	969	zglam(:,:,jvar) = glamt(:,:)
	970	zgphi(:,:,jvar) = gphit(:,:)
	971	END DO
	972
[7992]	973	SELECT CASE ( TRIM(cobstypesprof(jtype)) )
[9306]	974
[7992]	975	CASE('prof')
[9306]	976	zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
	977	zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
	978
[7992]	979	CASE('vel')
[9306]	980	zprofvar(:,:,:,1) = un(:,:,:)
	981	zprofvar(:,:,:,2) = vn(:,:,:)
	982	zprofmask(:,:,:,1) = umask(:,:,:)
	983	zprofmask(:,:,:,2) = vmask(:,:,:)
	984	zglam(:,:,1) = glamu(:,:)
	985	zglam(:,:,2) = glamv(:,:)
	986	zgphi(:,:,1) = gphiu(:,:)
	987	zgphi(:,:,2) = gphiv(:,:)
	988
	989	CASE('plchltot')
	990	#if defined key_hadocc
	991	! Chlorophyll from HadOCC
	992	zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
[10168]	993	#elif defined key_medusa
[9306]	994	! Add non-diatom and diatom chlorophyll from MEDUSA
	995	zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
	996	#elif defined key_fabm
	997	! Add all chlorophyll groups from ERSEM
[10729]	998	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
	999	& trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
[9306]	1000	#else
	1001	CALL ctl_stop( ' Trying to run plchltot observation operator', &
	1002	& ' but no biogeochemical model appears to have been defined' )
	1003	#endif
	1004	! Take the log10 where we can, otherwise exclude
	1005	tiny = 1.0e-20
	1006	WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )
	1007	zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:))
	1008	ELSEWHERE
	1009	zprofvar(:,:,:,:) = obfillflt
	1010	zprofmask(:,:,:,:) = 0
	1011	END WHERE
	1012	! Mask out model below any excluded values,
	1013	! to avoid interpolation issues
	1014	DO jvar = 1, profdataqc(jtype)%nvar
	1015	DO jj = 1, jpj
	1016	DO ji = 1, jpi
	1017	depth_loop: DO jk = 1, jpk
	1018	IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN
	1019	zprofmask(ji,jj,jk:jpk,jvar) = 0
	1020	EXIT depth_loop
	1021	ENDIF
	1022	END DO depth_loop
	1023	END DO
	1024	END DO
	1025	END DO
	1026
	1027	CASE('pchltot')
	1028	#if defined key_hadocc
	1029	! Chlorophyll from HadOCC
	1030	zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
[10168]	1031	#elif defined key_medusa
[9306]	1032	! Add non-diatom and diatom chlorophyll from MEDUSA
	1033	zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
	1034	#elif defined key_fabm
	1035	! Add all chlorophyll groups from ERSEM
[10729]	1036	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
	1037	& trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
[9306]	1038	#else
	1039	CALL ctl_stop( ' Trying to run pchltot observation operator', &
	1040	& ' but no biogeochemical model appears to have been defined' )
	1041	#endif
	1042
	1043	CASE('pno3')
	1044	#if defined key_hadocc
	1045	! Dissolved inorganic nitrogen from HadOCC
	1046	zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)
[10168]	1047	#elif defined key_medusa
[9306]	1048	! Dissolved inorganic nitrogen from MEDUSA
	1049	zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
	1050	#elif defined key_fabm
	1051	! Nitrate from ERSEM
[10729]	1052	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)
[9306]	1053	#else
	1054	CALL ctl_stop( ' Trying to run pno3 observation operator', &
	1055	& ' but no biogeochemical model appears to have been defined' )
	1056	#endif
	1057
	1058	CASE('psi4')
	1059	#if defined key_hadocc
	1060	CALL ctl_stop( ' Trying to run psi4 observation operator', &
	1061	& ' but HadOCC does not simulate silicate' )
[10168]	1062	#elif defined key_medusa
[9306]	1063	! Silicate from MEDUSA
	1064	zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
	1065	#elif defined key_fabm
	1066	! Silicate from ERSEM
[10729]	1067	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)
[9306]	1068	#else
	1069	CALL ctl_stop( ' Trying to run psi4 observation operator', &
	1070	& ' but no biogeochemical model appears to have been defined' )
	1071	#endif
	1072
	1073	CASE('ppo4')
	1074	#if defined key_hadocc
	1075	CALL ctl_stop( ' Trying to run ppo4 observation operator', &
	1076	& ' but HadOCC does not simulate phosphate' )
[10168]	1077	#elif defined key_medusa
[9306]	1078	CALL ctl_stop( ' Trying to run ppo4 observation operator', &
	1079	& ' but MEDUSA does not simulate phosphate' )
	1080	#elif defined key_fabm
	1081	! Phosphate from ERSEM
[10729]	1082	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)
[9306]	1083	#else
	1084	CALL ctl_stop( ' Trying to run ppo4 observation operator', &
	1085	& ' but no biogeochemical model appears to have been defined' )
	1086	#endif
	1087
	1088	CASE('pdic')
	1089	#if defined key_hadocc
	1090	! Dissolved inorganic carbon from HadOCC
	1091	zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)
[10168]	1092	#elif defined key_medusa
[9306]	1093	! Dissolved inorganic carbon from MEDUSA
	1094	zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
	1095	#elif defined key_fabm
	1096	! Dissolved inorganic carbon from ERSEM
[10729]	1097	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
[9306]	1098	#else
	1099	CALL ctl_stop( ' Trying to run pdic observation operator', &
	1100	& ' but no biogeochemical model appears to have been defined' )
	1101	#endif
	1102
	1103	CASE('palk')
	1104	#if defined key_hadocc
	1105	! Alkalinity from HadOCC
	1106	zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)
[10168]	1107	#elif defined key_medusa
[9306]	1108	! Alkalinity from MEDUSA
	1109	zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
	1110	#elif defined key_fabm
	1111	! Alkalinity from ERSEM
[10729]	1112	zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
[9306]	1113	#else
	1114	CALL ctl_stop( ' Trying to run palk observation operator', &
	1115	& ' but no biogeochemical model appears to have been defined' )
	1116	#endif
	1117
	1118	CASE('pph')
	1119	#if defined key_hadocc
	1120	CALL ctl_stop( ' Trying to run pph observation operator', &
	1121	& ' but HadOCC has no pH diagnostic defined' )
[10168]	1122	#elif defined key_medusa && defined key_roam
[9306]	1123	! pH from MEDUSA
	1124	zprofvar(:,:,:,1) = f3_pH(:,:,:)
	1125	#elif defined key_fabm
	1126	! pH from ERSEM
[10729]	1127	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph)
[9306]	1128	#else
	1129	CALL ctl_stop( ' Trying to run pph observation operator', &
	1130	& ' but no biogeochemical model appears to have been defined' )
	1131	#endif
	1132
	1133	CASE('po2')
	1134	#if defined key_hadocc
	1135	CALL ctl_stop( ' Trying to run po2 observation operator', &
	1136	& ' but HadOCC does not simulate oxygen' )
[10168]	1137	#elif defined key_medusa
[9306]	1138	! Oxygen from MEDUSA
	1139	zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
	1140	#elif defined key_fabm
	1141	! Oxygen from ERSEM
[10729]	1142	zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o)
[9306]	1143	#else
	1144	CALL ctl_stop( ' Trying to run po2 observation operator', &
	1145	& ' but no biogeochemical model appears to have been defined' )
	1146	#endif
	1147
[7992]	1148	CASE DEFAULT
	1149	CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
[9306]	1150
[7992]	1151	END SELECT
	1152
[9306]	1153	DO jvar = 1, profdataqc(jtype)%nvar
	1154	CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, &
	1155	& nit000, idaystp, jvar, &
	1156	& zprofvar(:,:,:,jvar), &
	1157	& fsdept(:,:,:), fsdepw(:,:,:), &
	1158	& zprofmask(:,:,:,jvar), &
	1159	& zglam(:,:,jvar), zgphi(:,:,jvar), &
	1160	& nn_1dint, nn_2dint_default, &
	1161	& kdailyavtypes = nn_profdavtypes )
	1162	END DO
[7992]	1163
[9306]	1164	CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar )
	1165	CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
	1166	CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam )
	1167	CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi )
	1168
[2128]	1169	END DO
	1170
	1171	ENDIF
	1172
[7992]	1173	IF ( nsurftypes > 0 ) THEN
	1174
[9306]	1175	!Allocate local work arrays
	1176	CALL wrk_alloc( jpi, jpj, zsurfvar )
	1177	CALL wrk_alloc( jpi, jpj, zsurfmask )
[10729]	1178	#if defined key_fabm
	1179	CALL wrk_alloc( jpi, jpj, jpk, pco2_3d )
	1180	#endif
[9306]	1181
[7992]	1182	DO jtype = 1, nsurftypes
	1183
	1184	!Defaults which might be changed
	1185	zsurfmask(:,:) = tmask(:,:,1)
[9306]	1186	llog10 = .FALSE.
[7992]	1187
	1188	SELECT CASE ( TRIM(cobstypessurf(jtype)) )
	1189	CASE('sst')
	1190	zsurfvar(:,:) = tsn(:,:,1,jp_tem)
	1191	CASE('sla')
	1192	zsurfvar(:,:) = sshn(:,:)
	1193	CASE('sss')
	1194	zsurfvar(:,:) = tsn(:,:,1,jp_sal)
	1195	CASE('sic')
	1196	IF ( kstp == 0 ) THEN
	1197	IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
	1198	CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
	1199	& 'time-step but some obs are valid then.' )
	1200	WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
	1201	& ' sea-ice obs will be missed'
	1202	ENDIF
	1203	surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
	1204	& surfdataqc(jtype)%nsstp(1)
	1205	CYCLE
	1206	ELSE
	1207	#if defined key_cice
	1208	zsurfvar(:,:) = fr_i(:,:)
	1209	#elif defined key_lim2 \|\| defined key_lim3
	1210	zsurfvar(:,:) = 1._wp - frld(:,:)
	1211	#else
	1212	CALL ctl_stop( ' Trying to run sea-ice observation operator', &
	1213	& ' but no sea-ice model appears to have been defined' )
[2128]	1214	#endif
[7992]	1215	ENDIF
[2128]	1216
[9306]	1217	CASE('slchltot')
[7992]	1218	#if defined key_hadocc
[9306]	1219	! Surface chlorophyll from HadOCC
	1220	zsurfvar(:,:) = HADOCC_CHL(:,:,1)
[10168]	1221	#elif defined key_medusa
[8653]	1222	! Add non-diatom and diatom surface chlorophyll from MEDUSA
	1223	zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
[7992]	1224	#elif defined key_fabm
[9306]	1225	! Add all surface chlorophyll groups from ERSEM
[10729]	1226	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
	1227	& trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
[7992]	1228	#else
[9306]	1229	CALL ctl_stop( ' Trying to run slchltot observation operator', &
[7992]	1230	& ' but no biogeochemical model appears to have been defined' )
	1231	#endif
[9306]	1232	llog10 = .TRUE.
	1233
	1234	CASE('slchldia')
	1235	#if defined key_hadocc
	1236	CALL ctl_stop( ' Trying to run slchldia observation operator', &
	1237	& ' but HadOCC does not explicitly simulate diatoms' )
[10168]	1238	#elif defined key_medusa
[9306]	1239	! Diatom surface chlorophyll from MEDUSA
	1240	zsurfvar(:,:) = trn(:,:,1,jpchd)
	1241	#elif defined key_fabm
	1242	! Diatom surface chlorophyll from ERSEM
[10729]	1243	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)
[9306]	1244	#else
	1245	CALL ctl_stop( ' Trying to run slchldia observation operator', &
	1246	& ' but no biogeochemical model appears to have been defined' )
	1247	#endif
	1248	llog10 = .TRUE.
	1249
	1250	CASE('slchlnon')
	1251	#if defined key_hadocc
	1252	CALL ctl_stop( ' Trying to run slchlnon observation operator', &
	1253	& ' but HadOCC does not explicitly simulate non-diatoms' )
[10168]	1254	#elif defined key_medusa
[9306]	1255	! Non-diatom surface chlorophyll from MEDUSA
	1256	zsurfvar(:,:) = trn(:,:,1,jpchn)
	1257	#elif defined key_fabm
	1258	! Add all non-diatom surface chlorophyll groups from ERSEM
[10729]	1259	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
	1260	& trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
[9306]	1261	#else
	1262	CALL ctl_stop( ' Trying to run slchlnon observation operator', &
	1263	& ' but no biogeochemical model appears to have been defined' )
	1264	#endif
	1265	llog10 = .TRUE.
	1266
	1267	CASE('slchldin')
	1268	#if defined key_hadocc
	1269	CALL ctl_stop( ' Trying to run slchldin observation operator', &
	1270	& ' but HadOCC does not explicitly simulate dinoflagellates' )
[10168]	1271	#elif defined key_medusa
[9306]	1272	CALL ctl_stop( ' Trying to run slchldin observation operator', &
	1273	& ' but MEDUSA does not explicitly simulate dinoflagellates' )
	1274	#elif defined key_fabm
	1275	! Dinoflagellate surface chlorophyll from ERSEM
[10729]	1276	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
[9306]	1277	#else
	1278	CALL ctl_stop( ' Trying to run slchldin observation operator', &
	1279	& ' but no biogeochemical model appears to have been defined' )
	1280	#endif
	1281	llog10 = .TRUE.
	1282
	1283	CASE('slchlmic')
	1284	#if defined key_hadocc
	1285	CALL ctl_stop( ' Trying to run slchlmic observation operator', &
	1286	& ' but HadOCC does not explicitly simulate microphytoplankton' )
[10168]	1287	#elif defined key_medusa
[9306]	1288	CALL ctl_stop( ' Trying to run slchlmic observation operator', &
	1289	& ' but MEDUSA does not explicitly simulate microphytoplankton' )
	1290	#elif defined key_fabm
	1291	! Add diatom and dinoflagellate surface chlorophyll from ERSEM
[10729]	1292	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
[9306]	1293	#else
	1294	CALL ctl_stop( ' Trying to run slchlmic observation operator', &
	1295	& ' but no biogeochemical model appears to have been defined' )
	1296	#endif
	1297	llog10 = .TRUE.
	1298
	1299	CASE('slchlnan')
	1300	#if defined key_hadocc
	1301	CALL ctl_stop( ' Trying to run slchlnan observation operator', &
	1302	& ' but HadOCC does not explicitly simulate nanophytoplankton' )
[10168]	1303	#elif defined key_medusa
[9306]	1304	CALL ctl_stop( ' Trying to run slchlnan observation operator', &
	1305	& ' but MEDUSA does not explicitly simulate nanophytoplankton' )
	1306	#elif defined key_fabm
	1307	! Nanophytoplankton surface chlorophyll from ERSEM
[10729]	1308	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)
[9306]	1309	#else
	1310	CALL ctl_stop( ' Trying to run slchlnan observation operator', &
	1311	& ' but no biogeochemical model appears to have been defined' )
	1312	#endif
	1313	llog10 = .TRUE.
	1314
	1315	CASE('slchlpic')
	1316	#if defined key_hadocc
	1317	CALL ctl_stop( ' Trying to run slchlpic observation operator', &
	1318	& ' but HadOCC does not explicitly simulate picophytoplankton' )
[10168]	1319	#elif defined key_medusa
[9306]	1320	CALL ctl_stop( ' Trying to run slchlpic observation operator', &
	1321	& ' but MEDUSA does not explicitly simulate picophytoplankton' )
	1322	#elif defined key_fabm
	1323	! Picophytoplankton surface chlorophyll from ERSEM
[10729]	1324	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)
[9306]	1325	#else
	1326	CALL ctl_stop( ' Trying to run slchlpic observation operator', &
	1327	& ' but no biogeochemical model appears to have been defined' )
	1328	#endif
	1329	llog10 = .TRUE.
	1330
	1331	CASE('schltot')
	1332	#if defined key_hadocc
	1333	! Surface chlorophyll from HadOCC
	1334	zsurfvar(:,:) = HADOCC_CHL(:,:,1)
[10168]	1335	#elif defined key_medusa
[9306]	1336	! Add non-diatom and diatom surface chlorophyll from MEDUSA
	1337	zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
	1338	#elif defined key_fabm
	1339	! Add all surface chlorophyll groups from ERSEM
[10729]	1340	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
	1341	& trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
[9306]	1342	#else
	1343	CALL ctl_stop( ' Trying to run schltot observation operator', &
	1344	& ' but no biogeochemical model appears to have been defined' )
	1345	#endif
	1346
	1347	CASE('slphytot')
	1348	#if defined key_hadocc
	1349	! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio
	1350	zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
[10168]	1351	#elif defined key_medusa
[9306]	1352	! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA
	1353	! multiplied by C:N ratio for each
	1354	zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
	1355	#elif defined key_fabm
	1356	! Add all surface phytoplankton carbon groups from ERSEM
[10729]	1357	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
	1358	& trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
[9306]	1359	#else
	1360	CALL ctl_stop( ' Trying to run slphytot observation operator', &
	1361	& ' but no biogeochemical model appears to have been defined' )
	1362	#endif
	1363	llog10 = .TRUE.
	1364
	1365	CASE('slphydia')
	1366	#if defined key_hadocc
	1367	CALL ctl_stop( ' Trying to run slphydia observation operator', &
	1368	& ' but HadOCC does not explicitly simulate diatoms' )
[10168]	1369	#elif defined key_medusa
[9306]	1370	! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
	1371	zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd
	1372	#elif defined key_fabm
	1373	! Diatom surface phytoplankton carbon from ERSEM
[10729]	1374	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)
[9306]	1375	#else
	1376	CALL ctl_stop( ' Trying to run slphydia observation operator', &
	1377	& ' but no biogeochemical model appears to have been defined' )
	1378	#endif
	1379	llog10 = .TRUE.
	1380
	1381	CASE('slphynon')
	1382	#if defined key_hadocc
	1383	CALL ctl_stop( ' Trying to run slphynon observation operator', &
	1384	& ' but HadOCC does not explicitly simulate non-diatoms' )
[10168]	1385	#elif defined key_medusa
[9306]	1386	! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
	1387	zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn
	1388	#elif defined key_fabm
	1389	! Add all non-diatom surface phytoplankton carbon groups from ERSEM
[10729]	1390	zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
	1391	& trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
[9306]	1392	#else
	1393	CALL ctl_stop( ' Trying to run slphynon observation operator', &
	1394	& ' but no biogeochemical model appears to have been defined' )
	1395	#endif
	1396	llog10 = .TRUE.
	1397
	1398	CASE('sspm')
[7992]	1399	#if defined key_spm
	1400	zsurfvar(:,:) = 0.0
	1401	DO jn = 1, jp_spm
	1402	zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes
	1403	END DO
	1404	#else
[9306]	1405	CALL ctl_stop( ' Trying to run sspm observation operator', &
[7992]	1406	& ' but no spm model appears to have been defined' )
	1407	#endif
[9306]	1408
	1409	CASE('sfco2')
[7992]	1410	#if defined key_hadocc
	1411	zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC
	1412	IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
	1413	& ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
	1414	zsurfvar(:,:) = obfillflt
	1415	zsurfmask(:,:) = 0
	1416	CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
	1417	& ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
	1418	ENDIF
[10168]	1419	#elif defined key_medusa && defined key_roam
[8653]	1420	zsurfvar(:,:) = f2_fco2w(:,:)
[7992]	1421	#elif defined key_fabm
	1422	! First, get pCO2 from FABM
[10729]	1423	pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
[7992]	1424	zsurfvar(:,:) = pco2_3d(:,:,1)
	1425	! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
	1426	! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
	1427	! and data reduction routines, Deep-Sea Research II, 56: 512-522.
	1428	! and
	1429	! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
	1430	! Marine Chemistry, 2: 203-215.
	1431	! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
	1432	! not explicitly included - atmospheric pressure is not necessarily available so this is
	1433	! the best assumption.
	1434	! Further, the (1-xCO2)^2 term has been neglected. This is common practice
	1435	! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
	1436	! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
	1437	! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
	1438	zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + &
	1439	& 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - &
	1440	& 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
	1441	& 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0) + &
	1442	& 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / &
	1443	& (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
	1444	#else
[9306]	1445	CALL ctl_stop( ' Trying to run sfco2 observation operator', &
[7992]	1446	& ' but no biogeochemical model appears to have been defined' )
	1447	#endif
[9306]	1448
	1449	CASE('spco2')
[7992]	1450	#if defined key_hadocc
	1451	zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC
	1452	IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
	1453	& ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
	1454	zsurfvar(:,:) = obfillflt
	1455	zsurfmask(:,:) = 0
	1456	CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
	1457	& ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
	1458	ENDIF
[10168]	1459	#elif defined key_medusa && defined key_roam
[8653]	1460	zsurfvar(:,:) = f2_pco2w(:,:)
[7992]	1461	#elif defined key_fabm
[10729]	1462	pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
[7992]	1463	zsurfvar(:,:) = pco2_3d(:,:,1)
	1464	#else
[9306]	1465	CALL ctl_stop( ' Trying to run spco2 observation operator', &
[7992]	1466	& ' but no biogeochemical model appears to have been defined' )
	1467	#endif
	1468
	1469	CASE DEFAULT
	1470
	1471	CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
	1472
	1473	END SELECT
[9306]	1474
	1475	IF ( llog10 ) THEN
	1476	! Take the log10 where we can, otherwise exclude
	1477	tiny = 1.0e-20
	1478	WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
	1479	zsurfvar(:,:) = LOG10(zsurfvar(:,:))
	1480	ELSEWHERE
	1481	zsurfvar(:,:) = obfillflt
	1482	zsurfmask(:,:) = 0
	1483	END WHERE
	1484	ENDIF
[7992]	1485
	1486	CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, &
	1487	& nit000, idaystp, zsurfvar, zsurfmask, &
	1488	& n2dintsurf(jtype), llnightav(jtype), &
	1489	& ravglamscl(jtype), ravgphiscl(jtype), &
	1490	& lfpindegs(jtype) )
	1491
[2128]	1492	END DO
[7992]	1493
[9306]	1494	CALL wrk_dealloc( jpi, jpj, zsurfvar )
	1495	CALL wrk_dealloc( jpi, jpj, zsurfmask )
[10729]	1496	#if defined key_fabm
	1497	CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d )
	1498	#endif
[9306]	1499
[2128]	1500	ENDIF
	1501
	1502	END SUBROUTINE dia_obs
[7992]	1503
	1504	SUBROUTINE dia_obs_wri
[2128]	1505	!!----------------------------------------------------------------------
	1506	!! * ROUTINE dia_obs_wri *
	1507	!!
	1508	!! ** Purpose : Call observation diagnostic output routines
	1509	!!
	1510	!! ** Method : Call observation diagnostic output routines
	1511	!!
[7992]	1512	!! ** Action :
[2128]	1513	!!
	1514	!! History :
	1515	!! ! 06-03 (K. Mogensen) Original code
	1516	!! ! 06-05 (K. Mogensen) Reformatted
	1517	!! ! 06-10 (A. Weaver) Cleaning
	1518	!! ! 07-03 (K. Mogensen) General handling of profiles
	1519	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
[7992]	1520	!! ! 15-08 (M. Martin) Combined writing for prof and surf types
[2128]	1521	!!----------------------------------------------------------------------
[7992]	1522	!! * Modules used
	1523	USE obs_rot_vel ! Rotation of velocities
	1524
[2128]	1525	IMPLICIT NONE
	1526
	1527	!! * Local declarations
[7992]	1528	INTEGER :: jtype ! Data set loop variable
	1529	INTEGER :: jo, jvar, jk
	1530	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	1531	& zu, &
	1532	& zv
[2128]	1533
	1534	!-----------------------------------------------------------------------
	1535	! Depending on switches call various observation output routines
	1536	!-----------------------------------------------------------------------
	1537
[7992]	1538	IF ( nproftypes > 0 ) THEN
[2128]	1539
[7992]	1540	DO jtype = 1, nproftypes
[2128]	1541
[7992]	1542	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
[2128]	1543
[7992]	1544	! For velocity data, rotate the model velocities to N/S, E/W
	1545	! using the compressed data structure.
	1546	ALLOCATE( &
	1547	& zu(profdataqc(jtype)%nvprot(1)), &
	1548	& zv(profdataqc(jtype)%nvprot(2)) &
	1549	& )
[2128]	1550
[9306]	1551	CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
[2128]	1552
[7992]	1553	DO jo = 1, profdataqc(jtype)%nprof
	1554	DO jvar = 1, 2
	1555	DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
[2128]	1556
[7992]	1557	IF ( jvar == 1 ) THEN
	1558	profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
	1559	ELSE
	1560	profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
	1561	ENDIF
[2128]	1562
[7992]	1563	END DO
	1564	END DO
	1565	END DO
[2128]	1566
[7992]	1567	DEALLOCATE( zu )
	1568	DEALLOCATE( zv )
[2128]	1569
[7992]	1570	END IF
[2128]	1571
[7992]	1572	CALL obs_prof_decompress( profdataqc(jtype), &
	1573	& profdata(jtype), .TRUE., numout )
[2128]	1574
[7992]	1575	CALL obs_wri_prof( profdata(jtype) )
[2128]	1576
	1577	END DO
	1578
	1579	ENDIF
	1580
[7992]	1581	IF ( nsurftypes > 0 ) THEN
[2128]	1582
[7992]	1583	DO jtype = 1, nsurftypes
[2128]	1584
[7992]	1585	CALL obs_surf_decompress( surfdataqc(jtype), &
	1586	& surfdata(jtype), .TRUE., numout )
[2128]	1587
[7992]	1588	CALL obs_wri_surf( surfdata(jtype) )
[2128]	1589
	1590	END DO
	1591
	1592	ENDIF
	1593
	1594	END SUBROUTINE dia_obs_wri
	1595
[4245]	1596	SUBROUTINE dia_obs_dealloc
	1597	IMPLICIT NONE
	1598	!!----------------------------------------------------------------------
	1599	!! * ROUTINE dia_obs_dealloc *
	1600	!!
	1601	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
	1602	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
	1603	!!
	1604	!! ** Method : Clean up various arrays left behind by the obs_oper.
	1605	!!
	1606	!! ** Action :
	1607	!!
	1608	!!----------------------------------------------------------------------
[7992]	1609	! obs_grid deallocation
[4245]	1610	CALL obs_grid_deallocate
	1611
[7992]	1612	! diaobs deallocation
	1613	IF ( nproftypes > 0 ) &
	1614	& DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
	1615
	1616	IF ( nsurftypes > 0 ) &
	1617	& DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
	1618	& n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
	1619
[4245]	1620	END SUBROUTINE dia_obs_dealloc
	1621
[2128]	1622	SUBROUTINE ini_date( ddobsini )
	1623	!!----------------------------------------------------------------------
	1624	!! * ROUTINE ini_date *
	1625	!!
[7992]	1626	!! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
[2128]	1627	!!
[7992]	1628	!! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format
[2128]	1629	!!
[7992]	1630	!! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format
	1631	!!
[2128]	1632	!! History :
	1633	!! ! 06-03 (K. Mogensen) Original code
	1634	!! ! 06-05 (K. Mogensen) Reformatted
	1635	!! ! 06-10 (A. Weaver) Cleaning
	1636	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
	1637	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
	1638	!!----------------------------------------------------------------------
	1639	USE phycst, ONLY : & ! Physical constants
	1640	& rday
	1641	USE dom_oce, ONLY : & ! Ocean space and time domain variables
	1642	& rdt
	1643
	1644	IMPLICIT NONE
	1645
	1646	!! * Arguments
[7992]	1647	REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
[2128]	1648
	1649	!! * Local declarations
	1650	INTEGER :: iyea ! date - (year, month, day, hour, minute)
	1651	INTEGER :: imon
	1652	INTEGER :: iday
	1653	INTEGER :: ihou
	1654	INTEGER :: imin
[7992]	1655	INTEGER :: imday ! Number of days in month.
	1656	INTEGER, DIMENSION(12) :: &
	1657	& imonth_len ! Length in days of the months of the current year
	1658	REAL(wp) :: zdayfrc ! Fraction of day
[2128]	1659
[7992]	1660	!----------------------------------------------------------------------
	1661	! Initial date initialization (year, month, day, hour, minute)
	1662	! (This assumes that the initial date is for 00z))
	1663	!----------------------------------------------------------------------
[2128]	1664	iyea = ndate0 / 10000
	1665	imon = ( ndate0 - iyea * 10000 ) / 100
	1666	iday = ndate0 - iyea * 10000 - imon * 100
	1667	ihou = 0
	1668	imin = 0
	1669
[7992]	1670	!----------------------------------------------------------------------
	1671	! Compute number of days + number of hours + min since initial time
	1672	!----------------------------------------------------------------------
[2128]	1673	iday = iday + ( nit000 -1 ) * rdt / rday
	1674	zdayfrc = ( nit000 -1 ) * rdt / rday
	1675	zdayfrc = zdayfrc - aint(zdayfrc)
	1676	ihou = int( zdayfrc * 24 )
	1677	imin = int( (zdayfrc * 24 - ihou) * 60 )
	1678
[7992]	1679	!-----------------------------------------------------------------------
	1680	! Convert number of days (iday) into a real date
	1681	!----------------------------------------------------------------------
[2128]	1682
	1683	CALL calc_month_len( iyea, imonth_len )
[7992]	1684
[2128]	1685	DO WHILE ( iday > imonth_len(imon) )
	1686	iday = iday - imonth_len(imon)
	1687	imon = imon + 1
	1688	IF ( imon > 12 ) THEN
	1689	imon = 1
	1690	iyea = iyea + 1
	1691	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
	1692	ENDIF
	1693	END DO
	1694
[7992]	1695	!----------------------------------------------------------------------
	1696	! Convert it into YYYYMMDD.HHMMSS format.
	1697	!----------------------------------------------------------------------
[2128]	1698	ddobsini = iyea * 10000_dp + imon * 100_dp + &
	1699	& iday + ihou * 0.01_dp + imin * 0.0001_dp
	1700
	1701
	1702	END SUBROUTINE ini_date
	1703
	1704	SUBROUTINE fin_date( ddobsfin )
	1705	!!----------------------------------------------------------------------
	1706	!! * ROUTINE fin_date *
	1707	!!
[7992]	1708	!! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
[2128]	1709	!!
[7992]	1710	!! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format
[2128]	1711	!!
[7992]	1712	!! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format
	1713	!!
[2128]	1714	!! History :
	1715	!! ! 06-03 (K. Mogensen) Original code
	1716	!! ! 06-05 (K. Mogensen) Reformatted
	1717	!! ! 06-10 (A. Weaver) Cleaning
	1718	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
	1719	!!----------------------------------------------------------------------
	1720	USE phycst, ONLY : & ! Physical constants
	1721	& rday
	1722	USE dom_oce, ONLY : & ! Ocean space and time domain variables
	1723	& rdt
	1724
	1725	IMPLICIT NONE
	1726
	1727	!! * Arguments
[7992]	1728	REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
[2128]	1729
	1730	!! * Local declarations
	1731	INTEGER :: iyea ! date - (year, month, day, hour, minute)
	1732	INTEGER :: imon
	1733	INTEGER :: iday
	1734	INTEGER :: ihou
	1735	INTEGER :: imin
[7992]	1736	INTEGER :: imday ! Number of days in month.
	1737	INTEGER, DIMENSION(12) :: &
	1738	& imonth_len ! Length in days of the months of the current year
	1739	REAL(wp) :: zdayfrc ! Fraction of day
	1740
[2128]	1741	!-----------------------------------------------------------------------
	1742	! Initial date initialization (year, month, day, hour, minute)
	1743	! (This assumes that the initial date is for 00z)
	1744	!-----------------------------------------------------------------------
	1745	iyea = ndate0 / 10000
	1746	imon = ( ndate0 - iyea * 10000 ) / 100
	1747	iday = ndate0 - iyea * 10000 - imon * 100
	1748	ihou = 0
	1749	imin = 0
[7992]	1750
[2128]	1751	!-----------------------------------------------------------------------
	1752	! Compute number of days + number of hours + min since initial time
	1753	!-----------------------------------------------------------------------
	1754	iday = iday + nitend * rdt / rday
	1755	zdayfrc = nitend * rdt / rday
	1756	zdayfrc = zdayfrc - AINT( zdayfrc )
	1757	ihou = INT( zdayfrc * 24 )
	1758	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
	1759
	1760	!-----------------------------------------------------------------------
	1761	! Convert number of days (iday) into a real date
	1762	!----------------------------------------------------------------------
	1763
	1764	CALL calc_month_len( iyea, imonth_len )
[7992]	1765
[2128]	1766	DO WHILE ( iday > imonth_len(imon) )
	1767	iday = iday - imonth_len(imon)
	1768	imon = imon + 1
	1769	IF ( imon > 12 ) THEN
	1770	imon = 1
	1771	iyea = iyea + 1
	1772	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
	1773	ENDIF
	1774	END DO
	1775
	1776	!-----------------------------------------------------------------------
	1777	! Convert it into YYYYMMDD.HHMMSS format
	1778	!-----------------------------------------------------------------------
	1779	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
	1780	& + ihou * 0.01_dp + imin * 0.0001_dp
	1781
	1782	END SUBROUTINE fin_date
[7992]	1783
[9306]	1784	SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
[7992]	1785
[9306]	1786	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
	1787	INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
	1788	INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
	1789	& ifiles ! Out number of files for each type
	1790	CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
	1791	& cobstypes ! List of obs types
	1792	CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
	1793	& cfiles ! List of files for all types
[7992]	1794
[9306]	1795	!Local variables
	1796	INTEGER :: jfile
	1797	INTEGER :: jtype
[7992]	1798
[9306]	1799	DO jtype = 1, ntypes
[7992]	1800
[9306]	1801	ifiles(jtype) = 0
	1802	DO jfile = 1, jpmaxnfiles
	1803	IF ( trim(cfiles(jtype,jfile)) /= '' ) &
	1804	ifiles(jtype) = ifiles(jtype) + 1
	1805	END DO
[7992]	1806
[9306]	1807	IF ( ifiles(jtype) == 0 ) THEN
	1808	CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))// &
	1809	& ' set to true but no files available to read' )
	1810	ENDIF
[7992]	1811
[9306]	1812	IF(lwp) THEN
	1813	WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:'
	1814	DO jfile = 1, ifiles(jtype)
	1815	WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile))
	1816	END DO
	1817	ENDIF
	1818
[7992]	1819	END DO
	1820
	1821	END SUBROUTINE obs_settypefiles
	1822
	1823	SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, &
	1824	& n2dint_default, n2dint_type, &
	1825	& ravglamscl_type, ravgphiscl_type, &
	1826	& lfp_indegs_type, lavnight_type, &
	1827	& n2dint, ravglamscl, ravgphiscl, &
	1828	& lfpindegs, lavnight )
	1829
[9306]	1830	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
	1831	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
	1832	INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type
	1833	INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type
	1834	REAL(wp), INTENT(IN) :: &
	1835	& ravglamscl_type, & !E/W diameter of obs footprint for this type
	1836	& ravgphiscl_type !N/S diameter of obs footprint for this type
	1837	LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres
	1838	LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average
	1839	CHARACTER(len=8), INTENT(IN) :: ctypein
[7992]	1840
[9306]	1841	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
	1842	& n2dint
	1843	REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
	1844	& ravglamscl, ravgphiscl
	1845	LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
	1846	& lfpindegs, lavnight
[7992]	1847
[9306]	1848	lavnight(jtype) = lavnight_type
[7992]	1849
[9306]	1850	IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN
	1851	n2dint(jtype) = n2dint_type
	1852	ELSE IF ( n2dint_type == -1 ) THEN
	1853	n2dint(jtype) = n2dint_default
[7992]	1854	ELSE
[9306]	1855	CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
	1856	& ' is not available')
[7992]	1857	ENDIF
	1858
[9306]	1859	! For averaging observation footprints set options for size of footprint
	1860	IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
	1861	IF ( ravglamscl_type > 0._wp ) THEN
	1862	ravglamscl(jtype) = ravglamscl_type
	1863	ELSE
	1864	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
	1865	'scale (ravglamscl) for observation type '//TRIM(ctypein) )
	1866	ENDIF
[7992]	1867
[9306]	1868	IF ( ravgphiscl_type > 0._wp ) THEN
	1869	ravgphiscl(jtype) = ravgphiscl_type
	1870	ELSE
	1871	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
	1872	'scale (ravgphiscl) for observation type '//TRIM(ctypein) )
	1873	ENDIF
[7992]	1874
[9306]	1875	lfpindegs(jtype) = lfp_indegs_type
[7992]	1876
[9306]	1877	ENDIF
	1878
	1879	! Write out info
	1880	IF(lwp) THEN
	1881	IF ( n2dint(jtype) <= 4 ) THEN
	1882	WRITE(numout,*) ' '//TRIM(ctypein)// &
	1883	& ' model counterparts will be interpolated horizontally'
	1884	ELSE IF ( n2dint(jtype) <= 6 ) THEN
	1885	WRITE(numout,*) ' '//TRIM(ctypein)// &
	1886	& ' model counterparts will be averaged horizontally'
	1887	WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype)
	1888	WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype)
	1889	IF ( lfpindegs(jtype) ) THEN
	1890	WRITE(numout,*) ' '//' (in degrees)'
	1891	ELSE
	1892	WRITE(numout,*) ' '//' (in metres)'
	1893	ENDIF
[7992]	1894	ENDIF
	1895	ENDIF
	1896
	1897	END SUBROUTINE obs_setinterpopts
	1898
[2128]	1899	END MODULE diaobs

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

Download in other formats: