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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_medusa/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8646

Last change on this file since 8646 was 8646, checked in by dford, 7 years ago
Add comments in response to review.
File size: 58.5 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 :: &
	47	& lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6.
	48	LOGICAL :: ln_diaobs !: Logical switch for the obs operator
	49	LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs
	50	LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres
	51	LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres
	52	LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres
	53	LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
[2128]	54
[7992]	55	REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
	56	REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
	57	REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
	58	REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
	59	REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
	60	REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
	61	REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
	62	REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
[2128]	63
[7992]	64	INTEGER :: nn_1dint !: Vertical interpolation method
	65	INTEGER :: nn_2dint !: Default horizontal interpolation method
	66	INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method
	67	INTEGER :: nn_2dint_sst !: SST horizontal interpolation method
	68	INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method
	69	INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method
	70
[2128]	71	INTEGER, DIMENSION(imaxavtypes) :: &
[7992]	72	& nn_profdavtypes !: Profile data types representing a daily average
	73	INTEGER :: nproftypes !: Number of profile obs types
	74	INTEGER :: nsurftypes !: Number of surface obs types
	75	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	76	& nvarsprof, & !: Number of profile variables
	77	& nvarssurf !: Number of surface variables
	78	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	79	& nextrprof, & !: Number of profile extra variables
	80	& nextrsurf !: Number of surface extra variables
	81	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	82	& n2dintsurf !: Interpolation option for surface variables
	83	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	84	& ravglamscl, & !: E/W diameter of averaging footprint for surface variables
	85	& ravgphiscl !: N/S diameter of averaging footprint for surface variables
[2128]	86	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
[7992]	87	& lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres
	88	& llnightav !: Logical for calculating night-time averages
[2128]	89
[7992]	90	TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
	91	& surfdata, & !: Initial surface data
	92	& surfdataqc !: Surface data after quality control
	93	TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
	94	& profdata, & !: Initial profile data
	95	& profdataqc !: Profile data after quality control
	96
	97	CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
	98	& cobstypesprof, & !: Profile obs types
	99	& cobstypessurf !: Surface obs types
	100
[2287]	101	!!----------------------------------------------------------------------
	102	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	103	!! $Id$
	104	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	105	!!----------------------------------------------------------------------
	106
[7992]	107	!! * Substitutions
	108	# include "domzgr_substitute.h90"
[2128]	109	CONTAINS
	110
	111	SUBROUTINE dia_obs_init
	112	!!----------------------------------------------------------------------
	113	!! * ROUTINE dia_obs_init *
	114	!!
	115	!! ** Purpose : Initialize and read observations
	116	!!
	117	!! ** Method : Read the namelist and call reading routines
	118	!!
	119	!! ** Action : Read the namelist and call reading routines
	120	!!
	121	!! History :
	122	!! ! 06-03 (K. Mogensen) Original code
	123	!! ! 06-05 (A. Weaver) Reformatted
	124	!! ! 06-10 (A. Weaver) Cleaning and add controls
	125	!! ! 07-03 (K. Mogensen) General handling of profiles
[7992]	126	!! ! 14-08 (J.While) Incorporated SST bias correction
	127	!! ! 15-02 (M. Martin) Simplification of namelist and code
[2128]	128	!!----------------------------------------------------------------------
	129
	130	IMPLICIT NONE
	131
	132	!! * Local declarations
[7992]	133	INTEGER, PARAMETER :: &
	134	& jpmaxnfiles = 1000 ! Maximum number of files for each obs type
	135	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	136	& ifilesprof, & ! Number of profile files
	137	& ifilessurf ! Number of surface files
	138	INTEGER :: ios ! Local integer output status for namelist read
	139	INTEGER :: jtype ! Counter for obs types
	140	INTEGER :: jvar ! Counter for variables
	141	INTEGER :: jfile ! Counter for files
	142	INTEGER :: jnumsstbias ! Number of SST bias files to read and apply
	143
	144	CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
	145	& cn_profbfiles, & ! T/S profile input filenames
	146	& cn_sstfbfiles, & ! Sea surface temperature input filenames
	147	& cn_slafbfiles, & ! Sea level anomaly input filenames
	148	& cn_sicfbfiles, & ! Seaice concentration input filenames
	149	& cn_velfbfiles, & ! Velocity profile input filenames
	150	& cn_sssfbfiles, & ! Sea surface salinity input filenames
	151	& cn_logchlfbfiles, & ! Log(Chl) input filenames
	152	& cn_spmfbfiles, & ! Sediment input filenames
	153	& cn_fco2fbfiles, & ! fco2 input filenames
	154	& cn_pco2fbfiles, & ! pco2 input filenames
	155	& cn_sstbiasfiles ! SST bias input filenames
	156
	157	CHARACTER(LEN=128) :: &
	158	& cn_altbiasfile ! Altimeter bias input filename
	159
	160
	161	LOGICAL :: ln_t3d ! Logical switch for temperature profiles
	162	LOGICAL :: ln_s3d ! Logical switch for salinity profiles
	163	LOGICAL :: ln_sla ! Logical switch for sea level anomalies
	164	LOGICAL :: ln_sst ! Logical switch for sea surface temperature
	165	LOGICAL :: ln_sic ! Logical switch for sea ice concentration
	166	LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs
	167	LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs
	168	LOGICAL :: ln_logchl ! Logical switch for log(Chl) obs
	169	LOGICAL :: ln_spm ! Logical switch for sediment obs
	170	LOGICAL :: ln_fco2 ! Logical switch for fco2 obs
	171	LOGICAL :: ln_pco2 ! Logical switch for pco2 obs
	172	LOGICAL :: ln_nea ! Logical switch to remove obs near land
	173	LOGICAL :: ln_altbias ! Logical switch for altimeter bias
	174	LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST
	175	LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files
	176	LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs
	177	LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
	178
	179	REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS
	180	REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS
	181
	182	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
	183	& clproffiles, & ! Profile filenames
	184	& clsurffiles ! Surface filenames
	185
	186	LOGICAL :: llvar1 ! Logical for profile variable 1
	187	LOGICAL :: llvar2 ! Logical for profile variable 1
	188
	189	REAL(wp), POINTER, DIMENSION(:,:) :: &
	190	& zglam1, & ! Model longitudes for profile variable 1
	191	& zglam2 ! Model longitudes for profile variable 2
	192	REAL(wp), POINTER, DIMENSION(:,:) :: &
	193	& zgphi1, & ! Model latitudes for profile variable 1
	194	& zgphi2 ! Model latitudes for profile variable 2
	195	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
	196	& zmask1, & ! Model land/sea mask associated with variable 1
	197	& zmask2 ! Model land/sea mask associated with variable 2
	198
	199
	200	NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, &
	201	& ln_sst, ln_sic, ln_sss, ln_vel3d, &
	202	& ln_logchl, ln_spm, ln_fco2, ln_pco2, &
	203	& ln_altbias, ln_sstbias, ln_nea, &
	204	& ln_grid_global, ln_grid_search_lookup, &
	205	& ln_ignmis, ln_s_at_t, ln_bound_reject, &
[4245]	206	& ln_sstnight, &
[7992]	207	& ln_sla_fp_indegs, ln_sst_fp_indegs, &
	208	& ln_sss_fp_indegs, ln_sic_fp_indegs, &
	209	& cn_profbfiles, cn_slafbfiles, &
	210	& cn_sstfbfiles, cn_sicfbfiles, &
	211	& cn_velfbfiles, cn_sssfbfiles, &
	212	& cn_logchlfbfiles, cn_spmfbfiles, &
	213	& cn_fco2fbfiles, cn_pco2fbfiles, &
	214	& cn_sstbiasfiles, cn_altbiasfile, &
	215	& cn_gridsearchfile, rn_gridsearchres, &
	216	& rn_dobsini, rn_dobsend, &
	217	& rn_sla_avglamscl, rn_sla_avgphiscl, &
	218	& rn_sst_avglamscl, rn_sst_avgphiscl, &
	219	& rn_sss_avglamscl, rn_sss_avgphiscl, &
	220	& rn_sic_avglamscl, rn_sic_avgphiscl, &
	221	& nn_1dint, nn_2dint, &
	222	& nn_2dint_sla, nn_2dint_sst, &
	223	& nn_2dint_sss, nn_2dint_sic, &
	224	& nn_msshc, rn_mdtcorr, rn_mdtcutoff, &
	225	& nn_profdavtypes
[2128]	226
[7992]	227	CALL wrk_alloc( jpi, jpj, zglam1 )
	228	CALL wrk_alloc( jpi, jpj, zglam2 )
	229	CALL wrk_alloc( jpi, jpj, zgphi1 )
	230	CALL wrk_alloc( jpi, jpj, zgphi2 )
	231	CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
	232	CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
[2128]	233
	234	!-----------------------------------------------------------------------
	235	! Read namelist parameters
	236	!-----------------------------------------------------------------------
	237
[7992]	238	! Some namelist arrays need initialising
	239	cn_profbfiles(:) = ''
	240	cn_slafbfiles(:) = ''
	241	cn_sstfbfiles(:) = ''
	242	cn_sicfbfiles(:) = ''
	243	cn_velfbfiles(:) = ''
	244	cn_sssfbfiles(:) = ''
	245	cn_logchlfbfiles(:) = ''
	246	cn_spmfbfiles(:) = ''
	247	cn_fco2fbfiles(:) = ''
	248	cn_pco2fbfiles(:) = ''
	249	cn_sstbiasfiles(:) = ''
	250	nn_profdavtypes(:) = -1
	251
	252	CALL ini_date( rn_dobsini )
	253	CALL fin_date( rn_dobsend )
	254
	255	! Read namelist namobs : control observation diagnostics
	256	REWIND( numnam_ref ) ! Namelist namobs in reference namelist
[4147]	257	READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
	258	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
[2128]	259
[7992]	260	REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist
[4147]	261	READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
	262	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
[4624]	263	IF(lwm) WRITE ( numond, namobs )
[4147]	264
[7992]	265	lk_diaobs = .FALSE.
	266	#if defined key_diaobs
	267	IF ( ln_diaobs ) lk_diaobs = .TRUE.
	268	#endif
	269
	270	IF ( .NOT. lk_diaobs ) THEN
	271	IF(lwp) WRITE(numout,cform_war)
	272	IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs'
	273	RETURN
[2128]	274	ENDIF
[7992]	275
[2128]	276	IF(lwp) THEN
	277	WRITE(numout,*)
	278	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
	279	WRITE(numout,*) '~~~~~~~~~~~~'
	280	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
[7992]	281	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
	282	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
	283	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
	284	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
	285	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic
	286	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
	287	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
	288	WRITE(numout,*) ' Logical switch for log(Chl) observations ln_logchl = ', ln_logchl
	289	WRITE(numout,*) ' Logical switch for SPM observations ln_spm = ', ln_spm
	290	WRITE(numout,*) ' Logical switch for FCO2 observations ln_fco2 = ', ln_fco2
	291	WRITE(numout,*) ' Logical switch for PCO2 observations ln_pco2 = ', ln_pco2
	292	WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global
	293	WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
[2128]	294	IF (ln_grid_search_lookup) &
[7992]	295	WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile
	296	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini
	297	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend
	298	WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint
	299	WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint
	300	WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea
	301	WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject
	302	WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc
	303	WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr
	304	WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff
	305	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
	306	WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias
	307	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
	308	WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes
	309	WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight
[2128]	310	ENDIF
[7992]	311	!-----------------------------------------------------------------------
	312	! Set up list of observation types to be used
	313	! and the files associated with each type
	314	!-----------------------------------------------------------------------
[2128]	315
[7992]	316	nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
	317	nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, &
	318	& ln_logchl, ln_spm, ln_fco2, ln_pco2 /) )
	319
	320	IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
[2128]	321	IF(lwp) WRITE(numout,cform_war)
[7992]	322	IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
	323	& ' are set to .FALSE. so turning off calls to dia_obs'
[2128]	324	nwarn = nwarn + 1
[7992]	325	lk_diaobs = .FALSE.
	326	RETURN
[2128]	327	ENDIF
	328
[7992]	329	IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes
	330	IF ( nproftypes > 0 ) THEN
[2128]	331
[7992]	332	ALLOCATE( cobstypesprof(nproftypes) )
	333	ALLOCATE( ifilesprof(nproftypes) )
	334	ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
[2128]	335
[7992]	336	jtype = 0
	337	IF (ln_t3d .OR. ln_s3d) THEN
	338	jtype = jtype + 1
	339	CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', &
	340	& cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
[2128]	341	ENDIF
[7992]	342	IF (ln_vel3d) THEN
	343	jtype = jtype + 1
	344	CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', &
	345	& cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
[2128]	346	ENDIF
	347
[7992]	348	ENDIF
[2128]	349
[7992]	350	IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes
	351	IF ( nsurftypes > 0 ) THEN
[2128]	352
[7992]	353	ALLOCATE( cobstypessurf(nsurftypes) )
	354	ALLOCATE( ifilessurf(nsurftypes) )
	355	ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
	356	ALLOCATE(n2dintsurf(nsurftypes))
	357	ALLOCATE(ravglamscl(nsurftypes))
	358	ALLOCATE(ravgphiscl(nsurftypes))
	359	ALLOCATE(lfpindegs(nsurftypes))
	360	ALLOCATE(llnightav(nsurftypes))
[2128]	361
[7992]	362	jtype = 0
	363	IF (ln_sla) THEN
	364	jtype = jtype + 1
	365	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', &
	366	& cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
	367	CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', &
	368	& nn_2dint, nn_2dint_sla, &
	369	& rn_sla_avglamscl, rn_sla_avgphiscl, &
	370	& ln_sla_fp_indegs, .FALSE., &
	371	& n2dintsurf, ravglamscl, ravgphiscl, &
	372	& lfpindegs, llnightav )
[2128]	373	ENDIF
[7992]	374	IF (ln_sst) THEN
	375	jtype = jtype + 1
	376	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', &
	377	& cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
	378	CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', &
	379	& nn_2dint, nn_2dint_sst, &
	380	& rn_sst_avglamscl, rn_sst_avgphiscl, &
	381	& ln_sst_fp_indegs, ln_sstnight, &
	382	& n2dintsurf, ravglamscl, ravgphiscl, &
	383	& lfpindegs, llnightav )
[2128]	384	ENDIF
[7992]	385	#if defined key_lim2 \|\| defined key_lim3 \|\| defined key_cice
	386	IF (ln_sic) THEN
	387	jtype = jtype + 1
	388	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', &
	389	& cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
	390	CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', &
	391	& nn_2dint, nn_2dint_sic, &
	392	& rn_sic_avglamscl, rn_sic_avgphiscl, &
	393	& ln_sic_fp_indegs, .FALSE., &
	394	& n2dintsurf, ravglamscl, ravgphiscl, &
	395	& lfpindegs, llnightav )
[2128]	396	ENDIF
[7992]	397	#endif
	398	IF (ln_sss) THEN
	399	jtype = jtype + 1
	400	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', &
	401	& cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
	402	CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', &
	403	& nn_2dint, nn_2dint_sss, &
	404	& rn_sss_avglamscl, rn_sss_avgphiscl, &
	405	& ln_sss_fp_indegs, .FALSE., &
	406	& n2dintsurf, ravglamscl, ravgphiscl, &
	407	& lfpindegs, llnightav )
[2128]	408	ENDIF
	409
[7992]	410	IF (ln_logchl) THEN
	411	jtype = jtype + 1
	412	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', &
	413	& cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles )
	414	CALL obs_setinterpopts( nsurftypes, jtype, 'logchl', &
	415	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
	416	& n2dintsurf, ravglamscl, ravgphiscl, &
	417	& lfpindegs, llnightav )
	418	ENDIF
[2128]	419
[7992]	420	IF (ln_spm) THEN
	421	jtype = jtype + 1
	422	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm ', &
	423	& cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles )
	424	CALL obs_setinterpopts( nsurftypes, jtype, 'spm ', &
	425	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
	426	& n2dintsurf, ravglamscl, ravgphiscl, &
	427	& lfpindegs, llnightav )
	428	ENDIF
[2128]	429
[7992]	430	IF (ln_fco2) THEN
	431	jtype = jtype + 1
	432	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2 ', &
	433	& cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
	434	CALL obs_setinterpopts( nsurftypes, jtype, 'fco2 ', &
	435	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
	436	& n2dintsurf, ravglamscl, ravgphiscl, &
	437	& lfpindegs, llnightav )
[2128]	438	ENDIF
	439
[7992]	440	IF (ln_pco2) THEN
	441	jtype = jtype + 1
	442	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2 ', &
	443	& cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
	444	CALL obs_setinterpopts( nsurftypes, jtype, 'pco2 ', &
	445	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
	446	& n2dintsurf, ravglamscl, ravgphiscl, &
	447	& lfpindegs, llnightav )
[2128]	448	ENDIF
	449
	450	ENDIF
	451
[7992]	452	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[2128]	453
	454
[7992]	455	!-----------------------------------------------------------------------
	456	! Obs operator parameter checking and initialisations
	457	!-----------------------------------------------------------------------
[2128]	458
[7992]	459	IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
	460	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
	461	RETURN
	462	ENDIF
[2128]	463
[7992]	464	IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
	465	CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
	466	& ' is not available')
	467	ENDIF
[2128]	468
[7992]	469	IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN
	470	CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
	471	& ' is not available')
	472	ENDIF
[2128]	473
[7992]	474	CALL obs_typ_init
[2128]	475
[7992]	476	CALL mppmap_init
[2128]	477
[7992]	478	CALL obs_grid_setup( )
[2128]	479
[7992]	480	!-----------------------------------------------------------------------
	481	! Depending on switches read the various observation types
	482	!-----------------------------------------------------------------------
[3651]	483
[7992]	484	IF ( nproftypes > 0 ) THEN
[2128]	485
[7992]	486	ALLOCATE(profdata(nproftypes))
	487	ALLOCATE(profdataqc(nproftypes))
	488	ALLOCATE(nvarsprof(nproftypes))
	489	ALLOCATE(nextrprof(nproftypes))
[3651]	490
[7992]	491	DO jtype = 1, nproftypes
[2128]	492
[7992]	493	nvarsprof(jtype) = 2
	494	IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
	495	nextrprof(jtype) = 1
	496	llvar1 = ln_t3d
	497	llvar2 = ln_s3d
	498	zglam1 = glamt
	499	zgphi1 = gphit
	500	zmask1 = tmask
	501	zglam2 = glamt
	502	zgphi2 = gphit
	503	zmask2 = tmask
	504	ENDIF
	505	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
	506	nextrprof(jtype) = 2
	507	llvar1 = ln_vel3d
	508	llvar2 = ln_vel3d
	509	zglam1 = glamu
	510	zgphi1 = gphiu
	511	zmask1 = umask
	512	zglam2 = glamv
	513	zgphi2 = gphiv
	514	zmask2 = vmask
	515	ENDIF
[2128]	516
[7992]	517	!Read in profile or profile obs types
	518	CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), &
	519	& clproffiles(jtype,1:ifilesprof(jtype)), &
	520	& nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
	521	& rn_dobsini, rn_dobsend, llvar1, llvar2, &
	522	& ln_ignmis, ln_s_at_t, .FALSE., &
	523	& kdailyavtypes = nn_profdavtypes )
[2128]	524
[7992]	525	DO jvar = 1, nvarsprof(jtype)
	526	CALL obs_prof_staend( profdata(jtype), jvar )
	527	END DO
[3651]	528
[7992]	529	CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
	530	& llvar1, llvar2, &
	531	& jpi, jpj, jpk, &
	532	& zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, &
	533	& ln_nea, ln_bound_reject, &
	534	& kdailyavtypes = nn_profdavtypes )
[2128]	535
[7992]	536	END DO
[2128]	537
[7992]	538	DEALLOCATE( ifilesprof, clproffiles )
[2128]	539
	540	ENDIF
	541
[7992]	542	IF ( nsurftypes > 0 ) THEN
[2128]	543
[7992]	544	ALLOCATE(surfdata(nsurftypes))
	545	ALLOCATE(surfdataqc(nsurftypes))
	546	ALLOCATE(nvarssurf(nsurftypes))
	547	ALLOCATE(nextrsurf(nsurftypes))
[2128]	548
[7992]	549	DO jtype = 1, nsurftypes
[2128]	550
[7992]	551	nvarssurf(jtype) = 1
	552	nextrsurf(jtype) = 0
	553	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
[2128]	554
[7992]	555	!Read in surface obs types
	556	CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
	557	& clsurffiles(jtype,1:ifilessurf(jtype)), &
	558	& nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
	559	& rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
[2128]	560
[7992]	561	CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
[2128]	562
[7992]	563	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
	564	CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
	565	IF ( ln_altbias ) &
	566	& CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
	567	ENDIF
[2128]	568
[7992]	569	IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
	570	jnumsstbias = 0
	571	DO jfile = 1, jpmaxnfiles
	572	IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
	573	& jnumsstbias = jnumsstbias + 1
	574	END DO
	575	IF ( jnumsstbias == 0 ) THEN
	576	CALL ctl_stop("ln_sstbias set but no bias files to read in")
	577	ENDIF
[2128]	578
[7992]	579	CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &
	580	& jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )
[2128]	581
[7992]	582	ENDIF
[2128]	583
[7992]	584	END DO
[2128]	585
[7992]	586	DEALLOCATE( ifilessurf, clsurffiles )
[2128]	587
[7992]	588	ENDIF
[2128]	589
[7992]	590	CALL wrk_dealloc( jpi, jpj, zglam1 )
	591	CALL wrk_dealloc( jpi, jpj, zglam2 )
	592	CALL wrk_dealloc( jpi, jpj, zgphi1 )
	593	CALL wrk_dealloc( jpi, jpj, zgphi2 )
	594	CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
	595	CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
[2128]	596
	597	END SUBROUTINE dia_obs_init
	598
	599	SUBROUTINE dia_obs( kstp )
	600	!!----------------------------------------------------------------------
	601	!! * ROUTINE dia_obs *
	602	!!
	603	!! ** Purpose : Call the observation operators on each time step
	604	!!
	605	!! ** Method : Call the observation operators on each time step to
[7992]	606	!! compute the model equivalent of the following data:
	607	!! - Profile data, currently T/S or U/V
	608	!! - Surface data, currently SST, SLA or sea-ice concentration.
[2128]	609	!!
[7992]	610	!! ** Action :
[2128]	611	!!
	612	!! History :
	613	!! ! 06-03 (K. Mogensen) Original code
	614	!! ! 06-05 (K. Mogensen) Reformatted
	615	!! ! 06-10 (A. Weaver) Cleaning
	616	!! ! 07-03 (K. Mogensen) General handling of profiles
	617	!! ! 07-04 (G. Smith) Generalized surface operators
	618	!! ! 08-10 (M. Valdivieso) obs operator for velocity profiles
[7992]	619	!! ! 15-08 (M. Martin) Combined surface/profile routines.
[2128]	620	!!----------------------------------------------------------------------
	621	!! * Modules used
[7992]	622	USE phycst, ONLY : & ! Physical constants
	623	& rday
	624	USE oce, ONLY : & ! Ocean dynamics and tracers variables
	625	& tsn, &
	626	& un, &
	627	& vn, &
[2128]	628	& sshn
[3294]	629	#if defined key_lim3
[7992]	630	USE ice, ONLY : & ! LIM3 Ice model variables
[2128]	631	& frld
	632	#endif
[3294]	633	#if defined key_lim2
[7992]	634	USE ice_2, ONLY : & ! LIM2 Ice model variables
[3294]	635	& frld
[2715]	636	#endif
[7992]	637	#if defined key_cice
	638	USE sbc_oce, ONLY : fr_i ! ice fraction
	639	#endif
	640	#if defined key_hadocc
	641	USE trc, ONLY : & ! HadOCC chlorophyll, fCO2 and pCO2
	642	& HADOCC_CHL, &
	643	& HADOCC_FCO2, &
	644	& HADOCC_PCO2, &
	645	& HADOCC_FILL_FLT
	646	#elif defined key_medusa && defined key_foam_medusa
	647	USE trc, ONLY : & ! MEDUSA chlorophyll, fCO2 and pCO2
[8455]	648	& trn
	649	USE par_medusa, ONLY: &
	650	& jpchn, &
	651	& jpchd
	652	#if defined key_roam
	653	USE sms_medusa, ONLY: &
	654	& f2_pco2w, &
	655	& f2_fco2w
	656	#endif
[7992]	657	#elif defined key_fabm
	658	USE fabm
	659	USE par_fabm
	660	#endif
	661	#if defined key_spm
	662	USE par_spm, ONLY: & ! ERSEM/SPM sediments
	663	& jp_spm
	664	USE trc, ONLY : &
	665	& trn
	666	#endif
	667
[2128]	668	IMPLICIT NONE
	669
	670	!! * Arguments
[7992]	671	INTEGER, INTENT(IN) :: kstp ! Current timestep
[2128]	672	!! * Local declarations
[7992]	673	INTEGER :: idaystp ! Number of timesteps per day
	674	INTEGER :: jtype ! Data loop variable
	675	INTEGER :: jvar ! Variable number
	676	INTEGER :: ji, jj ! Loop counters
	677	REAL(wp) :: tiny ! small number
	678	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
	679	& zprofvar1, & ! Model values for 1st variable in a prof ob
	680	& zprofvar2 ! Model values for 2nd variable in a prof ob
	681	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
	682	& zprofmask1, & ! Mask associated with zprofvar1
	683	& zprofmask2 ! Mask associated with zprofvar2
	684	REAL(wp), POINTER, DIMENSION(:,:) :: &
	685	& zsurfvar, & ! Model values equivalent to surface ob.
	686	& zsurfmask ! Mask associated with surface variable
	687	REAL(wp), POINTER, DIMENSION(:,:) :: &
	688	& zglam1, & ! Model longitudes for prof variable 1
	689	& zglam2, & ! Model longitudes for prof variable 2
	690	& zgphi1, & ! Model latitudes for prof variable 1
	691	& zgphi2 ! Model latitudes for prof variable 2
[2715]	692
[7992]	693
	694	!Allocate local work arrays
	695	CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
	696	CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
	697	CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
	698	CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
	699	CALL wrk_alloc( jpi, jpj, zsurfvar )
	700	CALL wrk_alloc( jpi, jpj, zsurfmask )
	701	CALL wrk_alloc( jpi, jpj, zglam1 )
	702	CALL wrk_alloc( jpi, jpj, zglam2 )
	703	CALL wrk_alloc( jpi, jpj, zgphi1 )
	704	CALL wrk_alloc( jpi, jpj, zgphi2 )
	705
[2128]	706	IF(lwp) THEN
	707	WRITE(numout,*)
	708	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
	709	WRITE(numout,*) '~~~~~~~'
[7992]	710	CALL FLUSH(numout)
[2128]	711	ENDIF
	712
	713	idaystp = NINT( rday / rdt )
	714
	715	!-----------------------------------------------------------------------
[7992]	716	! Call the profile and surface observation operators
[2128]	717	!-----------------------------------------------------------------------
	718
[7992]	719	IF ( nproftypes > 0 ) THEN
[2128]	720
[7992]	721	DO jtype = 1, nproftypes
[2128]	722
[7992]	723	SELECT CASE ( TRIM(cobstypesprof(jtype)) )
	724	CASE('prof')
	725	zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
	726	zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
	727	zprofmask1(:,:,:) = tmask(:,:,:)
	728	zprofmask2(:,:,:) = tmask(:,:,:)
	729	zglam1(:,:) = glamt(:,:)
	730	zglam2(:,:) = glamt(:,:)
	731	zgphi1(:,:) = gphit(:,:)
	732	zgphi2(:,:) = gphit(:,:)
	733	CASE('vel')
	734	zprofvar1(:,:,:) = un(:,:,:)
	735	zprofvar2(:,:,:) = vn(:,:,:)
	736	zprofmask1(:,:,:) = umask(:,:,:)
	737	zprofmask2(:,:,:) = vmask(:,:,:)
	738	zglam1(:,:) = glamu(:,:)
	739	zglam2(:,:) = glamv(:,:)
	740	zgphi1(:,:) = gphiu(:,:)
	741	zgphi2(:,:) = gphiv(:,:)
	742	CASE DEFAULT
	743	CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
	744	END SELECT
	745
	746	CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, &
	747	& nit000, idaystp, &
	748	& zprofvar1, zprofvar2, &
	749	& fsdept(:,:,:), fsdepw(:,:,:), &
	750	& zprofmask1, zprofmask2, &
	751	& zglam1, zglam2, zgphi1, zgphi2, &
	752	& nn_1dint, nn_2dint, &
	753	& kdailyavtypes = nn_profdavtypes )
	754
[2128]	755	END DO
	756
	757	ENDIF
	758
[7992]	759	IF ( nsurftypes > 0 ) THEN
	760
	761	DO jtype = 1, nsurftypes
	762
	763	!Defaults which might be changed
	764	zsurfmask(:,:) = tmask(:,:,1)
	765
	766	SELECT CASE ( TRIM(cobstypessurf(jtype)) )
	767	CASE('sst')
	768	zsurfvar(:,:) = tsn(:,:,1,jp_tem)
	769	CASE('sla')
	770	zsurfvar(:,:) = sshn(:,:)
	771	CASE('sss')
	772	zsurfvar(:,:) = tsn(:,:,1,jp_sal)
	773	CASE('sic')
	774	IF ( kstp == 0 ) THEN
	775	IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
	776	CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
	777	& 'time-step but some obs are valid then.' )
	778	WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
	779	& ' sea-ice obs will be missed'
	780	ENDIF
	781	surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
	782	& surfdataqc(jtype)%nsstp(1)
	783	CYCLE
	784	ELSE
	785	#if defined key_cice
	786	zsurfvar(:,:) = fr_i(:,:)
	787	#elif defined key_lim2 \|\| defined key_lim3
	788	zsurfvar(:,:) = 1._wp - frld(:,:)
	789	#else
	790	CALL ctl_stop( ' Trying to run sea-ice observation operator', &
	791	& ' but no sea-ice model appears to have been defined' )
[2128]	792	#endif
[7992]	793	ENDIF
[2128]	794
[7992]	795	CASE('logchl')
	796	#if defined key_hadocc
	797	zsurfvar(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
	798	#elif defined key_medusa && defined key_foam_medusa
[8646]	799	! Add non-diatom and diatom surface chlorophyll from MEDUSA
[8455]	800	zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
[7992]	801	#elif defined key_fabm
	802	chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
	803	zsurfvar(:,:) = chl_3d(:,:,1)
	804	#else
	805	CALL ctl_stop( ' Trying to run logchl observation operator', &
	806	& ' but no biogeochemical model appears to have been defined' )
	807	#endif
	808	zsurfmask(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
	809	! Take the log10 where we can, otherwise exclude
	810	tiny = 1.0e-20
	811	WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
	812	zsurfvar(:,:) = LOG10(zsurfvar(:,:))
	813	ELSEWHERE
	814	zsurfvar(:,:) = obfillflt
	815	zsurfmask(:,:) = 0
	816	END WHERE
	817	CASE('spm')
	818	#if defined key_spm
	819	zsurfvar(:,:) = 0.0
	820	DO jn = 1, jp_spm
	821	zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes
	822	END DO
	823	#else
	824	CALL ctl_stop( ' Trying to run spm observation operator', &
	825	& ' but no spm model appears to have been defined' )
	826	#endif
	827	CASE('fco2')
	828	#if defined key_hadocc
	829	zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC
	830	IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
	831	& ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
	832	zsurfvar(:,:) = obfillflt
	833	zsurfmask(:,:) = 0
	834	CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
	835	& ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
	836	ENDIF
[8455]	837	#elif defined key_medusa && defined key_foam_medusa && defined key_roam
	838	zsurfvar(:,:) = f2_fco2w(:,:)
[7992]	839	#elif defined key_fabm
	840	! First, get pCO2 from FABM
	841	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
	842	zsurfvar(:,:) = pco2_3d(:,:,1)
	843	! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
	844	! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
	845	! and data reduction routines, Deep-Sea Research II, 56: 512-522.
	846	! and
	847	! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
	848	! Marine Chemistry, 2: 203-215.
	849	! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
	850	! not explicitly included - atmospheric pressure is not necessarily available so this is
	851	! the best assumption.
	852	! Further, the (1-xCO2)^2 term has been neglected. This is common practice
	853	! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
	854	! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
	855	! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
	856	zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + &
	857	& 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - &
	858	& 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
	859	& 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0) + &
	860	& 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / &
	861	& (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
	862	#else
	863	CALL ctl_stop( ' Trying to run fco2 observation operator', &
	864	& ' but no biogeochemical model appears to have been defined' )
	865	#endif
	866	CASE('pco2')
	867	#if defined key_hadocc
	868	zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC
	869	IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
	870	& ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
	871	zsurfvar(:,:) = obfillflt
	872	zsurfmask(:,:) = 0
	873	CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
	874	& ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
	875	ENDIF
[8455]	876	#elif defined key_medusa && defined key_foam_medusa && defined key_roam
	877	zsurfvar(:,:) = f2_pco2w(:,:)
[7992]	878	#elif defined key_fabm
	879	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
	880	zsurfvar(:,:) = pco2_3d(:,:,1)
	881	#else
	882	CALL ctl_stop( ' Trying to run pCO2 observation operator', &
	883	& ' but no biogeochemical model appears to have been defined' )
	884	#endif
	885
	886	CASE DEFAULT
	887
	888	CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
	889
	890	END SELECT
	891
	892	CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, &
	893	& nit000, idaystp, zsurfvar, zsurfmask, &
	894	& n2dintsurf(jtype), llnightav(jtype), &
	895	& ravglamscl(jtype), ravgphiscl(jtype), &
	896	& lfpindegs(jtype) )
	897
[2128]	898	END DO
[7992]	899
[2128]	900	ENDIF
	901
[7992]	902	CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
	903	CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
	904	CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
	905	CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
	906	CALL wrk_dealloc( jpi, jpj, zsurfvar )
	907	CALL wrk_dealloc( jpi, jpj, zsurfmask )
	908	CALL wrk_dealloc( jpi, jpj, zglam1 )
	909	CALL wrk_dealloc( jpi, jpj, zglam2 )
	910	CALL wrk_dealloc( jpi, jpj, zgphi1 )
	911	CALL wrk_dealloc( jpi, jpj, zgphi2 )
[2715]	912
[2128]	913	END SUBROUTINE dia_obs
[7992]	914
	915	SUBROUTINE dia_obs_wri
[2128]	916	!!----------------------------------------------------------------------
	917	!! * ROUTINE dia_obs_wri *
	918	!!
	919	!! ** Purpose : Call observation diagnostic output routines
	920	!!
	921	!! ** Method : Call observation diagnostic output routines
	922	!!
[7992]	923	!! ** Action :
[2128]	924	!!
	925	!! History :
	926	!! ! 06-03 (K. Mogensen) Original code
	927	!! ! 06-05 (K. Mogensen) Reformatted
	928	!! ! 06-10 (A. Weaver) Cleaning
	929	!! ! 07-03 (K. Mogensen) General handling of profiles
	930	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
[7992]	931	!! ! 15-08 (M. Martin) Combined writing for prof and surf types
[2128]	932	!!----------------------------------------------------------------------
[7992]	933	!! * Modules used
	934	USE obs_rot_vel ! Rotation of velocities
	935
[2128]	936	IMPLICIT NONE
	937
	938	!! * Local declarations
[7992]	939	INTEGER :: jtype ! Data set loop variable
	940	INTEGER :: jo, jvar, jk
	941	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	942	& zu, &
	943	& zv
[2128]	944
	945	!-----------------------------------------------------------------------
	946	! Depending on switches call various observation output routines
	947	!-----------------------------------------------------------------------
	948
[7992]	949	IF ( nproftypes > 0 ) THEN
[2128]	950
[7992]	951	DO jtype = 1, nproftypes
[2128]	952
[7992]	953	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
[2128]	954
[7992]	955	! For velocity data, rotate the model velocities to N/S, E/W
	956	! using the compressed data structure.
	957	ALLOCATE( &
	958	& zu(profdataqc(jtype)%nvprot(1)), &
	959	& zv(profdataqc(jtype)%nvprot(2)) &
	960	& )
[2128]	961
[7992]	962	CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
[2128]	963
[7992]	964	DO jo = 1, profdataqc(jtype)%nprof
	965	DO jvar = 1, 2
	966	DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
[2128]	967
[7992]	968	IF ( jvar == 1 ) THEN
	969	profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
	970	ELSE
	971	profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
	972	ENDIF
[2128]	973
[7992]	974	END DO
	975	END DO
	976	END DO
[2128]	977
[7992]	978	DEALLOCATE( zu )
	979	DEALLOCATE( zv )
[2128]	980
[7992]	981	END IF
[2128]	982
[7992]	983	CALL obs_prof_decompress( profdataqc(jtype), &
	984	& profdata(jtype), .TRUE., numout )
[2128]	985
[7992]	986	CALL obs_wri_prof( profdata(jtype) )
[2128]	987
	988	END DO
	989
	990	ENDIF
	991
[7992]	992	IF ( nsurftypes > 0 ) THEN
[2128]	993
[7992]	994	DO jtype = 1, nsurftypes
[2128]	995
[7992]	996	CALL obs_surf_decompress( surfdataqc(jtype), &
	997	& surfdata(jtype), .TRUE., numout )
[2128]	998
[7992]	999	CALL obs_wri_surf( surfdata(jtype) )
[2128]	1000
	1001	END DO
	1002
	1003	ENDIF
	1004
	1005	END SUBROUTINE dia_obs_wri
	1006
[4245]	1007	SUBROUTINE dia_obs_dealloc
	1008	IMPLICIT NONE
	1009	!!----------------------------------------------------------------------
	1010	!! * ROUTINE dia_obs_dealloc *
	1011	!!
	1012	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
	1013	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
	1014	!!
	1015	!! ** Method : Clean up various arrays left behind by the obs_oper.
	1016	!!
	1017	!! ** Action :
	1018	!!
	1019	!!----------------------------------------------------------------------
[7992]	1020	! obs_grid deallocation
[4245]	1021	CALL obs_grid_deallocate
	1022
[7992]	1023	! diaobs deallocation
	1024	IF ( nproftypes > 0 ) &
	1025	& DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
	1026
	1027	IF ( nsurftypes > 0 ) &
	1028	& DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
	1029	& n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
	1030
[4245]	1031	END SUBROUTINE dia_obs_dealloc
	1032
[2128]	1033	SUBROUTINE ini_date( ddobsini )
	1034	!!----------------------------------------------------------------------
	1035	!! * ROUTINE ini_date *
	1036	!!
[7992]	1037	!! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
[2128]	1038	!!
[7992]	1039	!! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format
[2128]	1040	!!
[7992]	1041	!! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format
	1042	!!
[2128]	1043	!! History :
	1044	!! ! 06-03 (K. Mogensen) Original code
	1045	!! ! 06-05 (K. Mogensen) Reformatted
	1046	!! ! 06-10 (A. Weaver) Cleaning
	1047	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
	1048	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
	1049	!!----------------------------------------------------------------------
	1050	USE phycst, ONLY : & ! Physical constants
	1051	& rday
	1052	USE dom_oce, ONLY : & ! Ocean space and time domain variables
	1053	& rdt
	1054
	1055	IMPLICIT NONE
	1056
	1057	!! * Arguments
[7992]	1058	REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
[2128]	1059
	1060	!! * Local declarations
	1061	INTEGER :: iyea ! date - (year, month, day, hour, minute)
	1062	INTEGER :: imon
	1063	INTEGER :: iday
	1064	INTEGER :: ihou
	1065	INTEGER :: imin
[7992]	1066	INTEGER :: imday ! Number of days in month.
	1067	INTEGER, DIMENSION(12) :: &
	1068	& imonth_len ! Length in days of the months of the current year
	1069	REAL(wp) :: zdayfrc ! Fraction of day
[2128]	1070
[7992]	1071	!----------------------------------------------------------------------
	1072	! Initial date initialization (year, month, day, hour, minute)
	1073	! (This assumes that the initial date is for 00z))
	1074	!----------------------------------------------------------------------
[2128]	1075	iyea = ndate0 / 10000
	1076	imon = ( ndate0 - iyea * 10000 ) / 100
	1077	iday = ndate0 - iyea * 10000 - imon * 100
	1078	ihou = 0
	1079	imin = 0
	1080
[7992]	1081	!----------------------------------------------------------------------
	1082	! Compute number of days + number of hours + min since initial time
	1083	!----------------------------------------------------------------------
[2128]	1084	iday = iday + ( nit000 -1 ) * rdt / rday
	1085	zdayfrc = ( nit000 -1 ) * rdt / rday
	1086	zdayfrc = zdayfrc - aint(zdayfrc)
	1087	ihou = int( zdayfrc * 24 )
	1088	imin = int( (zdayfrc * 24 - ihou) * 60 )
	1089
[7992]	1090	!-----------------------------------------------------------------------
	1091	! Convert number of days (iday) into a real date
	1092	!----------------------------------------------------------------------
[2128]	1093
	1094	CALL calc_month_len( iyea, imonth_len )
[7992]	1095
[2128]	1096	DO WHILE ( iday > imonth_len(imon) )
	1097	iday = iday - imonth_len(imon)
	1098	imon = imon + 1
	1099	IF ( imon > 12 ) THEN
	1100	imon = 1
	1101	iyea = iyea + 1
	1102	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
	1103	ENDIF
	1104	END DO
	1105
[7992]	1106	!----------------------------------------------------------------------
	1107	! Convert it into YYYYMMDD.HHMMSS format.
	1108	!----------------------------------------------------------------------
[2128]	1109	ddobsini = iyea * 10000_dp + imon * 100_dp + &
	1110	& iday + ihou * 0.01_dp + imin * 0.0001_dp
	1111
	1112
	1113	END SUBROUTINE ini_date
	1114
	1115	SUBROUTINE fin_date( ddobsfin )
	1116	!!----------------------------------------------------------------------
	1117	!! * ROUTINE fin_date *
	1118	!!
[7992]	1119	!! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
[2128]	1120	!!
[7992]	1121	!! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format
[2128]	1122	!!
[7992]	1123	!! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format
	1124	!!
[2128]	1125	!! History :
	1126	!! ! 06-03 (K. Mogensen) Original code
	1127	!! ! 06-05 (K. Mogensen) Reformatted
	1128	!! ! 06-10 (A. Weaver) Cleaning
	1129	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
	1130	!!----------------------------------------------------------------------
	1131	USE phycst, ONLY : & ! Physical constants
	1132	& rday
	1133	USE dom_oce, ONLY : & ! Ocean space and time domain variables
	1134	& rdt
	1135
	1136	IMPLICIT NONE
	1137
	1138	!! * Arguments
[7992]	1139	REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
[2128]	1140
	1141	!! * Local declarations
	1142	INTEGER :: iyea ! date - (year, month, day, hour, minute)
	1143	INTEGER :: imon
	1144	INTEGER :: iday
	1145	INTEGER :: ihou
	1146	INTEGER :: imin
[7992]	1147	INTEGER :: imday ! Number of days in month.
	1148	INTEGER, DIMENSION(12) :: &
	1149	& imonth_len ! Length in days of the months of the current year
	1150	REAL(wp) :: zdayfrc ! Fraction of day
	1151
[2128]	1152	!-----------------------------------------------------------------------
	1153	! Initial date initialization (year, month, day, hour, minute)
	1154	! (This assumes that the initial date is for 00z)
	1155	!-----------------------------------------------------------------------
	1156	iyea = ndate0 / 10000
	1157	imon = ( ndate0 - iyea * 10000 ) / 100
	1158	iday = ndate0 - iyea * 10000 - imon * 100
	1159	ihou = 0
	1160	imin = 0
[7992]	1161
[2128]	1162	!-----------------------------------------------------------------------
	1163	! Compute number of days + number of hours + min since initial time
	1164	!-----------------------------------------------------------------------
	1165	iday = iday + nitend * rdt / rday
	1166	zdayfrc = nitend * rdt / rday
	1167	zdayfrc = zdayfrc - AINT( zdayfrc )
	1168	ihou = INT( zdayfrc * 24 )
	1169	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
	1170
	1171	!-----------------------------------------------------------------------
	1172	! Convert number of days (iday) into a real date
	1173	!----------------------------------------------------------------------
	1174
	1175	CALL calc_month_len( iyea, imonth_len )
[7992]	1176
[2128]	1177	DO WHILE ( iday > imonth_len(imon) )
	1178	iday = iday - imonth_len(imon)
	1179	imon = imon + 1
	1180	IF ( imon > 12 ) THEN
	1181	imon = 1
	1182	iyea = iyea + 1
	1183	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
	1184	ENDIF
	1185	END DO
	1186
	1187	!-----------------------------------------------------------------------
	1188	! Convert it into YYYYMMDD.HHMMSS format
	1189	!-----------------------------------------------------------------------
	1190	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
	1191	& + ihou * 0.01_dp + imin * 0.0001_dp
	1192
	1193	END SUBROUTINE fin_date
[7992]	1194
	1195	SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
	1196	& cfilestype, ifiles, cobstypes, cfiles )
	1197
	1198	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
	1199	INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
	1200	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
	1201	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
	1202	& ifiles ! Out appended number of files for this type
	1203
	1204	CHARACTER(len=6), INTENT(IN) :: ctypein
	1205	CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
	1206	& cfilestype ! In list of files for this obs type
	1207	CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
	1208	& cobstypes ! Out appended list of obs types
	1209	CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
	1210	& cfiles ! Out appended list of files for all types
	1211
	1212	!Local variables
	1213	INTEGER :: jfile
	1214
	1215	cfiles(jtype,:) = cfilestype(:)
	1216	cobstypes(jtype) = ctypein
	1217	ifiles(jtype) = 0
	1218	DO jfile = 1, jpmaxnfiles
	1219	IF ( trim(cfiles(jtype,jfile)) /= '' ) &
	1220	ifiles(jtype) = ifiles(jtype) + 1
	1221	END DO
	1222
	1223	IF ( ifiles(jtype) == 0 ) THEN
	1224	CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// &
	1225	& ' set to true but no files available to read' )
	1226	ENDIF
	1227
	1228	IF(lwp) THEN
	1229	WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:'
	1230	DO jfile = 1, ifiles(jtype)
	1231	WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile))
	1232	END DO
	1233	ENDIF
	1234
	1235	END SUBROUTINE obs_settypefiles
	1236
	1237	SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, &
	1238	& n2dint_default, n2dint_type, &
	1239	& ravglamscl_type, ravgphiscl_type, &
	1240	& lfp_indegs_type, lavnight_type, &
	1241	& n2dint, ravglamscl, ravgphiscl, &
	1242	& lfpindegs, lavnight )
	1243
	1244	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
	1245	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
	1246	INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type
	1247	INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type
	1248	REAL(wp), INTENT(IN) :: &
	1249	& ravglamscl_type, & !E/W diameter of obs footprint for this type
	1250	& ravgphiscl_type !N/S diameter of obs footprint for this type
	1251	LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres
	1252	LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average
	1253	CHARACTER(len=6), INTENT(IN) :: ctypein
	1254
	1255	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
	1256	& n2dint
	1257	REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
	1258	& ravglamscl, ravgphiscl
	1259	LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
	1260	& lfpindegs, lavnight
	1261
	1262	lavnight(jtype) = lavnight_type
	1263
	1264	IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
	1265	n2dint(jtype) = n2dint_type
	1266	ELSE
	1267	n2dint(jtype) = n2dint_default
	1268	ENDIF
	1269
	1270	! For averaging observation footprints set options for size of footprint
	1271	IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
	1272	IF ( ravglamscl_type > 0._wp ) THEN
	1273	ravglamscl(jtype) = ravglamscl_type
	1274	ELSE
	1275	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
	1276	'scale (ravglamscl) for observation type '//TRIM(ctypein) )
	1277	ENDIF
	1278
	1279	IF ( ravgphiscl_type > 0._wp ) THEN
	1280	ravgphiscl(jtype) = ravgphiscl_type
	1281	ELSE
	1282	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
	1283	'scale (ravgphiscl) for observation type '//TRIM(ctypein) )
	1284	ENDIF
	1285
	1286	lfpindegs(jtype) = lfp_indegs_type
	1287
	1288	ENDIF
	1289
	1290	! Write out info
	1291	IF(lwp) THEN
	1292	IF ( n2dint(jtype) <= 4 ) THEN
	1293	WRITE(numout,*) ' '//TRIM(ctypein)// &
	1294	& ' model counterparts will be interpolated horizontally'
	1295	ELSE IF ( n2dint(jtype) <= 6 ) THEN
	1296	WRITE(numout,*) ' '//TRIM(ctypein)// &
	1297	& ' model counterparts will be averaged horizontally'
	1298	WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype)
	1299	WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype)
	1300	IF ( lfpindegs(jtype) ) THEN
	1301	WRITE(numout,*) ' '//' (in degrees)'
	1302	ELSE
	1303	WRITE(numout,*) ' '//' (in metres)'
	1304	ENDIF
	1305	ENDIF
	1306	ENDIF
	1307
	1308	END SUBROUTINE obs_setinterpopts
	1309
[2128]	1310	END MODULE diaobs

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

Download in other formats: