New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaobs.F90 in branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

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

Last change on this file since 9016 was 9016, checked in by dford, 6 years ago
Change naming convention and fix a couple of bugs - code now compiles and runs.
File size: 74.8 KB

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

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

Download in other formats: