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

Context Navigation

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_pfts/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8105

Last change on this file since 8105 was 8105, checked in by dford, 7 years ago
Implement initial version of PFT logchl obs operator, developed by Jozef S.
File size: 99.5 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	!! 'key_diaobs' : Switch on the observation diagnostic computation
10	!!----------------------------------------------------------------------
11	!! dia_obs_init : Reading and prepare observations
12	!! dia_obs : Compute model equivalent to observations
13	!! dia_obs_wri : Write observational diagnostics
14	!! ini_date : Compute the initial date YYYYMMDD.HHMMSS
15	!! fin_date : Compute the final date YYYYMMDD.HHMMSS
16	!!----------------------------------------------------------------------
17	!! * Modules used
18	USE trc, ONLY: trn ! ERSEM state variables
19	USE wrk_nemo ! Memory Allocation
20	USE par_kind ! Precision variables
21	USE in_out_manager ! I/O manager
22	USE par_oce
23	USE dom_oce ! Ocean space and time domain variables
24	USE obs_const, ONLY: obfillflt ! Fill value
25	USE obs_fbm, ONLY: ln_cl4 ! Class 4 diagnostic switch
26	USE obs_read_prof ! Reading and allocation of observations (Coriolis)
27	USE obs_read_sla ! Reading and allocation of SLA observations
28	USE obs_read_sst ! Reading and allocation of SST observations
29	USE obs_sstbias ! Bias correction routine for SST
30	USE obs_readmdt ! Reading and allocation of MDT for SLA.
31	USE obs_read_seaice ! Reading and allocation of Sea Ice observations
32	USE obs_read_vel ! Reading and allocation of velocity component observations
33	USE obs_read_logchl ! Reading and allocation of logchl observations
34	USE obs_read_spm ! Reading and allocation of spm observations
35	USE obs_read_fco2 ! Reading and allocation of fco2 observations
36	USE obs_read_pco2 ! Reading and allocation of pco2 observations
37	USE obs_prep ! Preparation of obs. (grid search etc).
38	USE obs_oper ! Observation operators
39	USE obs_write ! Writing of observation related diagnostics
40	USE obs_grid ! Grid searching
41	USE obs_read_altbias ! Bias treatment for altimeter
42	USE obs_profiles_def ! Profile data definitions
43	USE obs_profiles ! Profile data storage
44	USE obs_surf_def ! Surface data definitions
45	USE obs_sla ! SLA data storage
46	USE obs_sst ! SST data storage
47	USE obs_seaice ! Sea Ice data storage
48	USE obs_logchl ! logchl data storage
49	USE obs_spm ! spm data storage
50	USE obs_fco2 ! fco2 data storage
51	USE obs_pco2 ! pco2 data storage
52	USE obs_types ! Definitions for observation types
53	USE mpp_map ! MPP mapping
54	USE lib_mpp ! For ctl_warn/stop
55
56	IMPLICIT NONE
57
58	!! * Routine accessibility
59	PRIVATE
60	PUBLIC dia_obs_init, & ! Initialize and read observations
61	& dia_obs, & ! Compute model equivalent to observations
62	& dia_obs_wri, & ! Write model equivalent to observations
63	& dia_obs_dealloc ! Deallocate dia_obs data
64
65	!! * Shared Module variables
66	LOGICAL, PUBLIC, PARAMETER :: &
67	#if defined key_diaobs
68	& lk_diaobs = .TRUE. !: Logical switch for observation diangostics
69	#else
70	& lk_diaobs = .FALSE. !: Logical switch for observation diangostics
71	#endif
72
73	!! * Module variables
74	LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles
75	LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles
76	LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set
77	LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set
78	LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles
79	LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies
80	LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files
81	LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files
82	LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature
83	LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature
84	LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data
85	LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files
86	LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration
87	LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations
88	LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data
89	LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data
90	LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data
91	LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data
92	LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files
93	LOGICAL, PUBLIC :: ln_logchl !: Logical switch for log10(chlorophyll)
94	LOGICAL, PUBLIC :: ln_logchlfb !: Logical switch for logchl from feedback files
95
96
97
98
99
100
101
102
103
104
105
106	! JOZEF ADDITION
107
108
109	LOGICAL, PUBLIC :: ln_logchlpft !: Logical switch for log10(chlorophyll) PFTs
110	LOGICAL, PUBLIC :: ln_logchlpftfb !: Logical switch for logchl PFTs from feedback files
111
112	INTEGER, PUBLIC :: nn_logchlpftscc = 4 ! Number of logchl PFTs
113
114
115
116	! END
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131	LOGICAL, PUBLIC :: ln_spm !: Logical switch for spm
132	LOGICAL, PUBLIC :: ln_spmfb !: Logical switch for spm from feedback files
133	LOGICAL, PUBLIC :: ln_fco2 !: Logical switch for fco2
134	LOGICAL, PUBLIC :: ln_fco2fb !: Logical switch for fco2 from feedback files
135	LOGICAL, PUBLIC :: ln_pco2 !: Logical switch for pco2
136	LOGICAL, PUBLIC :: ln_pco2fb !: Logical switch for pco2 from feedback files
137	LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height
138	LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity
139	LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations
140	LOGICAL, PUBLIC :: ln_nea !: Remove observations near land
141	LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias
142	LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files
143	LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations
144	LOGICAL, PUBLIC :: ln_sstbias !: Logical switch for bias corection of SST
145
146	REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS
147	REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS
148
149	INTEGER, PUBLIC :: n1dint !: Vertical interpolation method
150	INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method
151
152	INTEGER, DIMENSION(imaxavtypes) :: &
153	& endailyavtypes !: ENACT data types which are daily average
154
155	INTEGER, PARAMETER :: MaxNumFiles = 1000
156
157
158
159
160	LOGICAL, DIMENSION(MaxNumFiles) :: &
161	& ln_profb_ena, & !: Is the feedback files from ENACT data ?
162	! !: If so use endailyavtypes
163	& ln_profb_enatim !: Change tim for 820 enact data set.
164
165	INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type
166
167	LOGICAL, DIMENSION(MaxNumFiles) :: &
168	& ln_velfb_av !: Is the velocity feedback files daily average?
169	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
170	& ld_enact !: Profile data is ENACT so use endailyavtypes
171	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
172	& ld_velav !: Velocity data is daily averaged
173	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
174	& ld_sstnight !: SST observation corresponds to night mean
175
176	!!----------------------------------------------------------------------
177	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
178	!! $Id$
179	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
180	!!----------------------------------------------------------------------
181
182	!! * Substitutions
183	# include "domzgr_substitute.h90"
184	CONTAINS
185
186	SUBROUTINE dia_obs_init
187	!!----------------------------------------------------------------------
188	!! * ROUTINE dia_obs_init *
189	!!
190	!! ** Purpose : Initialize and read observations
191	!!
192	!! ** Method : Read the namelist and call reading routines
193	!!
194	!! ** Action : Read the namelist and call reading routines
195	!!
196	!! History :
197	!! ! 06-03 (K. Mogensen) Original code
198	!! ! 06-05 (A. Weaver) Reformatted
199	!! ! 06-10 (A. Weaver) Cleaning and add controls
200	!! ! 07-03 (K. Mogensen) General handling of profiles
201	!! ! 14-08 (J.While) Incorporated SST bias correction
202	!!----------------------------------------------------------------------
203
204	IMPLICIT NONE
205
206
207	!! * Local declarations
208	CHARACTER(len=128) :: enactfiles(MaxNumFiles)
209	CHARACTER(len=128) :: coriofiles(MaxNumFiles)
210	CHARACTER(len=128) :: profbfiles(MaxNumFiles)
211	CHARACTER(len=128) :: sstfiles(MaxNumFiles)
212	CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
213	CHARACTER(len=128) :: sstbias_files(MaxNumFiles)
214	CHARACTER(len=128) :: slafilesact(MaxNumFiles)
215	CHARACTER(len=128) :: slafilespas(MaxNumFiles)
216	CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
217	CHARACTER(len=128) :: seaicefiles(MaxNumFiles)
218	CHARACTER(len=128) :: velcurfiles(MaxNumFiles)
219	CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)
220	CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
221	CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
222	CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
223	CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
224	CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
225	CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
226	CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
227
228
229
230
231
232	! JOZEF ADDITION
233
234
235
236	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: logchlpftfiles
237	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: logchlpftfbfiles
238
239	! CHARACTER(len=128) :: logchlpftfiles(MaxNumFiles, nn_logchlpftscc)
240	! CHARACTER(len=128) :: logchlpftfbfiles(MaxNumFiles, nn_logchlpftscc)
241
242	! END
243
244
245
246
247
248
249
250	CHARACTER(len=128) :: spmfiles(MaxNumFiles)
251	CHARACTER(len=128) :: spmfbfiles(MaxNumFiles)
252	CHARACTER(len=128) :: fco2files(MaxNumFiles)
253	CHARACTER(len=128) :: fco2fbfiles(MaxNumFiles)
254	CHARACTER(len=128) :: pco2files(MaxNumFiles)
255	CHARACTER(len=128) :: pco2fbfiles(MaxNumFiles)
256	CHARACTER(LEN=128) :: reysstname
257	CHARACTER(LEN=12) :: reysstfmt
258	CHARACTER(LEN=128) :: bias_file
259	CHARACTER(LEN=20) :: datestr=" ", timestr=" "
260	NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, &
261	& ln_sla, ln_sladt, ln_slafb, &
262	& ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, &
263	& ln_bound_reject, &
264	& enactfiles, coriofiles, profbfiles, &
265	& slafilesact, slafilespas, slafbfiles, &
266	& sstfiles, sstfbfiles, &
267	& ln_seaice, seaicefiles, &
268	& dobsini, dobsend, n1dint, n2dint, &
269	& nmsshc, mdtcorr, mdtcutoff, &
270	& ln_reysst, ln_ghrsst, reysstname, reysstfmt, &
271	& ln_sstnight, &
272	& ln_grid_search_lookup, &
273	& grid_search_file, grid_search_res, &
274	& ln_grid_global, bias_file, ln_altbias, &
275	& endailyavtypes, ln_s_at_t, ln_profb_ena, &
276	& ln_vel3d, ln_velavcur, velavcurfiles, &
277	& ln_velhrcur, velhrcurfiles, &
278	& ln_velavadcp, velavadcpfiles, &
279	& ln_velhradcp, velhradcpfiles, &
280	& ln_velfb, velfbfiles, ln_velfb_av, &
281	& ln_logchl, ln_logchlfb, &
282	& logchlfiles, logchlfbfiles, &
283
284
285
286	! JOZEF ADDITION
287
288	& ln_logchlpft, ln_logchlpftfb, &
289	& logchlpftfiles, logchlpftfbfiles, &
290
291	! END
292
293
294
295
296	& ln_spm, ln_spmfb, &
297	& spmfiles, spmfbfiles, &
298	& ln_fco2, ln_fco2fb, &
299	& fco2files, fco2fbfiles, &
300	& ln_pco2, ln_pco2fb, &
301	& pco2files, pco2fbfiles, &
302	& ln_profb_enatim, ln_ignmis, ln_cl4, &
303	& ln_sstbias, sstbias_files
304
305	INTEGER :: jprofset
306	INTEGER :: jveloset
307	INTEGER :: jvar
308	INTEGER :: jnumenact
309	INTEGER :: jnumcorio
310	INTEGER :: jnumprofb
311	INTEGER :: jnumslaact
312	INTEGER :: jnumslapas
313	INTEGER :: jnumslafb
314	INTEGER :: jnumsst
315	INTEGER :: jnumsstfb
316	INTEGER :: jnumsstbias
317	INTEGER :: jnumseaice
318	INTEGER :: jnumvelavcur
319	INTEGER :: jnumvelhrcur
320	INTEGER :: jnumvelavadcp
321	INTEGER :: jnumvelhradcp
322	INTEGER :: jnumvelfb
323	INTEGER :: jnumlogchl
324	INTEGER :: jnumlogchlfb
325
326
327
328
329	! JOZEF ADDITION
330	! INTEGER :: nn_logchlpftscc
331
332	! INTEGER, DIMENSION(:), ALLOCATABLE :: jnumlogchlpft
333	! INTEGER, DIMENSION(:), ALLOCATABLE :: jnumlogchlpftfb
334
335	INTEGER :: jnumlogchlpft(nn_logchlpftscc)
336	INTEGER :: jnumlogchlpftfb(nn_logchlpftscc)
337
338
339	! END
340
341
342
343
344
345	INTEGER :: jnumspm
346	INTEGER :: jnumspmfb
347	INTEGER :: jnumfco2
348	INTEGER :: jnumfco2fb
349	INTEGER :: jnumpco2
350	INTEGER :: jnumpco2fb
351	INTEGER :: ji
352	INTEGER :: jset
353	INTEGER :: ios ! Local integer output status for namelist read
354	LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
355
356
357	! JOZEF ADDITION
358
359	INTEGER :: pft
360
361	!END
362
363
364
365	!-----------------------------------------------------------------------
366	! Read namelist parameters
367	!-----------------------------------------------------------------------
368
369	ln_logchl = .FALSE.
370	ln_logchlfb = .FALSE.
371
372
373
374
375
376
377	! JOZEF ADDITION
378	ln_logchlpft = .FALSE.
379	ln_logchlpftfb = .FALSE.
380
381	! #if defined key_fabm
382	! nn_logchlpftscc = 4
383	! #elif defined key_medusa
384	! nn_logchlpftscc = 2
385	! #elif defined key_hadocc
386	! nn_logchlpftscc = 1
387	! #else
388	! nn_logchlpftscc = 0
389	! #fi
390
391	! END
392
393
394
395
396
397	ln_spm = .FALSE.
398	ln_spmfb = .FALSE.
399	ln_fco2 = .FALSE.
400	ln_fco2fb = .FALSE.
401	ln_pco2 = .FALSE.
402	ln_pco2fb = .FALSE.
403
404	!Initalise all values in namelist arrays
405	enactfiles(:) = ''
406	coriofiles(:) = ''
407	profbfiles(:) = ''
408	slafilesact(:) = ''
409	slafilespas(:) = ''
410	slafbfiles(:) = ''
411	sstfiles(:) = ''
412	sstfbfiles(:) = ''
413	seaicefiles(:) = ''
414	velcurfiles(:) = ''
415	veladcpfiles(:) = ''
416	velavcurfiles(:) = ''
417	velhrcurfiles(:) = ''
418	velavadcpfiles(:) = ''
419	velhradcpfiles(:) = ''
420	velfbfiles(:) = ''
421	velcurfiles(:) = ''
422	veladcpfiles(:) = ''
423	logchlfiles(:) = ''
424	logchlfbfiles(:) = ''
425
426
427
428
429	! JOZEF ADDITION
430
431	ALLOCATE(logchlpftfiles(MaxNumFiles, nn_logchlpftscc))
432	ALLOCATE(logchlpftfbfiles(MaxNumFiles, nn_logchlpftscc))
433
434	DO pft = 1,nn_logchlpftscc
435	logchlpftfiles(:,pft) = ''
436	logchlpftfbfiles(:,pft) = ''
437	END DO
438	! END
439
440
441
442
443
444
445
446	spmfiles(:) = ''
447	spmfbfiles(:) = ''
448	fco2files(:) = ''
449	fco2fbfiles(:) = ''
450	pco2files(:) = ''
451	pco2fbfiles(:) = ''
452	sstbias_files(:) = ''
453	endailyavtypes(:) = -1
454	endailyavtypes(1) = 820
455	ln_profb_ena(:) = .FALSE.
456	ln_profb_enatim(:) = .TRUE.
457	ln_velfb_av(:) = .FALSE.
458	ln_ignmis = .FALSE.
459	ln_bound_reject = .TRUE.
460
461	! Read Namelist namobs : control observation diagnostics
462	REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation
463	READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
464	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
465
466	CALL ini_date( dobsini )
467	CALL fin_date( dobsend )
468
469	REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation
470	READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
471	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
472	IF(lwm) WRITE ( numond, namobs )
473
474	! Count number of files for each type
475	IF (ln_ena) THEN
476	lmask(:) = .FALSE.
477	WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
478	jnumenact = COUNT(lmask)
479	ENDIF
480	IF (ln_cor) THEN
481	lmask(:) = .FALSE.
482	WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
483	jnumcorio = COUNT(lmask)
484	ENDIF
485	IF (ln_profb) THEN
486	lmask(:) = .FALSE.
487	WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
488	jnumprofb = COUNT(lmask)
489	ENDIF
490	IF (ln_sladt) THEN
491	lmask(:) = .FALSE.
492	WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
493	jnumslaact = COUNT(lmask)
494	lmask(:) = .FALSE.
495	WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
496	jnumslapas = COUNT(lmask)
497	ENDIF
498	IF (ln_slafb) THEN
499	lmask(:) = .FALSE.
500	WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
501	jnumslafb = COUNT(lmask)
502	lmask(:) = .FALSE.
503	ENDIF
504	IF (ln_ghrsst) THEN
505	lmask(:) = .FALSE.
506	WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
507	jnumsst = COUNT(lmask)
508	ENDIF
509	IF (ln_sstfb) THEN
510	lmask(:) = .FALSE.
511	WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
512	jnumsstfb = COUNT(lmask)
513	lmask(:) = .FALSE.
514	ENDIF
515	IF (ln_sstbias) THEN
516	lmask(:) = .FALSE.
517	WHERE (sstbias_files(:) /= '') lmask(:) = .TRUE.
518	jnumsstbias = COUNT(lmask)
519	lmask(:) = .FALSE.
520	ENDIF
521	IF (ln_seaice) THEN
522	lmask(:) = .FALSE.
523	WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
524	jnumseaice = COUNT(lmask)
525	ENDIF
526	IF (ln_velavcur) THEN
527	lmask(:) = .FALSE.
528	WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
529	jnumvelavcur = COUNT(lmask)
530	ENDIF
531	IF (ln_velhrcur) THEN
532	lmask(:) = .FALSE.
533	WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
534	jnumvelhrcur = COUNT(lmask)
535	ENDIF
536	IF (ln_velavadcp) THEN
537	lmask(:) = .FALSE.
538	WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
539	jnumvelavadcp = COUNT(lmask)
540	ENDIF
541	IF (ln_velhradcp) THEN
542	lmask(:) = .FALSE.
543	WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
544	jnumvelhradcp = COUNT(lmask)
545	ENDIF
546	IF (ln_velfb) THEN
547	lmask(:) = .FALSE.
548	WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
549	jnumvelfb = COUNT(lmask)
550	lmask(:) = .FALSE.
551	ENDIF
552	IF (ln_logchl) THEN
553	lmask(:) = .FALSE.
554	WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE.
555	jnumlogchl = COUNT(lmask)
556	ENDIF
557	IF (ln_logchlfb) THEN
558	lmask(:) = .FALSE.
559	WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE.
560	jnumlogchlfb = COUNT(lmask)
561	ENDIF
562
563
564
565
566
567
568
569
570
571
572	! JOZEF ADDITION
573
574	DO pft=1,nn_logchlpftscc
575
576	IF (ln_logchlpft) THEN
577
578	lmask(:) = .FALSE.
579	WHERE (logchlpftfiles(:,pft) /= '') lmask(:) = .TRUE.
580	jnumlogchlpft(pft) = COUNT(lmask)
581
582	ENDIF
583	IF (ln_logchlpftfb) THEN
584	lmask(:) = .FALSE.
585	WHERE (logchlpftfbfiles(:,pft) /= '') lmask(:) = .TRUE.
586	jnumlogchlpftfb(pft) = COUNT(lmask)
587
588
589	ENDIF
590
591	!WRITE(numout,*)'jnumlogchlpftfb is', jnumlogchlpftfbfiles(:,pft)
592
593	END DO
594	! END
595
596
597
598	WRITE(numout,*)'jnumlogchlpftfb is', logchlpftfbfiles
599
600
601
602
603
604	IF (ln_spm) THEN
605	lmask(:) = .FALSE.
606	WHERE (spmfiles(:) /= '') lmask(:) = .TRUE.
607	jnumspm = COUNT(lmask)
608	ENDIF
609	IF (ln_spmfb) THEN
610	lmask(:) = .FALSE.
611	WHERE (spmfbfiles(:) /= '') lmask(:) = .TRUE.
612	jnumspmfb = COUNT(lmask)
613	ENDIF
614	IF (ln_fco2) THEN
615	lmask(:) = .FALSE.
616	WHERE (fco2files(:) /= '') lmask(:) = .TRUE.
617	jnumfco2 = COUNT(lmask)
618	ENDIF
619	IF (ln_fco2fb) THEN
620	lmask(:) = .FALSE.
621	WHERE (fco2fbfiles(:) /= '') lmask(:) = .TRUE.
622	jnumfco2fb = COUNT(lmask)
623	ENDIF
624	IF (ln_pco2) THEN
625	lmask(:) = .FALSE.
626	WHERE (pco2files(:) /= '') lmask(:) = .TRUE.
627	jnumpco2 = COUNT(lmask)
628	ENDIF
629	IF (ln_pco2fb) THEN
630	lmask(:) = .FALSE.
631	WHERE (pco2fbfiles(:) /= '') lmask(:) = .TRUE.
632	jnumpco2fb = COUNT(lmask)
633	ENDIF
634
635	! Control print
636	IF(lwp) THEN
637	WRITE(numout,*)
638	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
639	WRITE(numout,*) '~~~~~~~~~~~~'
640	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
641	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
642	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
643	WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena
644	WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor
645	WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb
646	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
647	WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt
648	WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb
649	WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh
650	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
651	WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst
652	WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst
653	WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb
654	WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias
655	WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight
656	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
657	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice
658	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
659	WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur
660	WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur
661	WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp
662	WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp
663	WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb
664	WRITE(numout,*) ' Logical switch for logchl observations ln_logchl = ', ln_logchl
665	WRITE(numout,*) ' Logical switch for feedback logchl data ln_logchlfb = ', ln_logchlfb
666
667
668
669
670
671	! JOZEF ADDITION
672
673	WRITE(numout,*) ' Logical switch for logchlpft observations ln_logchlpft = ', ln_logchlpft
674	WRITE(numout,*) ' Logical switch for feedback logchlpft data ln_logchlpftfb = ', ln_logchlpftfb
675
676	! END
677
678
679
680
681
682	WRITE(numout,*) ' Logical switch for spm observations ln_spm = ', ln_spm
683	WRITE(numout,*) ' Logical switch for feedback spm data ln_spmfb = ', ln_spmfb
684	WRITE(numout,*) ' Logical switch for fco2 observations ln_fco2 = ', ln_fco2
685	WRITE(numout,*) ' Logical switch for pco2 observations ln_pco2 = ', ln_pco2
686	WRITE(numout,*) ' Logical switch for feedback pco2 data ln_pco2fb = ', ln_pco2fb
687	WRITE(numout,*) ' Logical switch for feedback fco2 data ln_fco2fb = ', ln_fco2fb
688	WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global
689	WRITE(numout,*) &
690	' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup
691	IF (ln_grid_search_lookup) &
692	WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file
693	IF (ln_ena) THEN
694	DO ji = 1, jnumenact
695	WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', &
696	TRIM(enactfiles(ji))
697	END DO
698	ENDIF
699	IF (ln_cor) THEN
700	DO ji = 1, jnumcorio
701	WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', &
702	TRIM(coriofiles(ji))
703	END DO
704	ENDIF
705	IF (ln_profb) THEN
706	DO ji = 1, jnumprofb
707	IF (ln_profb_ena(ji)) THEN
708	WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', &
709	TRIM(profbfiles(ji))
710	ELSE
711	WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', &
712	TRIM(profbfiles(ji))
713	ENDIF
714	WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji)
715	END DO
716	ENDIF
717	IF (ln_sladt) THEN
718	DO ji = 1, jnumslaact
719	WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', &
720	TRIM(slafilesact(ji))
721	END DO
722	DO ji = 1, jnumslapas
723	WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', &
724	TRIM(slafilespas(ji))
725	END DO
726	ENDIF
727	IF (ln_slafb) THEN
728	DO ji = 1, jnumslafb
729	WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', &
730	TRIM(slafbfiles(ji))
731	END DO
732	ENDIF
733	IF (ln_ghrsst) THEN
734	DO ji = 1, jnumsst
735	WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', &
736	TRIM(sstfiles(ji))
737	END DO
738	ENDIF
739	IF (ln_sstfb) THEN
740	DO ji = 1, jnumsstfb
741	WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', &
742	TRIM(sstfbfiles(ji))
743	END DO
744	ENDIF
745	IF (ln_seaice) THEN
746	DO ji = 1, jnumseaice
747	WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', &
748	TRIM(seaicefiles(ji))
749	END DO
750	ENDIF
751	IF (ln_velavcur) THEN
752	DO ji = 1, jnumvelavcur
753	WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', &
754	TRIM(velavcurfiles(ji))
755	END DO
756	ENDIF
757	IF (ln_velhrcur) THEN
758	DO ji = 1, jnumvelhrcur
759	WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', &
760	TRIM(velhrcurfiles(ji))
761	END DO
762	ENDIF
763	IF (ln_velavadcp) THEN
764	DO ji = 1, jnumvelavadcp
765	WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', &
766	TRIM(velavadcpfiles(ji))
767	END DO
768	ENDIF
769	IF (ln_velhradcp) THEN
770	DO ji = 1, jnumvelhradcp
771	WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', &
772	TRIM(velhradcpfiles(ji))
773	END DO
774	ENDIF
775	IF (ln_velfb) THEN
776	DO ji = 1, jnumvelfb
777	IF (ln_velfb_av(ji)) THEN
778	WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', &
779	TRIM(velfbfiles(ji))
780	ELSE
781	WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', &
782	TRIM(velfbfiles(ji))
783	ENDIF
784	END DO
785	ENDIF
786
787
788	IF (ln_logchl) THEN
789	DO ji = 1, jnumlogchl
790	WRITE(numout,'(1X,2A)') ' logchl input observation file name logchlfiles = ', &
791	TRIM(logchlfiles(ji))
792	END DO
793	ENDIF
794	IF (ln_logchlfb) THEN
795	DO ji = 1, jnumlogchlfb
796	WRITE(numout,'(1X,2A)') ' Feedback logchl input observation file name logchlfbfiles = ', &
797	TRIM(logchlfbfiles(ji))
798	END DO
799	ENDIF
800
801
802
803
804
805	! JOZEF ADDITION
806
807	IF (ln_logchlpft) THEN
808
809	DO pft=1,nn_logchlpftscc
810
811	DO ji = 1, jnumlogchlpft(pft)
812	WRITE(numout,'(1X,2A)') ' logchlpft input observation file name logchlfiles = ', &
813	TRIM(logchlpftfiles(ji,pft))
814	END DO
815	END DO
816	ENDIF
817	IF (ln_logchlpftfb) THEN
818
819	DO pft=1,nn_logchlpftscc
820
821	DO ji = 1, jnumlogchlpftfb(pft)
822	WRITE(numout,'(1X,2A)') ' Feedback logchlpft input observation file name logchlfbfiles = ', &
823	TRIM(logchlpftfbfiles(ji,pft))
824	END DO
825	END DO
826
827
828	ENDIF
829
830
831	!END
832
833
834
835
836
837
838
839
840	IF (ln_spm) THEN
841	DO ji = 1, jnumspm
842	WRITE(numout,'(1X,2A)') ' spm input observation file name spmfiles = ', &
843	TRIM(spmfiles(ji))
844	END DO
845	ENDIF
846	IF (ln_spmfb) THEN
847	DO ji = 1, jnumspmfb
848	WRITE(numout,'(1X,2A)') ' Feedback spm input observation file name spmfbfiles = ', &
849	TRIM(spmfbfiles(ji))
850	END DO
851	ENDIF
852	IF (ln_fco2) THEN
853	DO ji = 1, jnumfco2
854	WRITE(numout,'(1X,2A)') ' fco2 input observation file name fco2files = ', &
855	TRIM(fco2files(ji))
856	END DO
857	ENDIF
858	IF (ln_fco2fb) THEN
859	DO ji = 1, jnumfco2fb
860	WRITE(numout,'(1X,2A)') ' Feedback fco2 input observation file name fco2fbfiles = ', &
861	TRIM(fco2fbfiles(ji))
862	END DO
863	ENDIF
864	IF (ln_pco2) THEN
865	DO ji = 1, jnumpco2
866	WRITE(numout,'(1X,2A)') ' pco2 input observation file name pco2files = ', &
867	TRIM(pco2files(ji))
868	END DO
869	ENDIF
870	IF (ln_pco2fb) THEN
871	DO ji = 1, jnumpco2fb
872	WRITE(numout,'(1X,2A)') ' Feedback pco2 input observation file name pco2fbfiles = ', &
873	TRIM(pco2fbfiles(ji))
874	END DO
875	ENDIF
876	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini
877	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend
878	WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint
879	WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint
880	WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea
881	WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject
882	WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc
883	WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr
884	WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff
885	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
886	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
887	WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes
888
889	ENDIF
890
891	IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
892	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
893	RETURN
894	ENDIF
895
896	IF ( ln_grid_global ) THEN
897	CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
898	ENDIF
899
900	CALL obs_typ_init
901
902	IF ( ln_grid_global ) THEN
903	CALL mppmap_init
904	ENDIF
905
906	! Parameter control
907	#if defined key_diaobs
908	IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
909	& ( .NOT. ln_vel3d ).AND. &
910	& ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
911	& ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
912	& ( .NOT. ln_spm ).AND.( .NOT. ln_fco2 ).AND.( .NOT. ln_pco2 ) ) THEN
913	IF(lwp) WRITE(numout,cform_war)
914	IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
915	& ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
916	& ' ln_logchl, ln_spm, ln_fco2, ln_pco2 are all set to .FALSE.'
917	nwarn = nwarn + 1
918	ENDIF
919	#endif
920
921	CALL obs_grid_setup( )
922	IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
923	CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
924	& ' is not available')
925	ENDIF
926	IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
927	CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
928	& ' is not available')
929	ENDIF
930
931	!-----------------------------------------------------------------------
932	! Depending on switches read the various observation types
933	!-----------------------------------------------------------------------
934	! - Temperature/salinity profiles
935
936	IF ( ln_t3d .OR. ln_s3d ) THEN
937
938	! Set the number of variables for profiles to 2 (T and S)
939	nprofvars = 2
940	! Set the number of extra variables for profiles to 1 (insitu temp).
941	nprofextr = 1
942
943	! Count how may insitu data sets we have and allocate data.
944	jprofset = 0
945	IF ( ln_ena ) jprofset = jprofset + 1
946	IF ( ln_cor ) jprofset = jprofset + 1
947	IF ( ln_profb ) jprofset = jprofset + jnumprofb
948	nprofsets = jprofset
949	IF ( nprofsets > 0 ) THEN
950	ALLOCATE(ld_enact(nprofsets))
951	ALLOCATE(profdata(nprofsets))
952	ALLOCATE(prodatqc(nprofsets))
953	ENDIF
954
955	jprofset = 0
956
957	! ENACT insitu data
958
959	IF ( ln_ena ) THEN
960
961	jprofset = jprofset + 1
962
963	ld_enact(jprofset) = .TRUE.
964
965	CALL obs_rea_pro_dri( 1, profdata(jprofset), &
966	& jnumenact, enactfiles(1:jnumenact), &
967	& nprofvars, nprofextr, &
968	& nitend-nit000+2, &
969	& dobsini, dobsend, ln_t3d, ln_s3d, &
970	& ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
971	& kdailyavtypes = endailyavtypes )
972
973	DO jvar = 1, 2
974
975	CALL obs_prof_staend( profdata(jprofset), jvar )
976
977	END DO
978
979	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
980	& ln_t3d, ln_s3d, ln_nea, &
981	& kdailyavtypes=endailyavtypes )
982
983	ENDIF
984
985	! Coriolis insitu data
986
987	IF ( ln_cor ) THEN
988
989	jprofset = jprofset + 1
990
991	ld_enact(jprofset) = .FALSE.
992
993	CALL obs_rea_pro_dri( 2, profdata(jprofset), &
994	& jnumcorio, coriofiles(1:jnumcorio), &
995	& nprofvars, nprofextr, &
996	& nitend-nit000+2, &
997	& dobsini, dobsend, ln_t3d, ln_s3d, &
998	& ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
999
1000	DO jvar = 1, 2
1001
1002	CALL obs_prof_staend( profdata(jprofset), jvar )
1003
1004	END DO
1005
1006	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
1007	& ln_t3d, ln_s3d, ln_nea )
1008
1009	ENDIF
1010
1011	! Feedback insitu data
1012
1013	IF ( ln_profb ) THEN
1014
1015	DO jset = 1, jnumprofb
1016
1017	jprofset = jprofset + 1
1018	ld_enact (jprofset) = ln_profb_ena(jset)
1019
1020	CALL obs_rea_pro_dri( 0, profdata(jprofset), &
1021	& 1, profbfiles(jset:jset), &
1022	& nprofvars, nprofextr, &
1023	& nitend-nit000+2, &
1024	& dobsini, dobsend, ln_t3d, ln_s3d, &
1025	& ln_ignmis, ln_s_at_t, &
1026	& ld_enact(jprofset).AND.&
1027	& ln_profb_enatim(jset), &
1028	& .FALSE., kdailyavtypes = endailyavtypes )
1029
1030	DO jvar = 1, 2
1031
1032	CALL obs_prof_staend( profdata(jprofset), jvar )
1033
1034	END DO
1035
1036	IF ( ld_enact(jprofset) ) THEN
1037	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
1038	& ln_t3d, ln_s3d, ln_nea, &
1039	& kdailyavtypes = endailyavtypes )
1040	ELSE
1041	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
1042	& ln_t3d, ln_s3d, ln_nea )
1043	ENDIF
1044
1045	END DO
1046
1047	ENDIF
1048
1049	ENDIF
1050
1051	! - Sea level anomalies
1052	IF ( ln_sla ) THEN
1053	! Set the number of variables for sla to 1
1054	nslavars = 1
1055
1056	! Set the number of extra variables for sla to 2
1057	nslaextr = 2
1058
1059	! Set the number of sla data sets to 2
1060	nslasets = 0
1061	IF ( ln_sladt ) THEN
1062	nslasets = nslasets + 2
1063	ENDIF
1064	IF ( ln_slafb ) THEN
1065	nslasets = nslasets + jnumslafb
1066	ENDIF
1067
1068	ALLOCATE(sladata(nslasets))
1069	ALLOCATE(sladatqc(nslasets))
1070	sladata(:)%nsurf=0
1071	sladatqc(:)%nsurf=0
1072
1073	nslasets = 0
1074
1075	! AVISO SLA data
1076
1077	IF ( ln_sladt ) THEN
1078
1079	! Active SLA observations
1080
1081	nslasets = nslasets + 1
1082
1083	CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
1084	& slafilesact(1:jnumslaact), &
1085	& nslavars, nslaextr, nitend-nit000+2, &
1086	& dobsini, dobsend, ln_ignmis, .FALSE. )
1087	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
1088	& ln_sla, ln_nea )
1089
1090	! Passive SLA observations
1091
1092	nslasets = nslasets + 1
1093
1094	CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
1095	& slafilespas(1:jnumslapas), &
1096	& nslavars, nslaextr, nitend-nit000+2, &
1097	& dobsini, dobsend, ln_ignmis, .FALSE. )
1098
1099	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
1100	& ln_sla, ln_nea )
1101
1102	ENDIF
1103
1104	! Feedback SLA data
1105
1106	IF ( ln_slafb ) THEN
1107
1108	DO jset = 1, jnumslafb
1109
1110	nslasets = nslasets + 1
1111
1112	CALL obs_rea_sla( 0, sladata(nslasets), 1, &
1113	& slafbfiles(jset:jset), &
1114	& nslavars, nslaextr, nitend-nit000+2, &
1115	& dobsini, dobsend, ln_ignmis, .FALSE. )
1116	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
1117	& ln_sla, ln_nea )
1118
1119	END DO
1120
1121	ENDIF
1122
1123	CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
1124
1125	! read in altimeter bias
1126
1127	IF ( ln_altbias ) THEN
1128	CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
1129	ENDIF
1130
1131	ENDIF
1132
1133	! - Sea surface height
1134	IF ( ln_ssh ) THEN
1135	IF(lwp) WRITE(numout,*) ' SSH currently not available'
1136	ENDIF
1137
1138	! - Sea surface temperature
1139	IF ( ln_sst ) THEN
1140
1141	! Set the number of variables for sst to 1
1142	nsstvars = 1
1143
1144	! Set the number of extra variables for sst to 0
1145	nsstextr = 0
1146
1147	nsstsets = 0
1148
1149	IF (ln_reysst) nsstsets = nsstsets + 1
1150	IF (ln_ghrsst) nsstsets = nsstsets + 1
1151	IF ( ln_sstfb ) THEN
1152	nsstsets = nsstsets + jnumsstfb
1153	ENDIF
1154
1155	ALLOCATE(sstdata(nsstsets))
1156	ALLOCATE(sstdatqc(nsstsets))
1157	ALLOCATE(ld_sstnight(nsstsets))
1158	sstdata(:)%nsurf=0
1159	sstdatqc(:)%nsurf=0
1160	ld_sstnight(:)=.false.
1161
1162	nsstsets = 0
1163
1164	IF (ln_reysst) THEN
1165
1166	nsstsets = nsstsets + 1
1167
1168	ld_sstnight(nsstsets) = ln_sstnight
1169
1170	CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
1171	& nsstvars, nsstextr, &
1172	& nitend-nit000+2, dobsini, dobsend )
1173	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
1174	& ln_nea )
1175
1176	ENDIF
1177
1178	IF (ln_ghrsst) THEN
1179
1180	nsstsets = nsstsets + 1
1181
1182	ld_sstnight(nsstsets) = ln_sstnight
1183
1184	CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
1185	& sstfiles(1:jnumsst), &
1186	& nsstvars, nsstextr, nitend-nit000+2, &
1187	& dobsini, dobsend, ln_ignmis, .FALSE. )
1188	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
1189	& ln_nea )
1190
1191	ENDIF
1192
1193	! Feedback SST data
1194
1195	IF ( ln_sstfb ) THEN
1196
1197	DO jset = 1, jnumsstfb
1198
1199	nsstsets = nsstsets + 1
1200
1201	ld_sstnight(nsstsets) = ln_sstnight
1202
1203	CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
1204	& sstfbfiles(jset:jset), &
1205	& nsstvars, nsstextr, nitend-nit000+2, &
1206	& dobsini, dobsend, ln_ignmis, .FALSE. )
1207	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
1208	& ln_sst, ln_nea )
1209
1210	END DO
1211
1212	ENDIF
1213
1214	!Read in bias field and correct SST.
1215	IF ( ln_sstbias ) THEN
1216	IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// &
1217	" but no bias"// &
1218	" files to read in")
1219	CALL obs_app_sstbias( nsstsets, sstdatqc, n2dint, &
1220	jnumsstbias, &
1221	sstbias_files(1:jnumsstbias) )
1222	ENDIF
1223
1224	ENDIF
1225
1226	! - Sea surface salinity
1227	IF ( ln_sss ) THEN
1228	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1229	ENDIF
1230
1231	! - Sea Ice Concentration
1232
1233	IF ( ln_seaice ) THEN
1234
1235	! Set the number of variables for seaice to 1
1236	nseaicevars = 1
1237
1238	! Set the number of extra variables for seaice to 0
1239	nseaiceextr = 0
1240
1241	! Set the number of data sets to 1
1242	nseaicesets = 1
1243
1244	ALLOCATE(seaicedata(nseaicesets))
1245	ALLOCATE(seaicedatqc(nseaicesets))
1246	seaicedata(:)%nsurf=0
1247	seaicedatqc(:)%nsurf=0
1248
1249	CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
1250	& seaicefiles(1:jnumseaice), &
1251	& nseaicevars, nseaiceextr, nitend-nit000+2, &
1252	& dobsini, dobsend, ln_ignmis, .FALSE. )
1253
1254	CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
1255	& ln_seaice, ln_nea )
1256
1257	ENDIF
1258
1259	IF (ln_vel3d) THEN
1260
1261	! Set the number of variables for profiles to 2 (U and V)
1262	nvelovars = 2
1263
1264	! Set the number of extra variables for profiles to 2 to store
1265	! rotation parameters
1266	nveloextr = 2
1267
1268	jveloset = 0
1269
1270	IF ( ln_velavcur ) jveloset = jveloset + 1
1271	IF ( ln_velhrcur ) jveloset = jveloset + 1
1272	IF ( ln_velavadcp ) jveloset = jveloset + 1
1273	IF ( ln_velhradcp ) jveloset = jveloset + 1
1274	IF (ln_velfb) jveloset = jveloset + jnumvelfb
1275
1276	nvelosets = jveloset
1277	IF ( nvelosets > 0 ) THEN
1278	ALLOCATE( velodata(nvelosets) )
1279	ALLOCATE( veldatqc(nvelosets) )
1280	ALLOCATE( ld_velav(nvelosets) )
1281	ENDIF
1282
1283	jveloset = 0
1284
1285	! Daily averaged data
1286
1287	IF ( ln_velavcur ) THEN
1288
1289	jveloset = jveloset + 1
1290
1291	ld_velav(jveloset) = .TRUE.
1292
1293	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
1294	& velavcurfiles(1:jnumvelavcur), &
1295	& nvelovars, nveloextr, &
1296	& nitend-nit000+2, &
1297	& dobsini, dobsend, ln_ignmis, &
1298	& ld_velav(jveloset), &
1299	& .FALSE. )
1300
1301	DO jvar = 1, 2
1302	CALL obs_prof_staend( velodata(jveloset), jvar )
1303	END DO
1304
1305	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1306	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1307
1308	ENDIF
1309
1310	! High frequency data
1311
1312	IF ( ln_velhrcur ) THEN
1313
1314	jveloset = jveloset + 1
1315
1316	ld_velav(jveloset) = .FALSE.
1317
1318	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
1319	& velhrcurfiles(1:jnumvelhrcur), &
1320	& nvelovars, nveloextr, &
1321	& nitend-nit000+2, &
1322	& dobsini, dobsend, ln_ignmis, &
1323	& ld_velav(jveloset), &
1324	& .FALSE. )
1325
1326	DO jvar = 1, 2
1327	CALL obs_prof_staend( velodata(jveloset), jvar )
1328	END DO
1329
1330	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1331	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1332
1333	ENDIF
1334
1335	! Daily averaged data
1336
1337	IF ( ln_velavadcp ) THEN
1338
1339	jveloset = jveloset + 1
1340
1341	ld_velav(jveloset) = .TRUE.
1342
1343	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
1344	& velavadcpfiles(1:jnumvelavadcp), &
1345	& nvelovars, nveloextr, &
1346	& nitend-nit000+2, &
1347	& dobsini, dobsend, ln_ignmis, &
1348	& ld_velav(jveloset), &
1349	& .FALSE. )
1350
1351	DO jvar = 1, 2
1352	CALL obs_prof_staend( velodata(jveloset), jvar )
1353	END DO
1354
1355	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1356	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1357
1358	ENDIF
1359
1360	! High frequency data
1361
1362	IF ( ln_velhradcp ) THEN
1363
1364	jveloset = jveloset + 1
1365
1366	ld_velav(jveloset) = .FALSE.
1367
1368	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
1369	& velhradcpfiles(1:jnumvelhradcp), &
1370	& nvelovars, nveloextr, &
1371	& nitend-nit000+2, &
1372	& dobsini, dobsend, ln_ignmis, &
1373	& ld_velav(jveloset), &
1374	& .FALSE. )
1375
1376	DO jvar = 1, 2
1377	CALL obs_prof_staend( velodata(jveloset), jvar )
1378	END DO
1379
1380	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1381	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1382
1383	ENDIF
1384
1385	IF ( ln_velfb ) THEN
1386
1387	DO jset = 1, jnumvelfb
1388
1389	jveloset = jveloset + 1
1390
1391	ld_velav(jveloset) = ln_velfb_av(jset)
1392
1393	CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
1394	& velfbfiles(jset:jset), &
1395	& nvelovars, nveloextr, &
1396	& nitend-nit000+2, &
1397	& dobsini, dobsend, ln_ignmis, &
1398	& ld_velav(jveloset), &
1399	& .FALSE. )
1400
1401	DO jvar = 1, 2
1402	CALL obs_prof_staend( velodata(jveloset), jvar )
1403	END DO
1404
1405	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1406	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1407
1408
1409	END DO
1410
1411	ENDIF
1412
1413	ENDIF
1414
1415	! - log10(chlorophyll)
1416
1417	IF ( ln_logchl ) THEN
1418
1419	! Set the number of variables for logchl to 1
1420	nlogchlvars = 1
1421
1422	! Set the number of extra variables for logchl to 0
1423	nlogchlextr = 0
1424
1425	IF ( ln_logchlfb ) THEN
1426	nlogchlsets = jnumlogchlfb
1427	ELSE
1428	nlogchlsets = 1
1429	ENDIF
1430
1431	ALLOCATE(logchldata(nlogchlsets))
1432	ALLOCATE(logchldatqc(nlogchlsets))
1433	logchldata(:)%nsurf=0
1434	logchldatqc(:)%nsurf=0
1435
1436	nlogchlsets = 0
1437
1438
1439
1440
1441
1442
1443	IF ( ln_logchlfb ) THEN ! Feedback file format
1444
1445	DO jset = 1, jnumlogchlfb
1446
1447	nlogchlsets = nlogchlsets + 1
1448
1449	CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, &
1450	& logchlfbfiles(jset:jset), &
1451	& nlogchlvars, nlogchlextr, nitend-nit000+2, &
1452	& dobsini, dobsend, ln_ignmis, .FALSE. )
1453
1454	CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1455	& ln_logchl, ln_nea )
1456
1457	ENDDO
1458
1459	ELSE ! Original file format
1460
1461	nlogchlsets = nlogchlsets + 1
1462
1463	CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, &
1464	& logchlfiles(1:jnumlogchl), &
1465	& nlogchlvars, nlogchlextr, nitend-nit000+2, &
1466	& dobsini, dobsend, ln_ignmis, .FALSE. )
1467
1468	CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1469	& ln_logchl, ln_nea )
1470
1471	ENDIF
1472
1473	ENDIF
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484	! JOZEF ADDITION
1485
1486
1487	IF ( ln_logchlpft ) THEN
1488
1489	! Set the number of variables for logchl to 1
1490	nlogchlpftvars = 1 ! ??????? ASK DAVID
1491
1492	! Set the number of extra variables for logchl to 0
1493	nlogchlpftextr = 0
1494
1495	DO pft = 1, nn_logchlpftscc
1496
1497
1498
1499	IF ( ln_logchlpftfb ) THEN
1500	nlogchlpftsets = jnumlogchlpftfb(pft)
1501	ELSE
1502	nlogchlpftsets = 1
1503	ENDIF
1504
1505
1506
1507	END DO
1508
1509	ALLOCATE(logchlpftdata(nn_logchlpftscc))
1510	ALLOCATE(logchlpftdatqc(nn_logchlpftscc))
1511
1512
1513
1514	DO pft = 1, nn_logchlpftscc
1515
1516	logchlpftdata(pft)%nsurf=0
1517	logchlpftdatqc(pft)%nsurf=0
1518
1519	nlogchlpftsets = 0
1520
1521	WRITE(numout,*)'jnumlogchlpftfb is', jnumlogchlpftfb, ln_logchlpftfb
1522
1523	IF ( ln_logchlpftfb ) THEN ! Feedback file format
1524
1525	DO jset = 1, jnumlogchlpftfb(pft)
1526
1527	nlogchlpftsets = nlogchlpftsets + 1
1528
1529	CALL obs_rea_logchl( 0, logchlpftdata(pft), 1, &
1530	& logchlpftfbfiles(jset:jset,pft), &
1531	& nlogchlpftvars, nlogchlpftextr, nitend-nit000+2, &
1532	& dobsini, dobsend, ln_ignmis, .FALSE. )
1533
1534	CALL obs_pre_logchl( logchlpftdata(pft), logchlpftdatqc(pft), &
1535	& ln_logchlpft, ln_nea )
1536
1537	ENDDO
1538
1539	ELSE ! Original file format
1540
1541	nlogchlpftsets = 1
1542
1543	CALL obs_rea_logchl( 1, logchlpftdata(pft), jnumlogchlpft(pft), &
1544	& logchlpftfiles(1:jnumlogchlpft(pft),pft), &
1545	& nlogchlpftvars, nlogchlpftextr, nitend-nit000+2, &
1546	& dobsini, dobsend, ln_ignmis, .FALSE. )
1547
1548	CALL obs_pre_logchl( logchlpftdata(pft), logchlpftdatqc(pft), &
1549	& ln_logchlpft, ln_nea )
1550
1551	ENDIF
1552
1553	ENDDO
1554
1555
1556
1557
1558
1559	ENDIF
1560
1561
1562	nlogchlpftsets = 1 ! THIS SPECIAL CASE NEEDS TO BE GENERALIZED LATER
1563
1564	WRITE(numout,*)'Producing nlogchlpftsets', nlogchlpftsets
1565
1566	! END
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586	! - spm
1587
1588	IF ( ln_spm ) THEN
1589
1590	! Set the number of variables for spm to 1
1591	nspmvars = 1
1592
1593	! Set the number of extra variables for spm to 0
1594	nspmextr = 0
1595
1596	IF ( ln_spmfb ) THEN
1597	nspmsets = jnumspmfb
1598	ELSE
1599	nspmsets = 1
1600	ENDIF
1601
1602	ALLOCATE(spmdata(nspmsets))
1603	ALLOCATE(spmdatqc(nspmsets))
1604	spmdata(:)%nsurf=0
1605	spmdatqc(:)%nsurf=0
1606
1607	nspmsets = 0
1608
1609	IF ( ln_spmfb ) THEN ! Feedback file format
1610
1611	DO jset = 1, jnumspmfb
1612
1613	nspmsets = nspmsets + 1
1614
1615	CALL obs_rea_spm( 0, spmdata(nspmsets), 1, &
1616	& spmfbfiles(jset:jset), &
1617	& nspmvars, nspmextr, nitend-nit000+2, &
1618	& dobsini, dobsend, ln_ignmis, .FALSE. )
1619
1620	CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1621	& ln_spm, ln_nea )
1622
1623	ENDDO
1624
1625	ELSE ! Original file format
1626
1627	nspmsets = nspmsets + 1
1628
1629	CALL obs_rea_spm( 1, spmdata(nspmsets), jnumspm, &
1630	& spmfiles(1:jnumspm), &
1631	& nspmvars, nspmextr, nitend-nit000+2, &
1632	& dobsini, dobsend, ln_ignmis, .FALSE. )
1633
1634	CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1635	& ln_spm, ln_nea )
1636
1637	ENDIF
1638
1639	ENDIF
1640
1641	! - fco2
1642
1643	IF ( ln_fco2 ) THEN
1644
1645	! Set the number of variables for fco2 to 1
1646	nfco2vars = 1
1647
1648	! Set the number of extra variables for fco2 to 0
1649	nfco2extr = 0
1650
1651	IF ( ln_fco2fb ) THEN
1652	nfco2sets = jnumfco2fb
1653	ELSE
1654	nfco2sets = 1
1655	ENDIF
1656
1657	ALLOCATE(fco2data(nfco2sets))
1658	ALLOCATE(fco2datqc(nfco2sets))
1659	fco2data(:)%nsurf=0
1660	fco2datqc(:)%nsurf=0
1661
1662	nfco2sets = 0
1663
1664	IF ( ln_fco2fb ) THEN ! Feedback file format
1665
1666	DO jset = 1, jnumfco2fb
1667
1668	nfco2sets = nfco2sets + 1
1669
1670	CALL obs_rea_fco2( 0, fco2data(nfco2sets), 1, &
1671	& fco2fbfiles(jset:jset), &
1672	& nfco2vars, nfco2extr, nitend-nit000+2, &
1673	& dobsini, dobsend, ln_ignmis, .FALSE. )
1674
1675	CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
1676	& ln_fco2, ln_nea )
1677
1678	ENDDO
1679
1680	ELSE ! Original file format
1681
1682	nfco2sets = nfco2sets + 1
1683
1684	CALL obs_rea_fco2( 1, fco2data(nfco2sets), jnumfco2, &
1685	& fco2files(1:jnumfco2), &
1686	& nfco2vars, nfco2extr, nitend-nit000+2, &
1687	& dobsini, dobsend, ln_ignmis, .FALSE. )
1688
1689	CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
1690	& ln_fco2, ln_nea )
1691
1692	ENDIF
1693
1694	ENDIF
1695
1696	! - pco2
1697
1698	IF ( ln_pco2 ) THEN
1699
1700	! Set the number of variables for pco2 to 1
1701	npco2vars = 1
1702
1703	! Set the number of extra variables for pco2 to 0
1704	npco2extr = 0
1705
1706	IF ( ln_pco2fb ) THEN
1707	npco2sets = jnumpco2fb
1708	ELSE
1709	npco2sets = 1
1710	ENDIF
1711
1712	ALLOCATE(pco2data(npco2sets))
1713	ALLOCATE(pco2datqc(npco2sets))
1714	pco2data(:)%nsurf=0
1715	pco2datqc(:)%nsurf=0
1716
1717	npco2sets = 0
1718
1719	IF ( ln_pco2fb ) THEN ! Feedback file format
1720
1721	DO jset = 1, jnumpco2fb
1722
1723	npco2sets = npco2sets + 1
1724
1725	CALL obs_rea_pco2( 0, pco2data(npco2sets), 1, &
1726	& pco2fbfiles(jset:jset), &
1727	& npco2vars, npco2extr, nitend-nit000+2, &
1728	& dobsini, dobsend, ln_ignmis, .FALSE. )
1729
1730	CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
1731	& ln_pco2, ln_nea )
1732
1733	ENDDO
1734
1735	ELSE ! Original file format
1736
1737	npco2sets = npco2sets + 1
1738
1739	CALL obs_rea_pco2( 1, pco2data(npco2sets), jnumpco2, &
1740	& pco2files(1:jnumpco2), &
1741	& npco2vars, npco2extr, nitend-nit000+2, &
1742	& dobsini, dobsend, ln_ignmis, .FALSE. )
1743
1744	CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
1745	& ln_pco2, ln_nea )
1746
1747	ENDIF
1748
1749	ENDIF
1750
1751	END SUBROUTINE dia_obs_init
1752
1753	SUBROUTINE dia_obs( kstp )
1754	!!----------------------------------------------------------------------
1755	!! * ROUTINE dia_obs *
1756	!!
1757	!! ** Purpose : Call the observation operators on each time step
1758	!!
1759	!! ** Method : Call the observation operators on each time step to
1760	!! compute the model equivalent of the following date:
1761	!! - T profiles
1762	!! - S profiles
1763	!! - Sea surface height (referenced to a mean)
1764	!! - Sea surface temperature
1765	!! - Sea surface salinity
1766	!! - Velocity component (U,V) profiles
1767	!! - Sea surface log10(chlorophyll)
1768	!! - Sea surface spm
1769	!! - Sea surface fco2
1770	!! - Sea surface pco2
1771	!!
1772	!! ** Action :
1773	!!
1774	!! History :
1775	!! ! 06-03 (K. Mogensen) Original code
1776	!! ! 06-05 (K. Mogensen) Reformatted
1777	!! ! 06-10 (A. Weaver) Cleaning
1778	!! ! 07-03 (K. Mogensen) General handling of profiles
1779	!! ! 07-04 (G. Smith) Generalized surface operators
1780	!! ! 08-10 (M. Valdivieso) obs operator for velocity profiles
1781	!! ! 14-08 (J. While) observation operator for profiles in
1782	!! generalised vertical coordinates
1783	!!----------------------------------------------------------------------
1784	!! * Modules used
1785	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1786	& rdt, &
1787	& gdept_1d, &
1788	#if defined key_vvl
1789	& gdept_n, &
1790	#else
1791	& gdept_1d, &
1792	#endif
1793	& tmask, umask, vmask
1794	USE phycst, ONLY : & ! Physical constants
1795	& rday
1796	USE oce, ONLY : & ! Ocean dynamics and tracers variables
1797	& tsn, &
1798	& un, vn, &
1799	& sshn
1800	#if defined key_lim3
1801	USE ice, ONLY : & ! LIM Ice model variables
1802	& frld
1803	#endif
1804	#if defined key_lim2
1805	USE ice_2, ONLY : & ! LIM Ice model variables
1806	& frld
1807	#endif
1808	#if defined key_hadocc
1809	USE trc, ONLY : & ! HadOCC chlorophyll, fCO2 and pCO2
1810	& HADOCC_CHL, &
1811	& HADOCC_FCO2, &
1812	& HADOCC_PCO2, &
1813	& HADOCC_FILL_FLT
1814	#elif defined key_medusa && defined key_foam_medusa
1815	USE trc, ONLY : & ! MEDUSA chlorophyll, fCO2 and pCO2
1816	& MEDUSA_CHL, &
1817	& MEDUSA_FCO2, &
1818	& MEDUSA_PCO2, &
1819	& MEDUSA_FILL_FLT
1820	#elif defined key_fabm
1821	USE fabm
1822	USE par_fabm
1823	#endif
1824	#if defined key_spm
1825	USE par_spm, ONLY: & ! ERSEM/SPM sediments
1826	& jp_spm
1827	#endif
1828	IMPLICIT NONE
1829
1830	!! * Arguments
1831	INTEGER, INTENT(IN) :: kstp ! Current timestep
1832	!! * Local declarations
1833	INTEGER :: idaystp ! Number of timesteps per day
1834	INTEGER :: jprofset ! Profile data set loop variable
1835	INTEGER :: jslaset ! SLA data set loop variable
1836	INTEGER :: jsstset ! SST data set loop variable
1837	INTEGER :: jseaiceset ! sea ice data set loop variable
1838	INTEGER :: jveloset ! velocity profile data loop variable
1839	INTEGER :: jlogchlset ! logchl data set loop variable
1840	INTEGER :: jspmset ! spm data set loop variable
1841	INTEGER :: jfco2set ! fco2 data set loop variable
1842	INTEGER :: jpco2set ! pco2 data set loop variable
1843	INTEGER :: jvar ! Variable number
1844	#if ! defined key_lim2 && ! defined key_lim3
1845	REAL(wp), POINTER, DIMENSION(:,:) :: frld
1846	#endif
1847	REAL(wp) :: tiny ! small number
1848	REAL(wp), DIMENSION(jpi,jpj) :: &
1849	logchl ! array for log chlorophyll
1850
1851
1852
1853
1854
1855
1856
1857	! JOZEF ADDITION
1858
1859
1860	REAL(wp), DIMENSION(jpi,jpj,nn_logchlpftscc) :: &
1861	logchlpfts
1862
1863
1864	! END
1865
1866
1867
1868
1869
1870
1871
1872	! array for log chlorophyll pfts
1873	REAL(wp), DIMENSION(jpi,jpj) :: &
1874	maskchl ! array for special chlorophyll mask
1875	REAL(wp), DIMENSION(jpi,jpj) :: &
1876	spm ! array for spm
1877	REAL(wp), DIMENSION(jpi,jpj) :: &
1878	fco2 ! array for fco2
1879	REAL(wp), DIMENSION(jpi,jpj) :: &
1880	maskfco2 ! array for special fco2 mask
1881	REAL(wp), DIMENSION(jpi,jpj) :: &
1882	pco2 ! array for pco2
1883	REAL(wp), DIMENSION(jpi,jpj) :: &
1884	maskpco2 ! array for special pco2 mask
1885	INTEGER :: jn ! loop index
1886	#if defined key_fabm
1887	REAL(wp), DIMENSION(jpi,jpj,jpk) :: logchl_3d
1888
1889
1890
1891
1892
1893
1894
1895	! JOZEF ADDITION
1896
1897	REAL(wp), DIMENSION(jpi,jpj,jpk,nn_logchlpftscc) :: logchlpfts_3d
1898	#endif
1899
1900	INTEGER :: pft
1901	! END
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1912
1913	#if ! defined key_lim2 && ! defined key_lim3
1914	CALL wrk_alloc(jpi,jpj,frld)
1915	#endif
1916
1917	IF(lwp) THEN
1918	WRITE(numout,*)
1919	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1920	WRITE(numout,*) '~~~~~~~'
1921	ENDIF
1922
1923	idaystp = NINT( rday / rdt )
1924
1925	!-----------------------------------------------------------------------
1926	! No LIM => frld == 0.0_wp
1927	!-----------------------------------------------------------------------
1928	#if ! defined key_lim2 && ! defined key_lim3
1929	frld(:,:) = 0.0_wp
1930	#endif
1931	!-----------------------------------------------------------------------
1932	! Depending on switches call various observation operators
1933	!-----------------------------------------------------------------------
1934
1935	! - Temperature/salinity profiles
1936	IF ( ln_t3d .OR. ln_s3d ) THEN
1937	DO jprofset = 1, nprofsets
1938	IF( (.NOT. lk_vvl) .AND. (ln_zco .OR. ln_zps) ) THEN
1939	IF(lwp) THEN
1940	WRITE(numout,*) 'dia_obs : calling obs_pro_opt'
1941	ENDIF
1942	IF ( ld_enact(jprofset) ) THEN
1943	CALL obs_pro_opt( prodatqc(jprofset), &
1944	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1945	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1946	& gdept_1d, tmask, n1dint, n2dint, &
1947	& kdailyavtypes = endailyavtypes )
1948	ELSE
1949	CALL obs_pro_opt( prodatqc(jprofset), &
1950	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1951	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1952	& gdept_1d, tmask, n1dint, n2dint )
1953	ENDIF
1954	ELSE
1955	IF(lwp) THEN
1956	WRITE(numout,*) 'dia_obs : calling obs_pro_sco_opt'
1957	ENDIF
1958	IF ( ld_enact(jprofset) ) THEN
1959	CALL obs_pro_sco_opt( prodatqc(jprofset), &
1960	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1961	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1962	& fsdept(:,:,:), fsdepw(:,:,:), &
1963	& tmask, n1dint, n2dint, &
1964	& kdailyavtypes = endailyavtypes )
1965	ELSE
1966	CALL obs_pro_sco_opt( prodatqc(jprofset), &
1967	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1968	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1969	& fsdept(:,:,:), fsdepw(:,:,:), &
1970	& tmask, n1dint, n2dint )
1971	ENDIF
1972	ENDIF
1973	END DO
1974	ENDIF
1975
1976	! - Sea surface anomaly
1977	IF ( ln_sla ) THEN
1978	DO jslaset = 1, nslasets
1979	CALL obs_sla_opt( sladatqc(jslaset), &
1980	& kstp, jpi, jpj, nit000, sshn, &
1981	& tmask(:,:,1), n2dint )
1982	END DO
1983	ENDIF
1984
1985	! - Sea surface temperature
1986	IF ( ln_sst ) THEN
1987	DO jsstset = 1, nsstsets
1988	CALL obs_sst_opt( sstdatqc(jsstset), &
1989	& kstp, jpi, jpj, nit000, idaystp, &
1990	& tsn(:,:,1,jp_tem), tmask(:,:,1), &
1991	& n2dint, ld_sstnight(jsstset) )
1992	END DO
1993	ENDIF
1994
1995	! - Sea surface salinity
1996	IF ( ln_sss ) THEN
1997	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1998	ENDIF
1999
2000	#if defined key_lim2 \|\| defined key_lim3
2001	IF ( ln_seaice ) THEN
2002	DO jseaiceset = 1, nseaicesets
2003	CALL obs_seaice_opt( seaicedatqc(jseaiceset), &
2004	& kstp, jpi, jpj, nit000, 1.-frld, &
2005	& tmask(:,:,1), n2dint )
2006	END DO
2007	ENDIF
2008	#endif
2009
2010	! - Velocity profiles
2011	IF ( ln_vel3d ) THEN
2012	DO jveloset = 1, nvelosets
2013	! zonal component of velocity
2014	CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
2015	& nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
2016	n1dint, n2dint, ld_velav(jveloset) )
2017	END DO
2018	ENDIF
2019
2020	IF ( ln_logchl ) THEN
2021
2022	#if defined key_hadocc
2023	logchl(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
2024	#elif defined key_medusa && defined key_foam_medusa
2025	logchl(:,:) = MEDUSA_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
2026	#elif defined key_fabm
2027	logchl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
2028	logchl(:,:) = logchl_3d(:,:,1)
2029	#else
2030	CALL ctl_stop( ' Trying to run logchl observation operator', &
2031	& ' but no biogeochemical model appears to have been defined' )
2032	#endif
2033
2034	maskchl(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
2035
2036	! Take the log10 where we can, otherwise exclude
2037	tiny = 1.0e-20
2038	WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
2039	logchl(:,:) = LOG10(logchl(:,:))
2040	ELSEWHERE
2041	logchl(:,:) = obfillflt
2042	maskchl(:,:) = 0
2043	END WHERE
2044
2045	DO jlogchlset = 1, nlogchlsets
2046	CALL obs_logchl_opt( logchldatqc(jlogchlset), &
2047	& kstp, jpi, jpj, nit000, logchl(:,:), &
2048	& maskchl(:,:), n2dint )
2049	END DO
2050	ENDIF
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060	! JOZEF ADDITION
2061
2062
2063
2064	IF ( ln_logchlpft ) THEN
2065
2066	#if defined key_hadocc
2067	logchl(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
2068	#elif defined key_medusa && defined key_foam_medusa
2069	logchl(:,:) = MEDUSA_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
2070	#elif defined key_fabm
2071	logchlpfts_3d(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1)
2072	logchlpfts(:,:,1) = logchlpfts_3d(:,:,1,1)
2073	logchlpfts_3d(:,:,:,2) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl2)
2074	logchlpfts(:,:,2) = logchlpfts_3d(:,:,1,2)
2075	logchlpfts_3d(:,:,:,3) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl3)
2076	logchlpfts(:,:,3) = logchlpfts_3d(:,:,1,3)
2077	logchlpfts_3d(:,:,:,4) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
2078	logchlpfts(:,:,4) = logchlpfts_3d(:,:,1,4)
2079	#else
2080	CALL ctl_stop( ' Trying to run logchl observation operator', &
2081	& ' but no biogeochemical model appears to have been defined' )
2082	#endif
2083
2084	! create a special mask to exclude certain things
2085
2086	! Take the log10 where we can, otherwise exclude
2087	tiny = 1.0e-20
2088
2089	DO pft=1,nn_logchlpftscc
2090
2091	maskchl(:,:) = tmask(:,:,1)
2092
2093	WRITE(numout,*) 'Running pft: ',pft
2094
2095	WHERE(logchlpfts(:,:,pft) > tiny .AND. logchlpfts(:,:,pft) /= obfillflt )
2096	logchlpfts(:,:,pft) = LOG10(logchlpfts(:,:,pft))
2097	ELSEWHERE
2098	logchlpfts(:,:,pft) = obfillflt
2099	maskchl(:,:) = 0
2100	END WHERE
2101
2102	WRITE(numout,*) 'Max value of Pft: ', MINVAL(logchlpfts(:,:,pft))
2103
2104	DO jlogchlset = 1, nlogchlpftsets
2105	WRITE(numout,*) 'nlogchlpftsets, nn_logchlpftscc, obfillflt = ', nlogchlpftsets, nn_logchlpftscc, obfillflt
2106	WRITE(numout,*) 'kstp, jpi, jpj, nit000, n2dint = ', kstp, jpi, jpj, nit000, n2dint
2107	!WRITE(numout,*) 'ALLOCATED(logchlpftdatqc) = ', ALLOCATED(logchlpftdatqc)
2108	WRITE(numout,*) 'SHAPE(logchlpftdatqc) = ', SHAPE(logchlpftdatqc)
2109	WRITE(numout,*) 'SIZE(logchlpftdatqc) = ', SIZE(logchlpftdatqc)
2110	!WRITE(numout,*) 'ALLOCATED(logchlpfts) = ', ALLOCATED(logchlpfts)
2111	WRITE(numout,*) 'SHAPE(logchlpfts) = ', SHAPE(logchlpfts)
2112	WRITE(numout,*) 'SIZE(logchlpfts) = ', SIZE(logchlpfts)
2113	!WRITE(numout,*) 'ALLOCATED(maskchl) = ', ALLOCATED(maskchl)
2114	WRITE(numout,*) 'SHAPE(maskchl) = ', SHAPE(maskchl)
2115	WRITE(numout,*) 'SIZE(maskchl) = ', SIZE(maskchl)
2116	CALL flush(numout)
2117
2118	CALL obs_logchl_opt( logchlpftdatqc(pft), &
2119	& kstp, jpi, jpj, nit000, logchlpfts(:,:,pft), &
2120	& maskchl(:,:), n2dint )
2121
2122	WRITE(numout,*) 'Max value of Pft: ', MINVAL(logchlpfts(:,:,pft))
2123	END DO
2124
2125
2126	END DO
2127
2128	ENDIF
2129
2130	! END
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147	IF ( ln_spm ) THEN
2148	#if defined key_spm
2149	spm(:,:) = 0.0
2150	DO jn = 1, jp_spm
2151	spm(:,:) = spm(:,:) + trn(:,:,1,jn) ! sum SPM sizes
2152	END DO
2153	#else
2154	CALL ctl_stop( ' Trying to run spm observation operator', &
2155	& ' but no spm model appears to have been defined' )
2156	#endif
2157
2158	DO jspmset = 1, nspmsets
2159	CALL obs_spm_opt( spmdatqc(jspmset), &
2160	& kstp, jpi, jpj, nit000, spm(:,:), &
2161	& tmask(:,:,1), n2dint )
2162	END DO
2163	ENDIF
2164
2165	IF ( ln_fco2 ) THEN
2166	maskfco2(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
2167	#if defined key_hadocc
2168	fco2(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC
2169	IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
2170	fco2(:,:) = obfillflt
2171	maskfco2(:,:) = 0
2172	CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
2173	& ' on timestep ' // TRIM(STR(kstp)), &
2174	& ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
2175	ENDIF
2176	#elif defined key_medusa && defined key_foam_medusa
2177	fco2(:,:) = MEDUSA_FCO2(:,:) ! fCO2 from MEDUSA
2178	IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
2179	fco2(:,:) = obfillflt
2180	maskfco2(:,:) = 0
2181	CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
2182	& ' on timestep ' // TRIM(STR(kstp)), &
2183	& ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
2184	ENDIF
2185	#elif defined key_fabm
2186	! First, get pCO2 from FABM
2187	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
2188	pco2(:,:) = pco2_3d(:,:,1)
2189	! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
2190	! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
2191	! and data reduction routines, Deep-Sea Research II, 56: 512-522.
2192	! and
2193	! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
2194	! Marine Chemistry, 2: 203-215.
2195	! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
2196	! not explicitly included - atmospheric pressure is not necessarily available so this is
2197	! the best assumption.
2198	! Further, the (1-xCO2)^2 term has been neglected. This is common practice
2199	! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
2200	! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
2201	! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
2202	fco2(:,:) = pco2(:,:) * EXP((-1636.75 + &
2203	& 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - &
2204	& 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
2205	& 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0) + &
2206	& 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / &
2207	& (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
2208	#else
2209	CALL ctl_stop( ' Trying to run fco2 observation operator', &
2210	& ' but no biogeochemical model appears to have been defined' )
2211	#endif
2212
2213	DO jfco2set = 1, nfco2sets
2214	CALL obs_fco2_opt( fco2datqc(jfco2set), &
2215	& kstp, jpi, jpj, nit000, fco2(:,:), &
2216	& maskfco2(:,:), n2dint )
2217	END DO
2218	ENDIF
2219
2220	IF ( ln_pco2 ) THEN
2221	maskpco2(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
2222	#if defined key_hadocc
2223	pco2(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC
2224	IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
2225	pco2(:,:) = obfillflt
2226	maskpco2(:,:) = 0
2227	CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
2228	& ' on timestep ' // TRIM(STR(kstp)), &
2229	& ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
2230	ENDIF
2231	#elif defined key_medusa && defined key_foam_medusa
2232	pco2(:,:) = MEDUSA_PCO2(:,:) ! pCO2 from MEDUSA
2233	IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
2234	pco2(:,:) = obfillflt
2235	maskpco2(:,:) = 0
2236	CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
2237	& ' on timestep ' // TRIM(STR(kstp)), &
2238	& ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
2239	ENDIF
2240	#elif defined key_fabm
2241	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
2242	pco2(:,:) = pco2_3d(:,:,1)
2243	#else
2244	CALL ctl_stop( ' Trying to run pCO2 observation operator', &
2245	& ' but no biogeochemical model appears to have been defined' )
2246	#endif
2247
2248	DO jpco2set = 1, npco2sets
2249	CALL obs_pco2_opt( pco2datqc(jpco2set), &
2250	& kstp, jpi, jpj, nit000, pco2(:,:), &
2251	& maskpco2(:,:), n2dint )
2252	END DO
2253	ENDIF
2254
2255	#if ! defined key_lim2 && ! defined key_lim3
2256	CALL wrk_dealloc(jpi,jpj,frld)
2257	#endif
2258
2259	END SUBROUTINE dia_obs
2260
2261	SUBROUTINE dia_obs_wri
2262	!!----------------------------------------------------------------------
2263	!! * ROUTINE dia_obs_wri *
2264	!!
2265	!! ** Purpose : Call observation diagnostic output routines
2266	!!
2267	!! ** Method : Call observation diagnostic output routines
2268	!!
2269	!! ** Action :
2270	!!
2271	!! History :
2272	!! ! 06-03 (K. Mogensen) Original code
2273	!! ! 06-05 (K. Mogensen) Reformatted
2274	!! ! 06-10 (A. Weaver) Cleaning
2275	!! ! 07-03 (K. Mogensen) General handling of profiles
2276	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
2277	!!----------------------------------------------------------------------
2278	IMPLICIT NONE
2279
2280	!! * Local declarations
2281
2282	INTEGER :: jprofset ! Profile data set loop variable
2283	INTEGER :: jveloset ! Velocity data set loop variable
2284	INTEGER :: jslaset ! SLA data set loop variable
2285	INTEGER :: jsstset ! SST data set loop variable
2286	INTEGER :: jseaiceset ! Sea Ice data set loop variable
2287	INTEGER :: jlogchlset ! logchl data set loop variable
2288	INTEGER :: jspmset ! spm data set loop variable
2289	INTEGER :: jfco2set ! fco2 data set loop variable
2290	INTEGER :: jpco2set ! pco2 data set loop variable
2291	INTEGER :: jset
2292	INTEGER :: jfbini
2293	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
2294	CHARACTER(LEN=20) :: cdtmp
2295	INTEGER :: pft
2296	!-----------------------------------------------------------------------
2297	! Depending on switches call various observation output routines
2298	!-----------------------------------------------------------------------
2299
2300	! - Temperature/salinity profiles
2301
2302	IF( ln_t3d .OR. ln_s3d ) THEN
2303
2304	! Copy data from prodatqc to profdata structures
2305	DO jprofset = 1, nprofsets
2306
2307	CALL obs_prof_decompress( prodatqc(jprofset), &
2308	& profdata(jprofset), .TRUE., numout )
2309
2310	END DO
2311
2312	! Write the profiles.
2313
2314	jprofset = 0
2315
2316	! ENACT insitu data
2317
2318	IF ( ln_ena ) THEN
2319
2320	jprofset = jprofset + 1
2321
2322	CALL obs_wri_p3d( 'enact', profdata(jprofset) )
2323
2324	ENDIF
2325
2326	! Coriolis insitu data
2327
2328	IF ( ln_cor ) THEN
2329
2330	jprofset = jprofset + 1
2331
2332	CALL obs_wri_p3d( 'corio', profdata(jprofset) )
2333
2334	ENDIF
2335
2336	! Feedback insitu data
2337
2338	IF ( ln_profb ) THEN
2339
2340	jfbini = jprofset + 1
2341
2342	DO jprofset = jfbini, nprofsets
2343
2344	jset = jprofset - jfbini + 1
2345	WRITE(cdtmp,'(A,I2.2)')'profb_',jset
2346	CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
2347
2348	END DO
2349
2350	ENDIF
2351
2352	ENDIF
2353
2354	! - Sea surface anomaly
2355	IF ( ln_sla ) THEN
2356
2357	! Copy data from sladatqc to sladata structures
2358	DO jslaset = 1, nslasets
2359
2360	CALL obs_surf_decompress( sladatqc(jslaset), &
2361	& sladata(jslaset), .TRUE., numout )
2362
2363	END DO
2364
2365	jslaset = 0
2366
2367	! Write the AVISO SLA data
2368
2369	IF ( ln_sladt ) THEN
2370
2371	jslaset = 1
2372	CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
2373	jslaset = 2
2374	CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
2375
2376	ENDIF
2377
2378	IF ( ln_slafb ) THEN
2379
2380	jfbini = jslaset + 1
2381
2382	DO jslaset = jfbini, nslasets
2383
2384	jset = jslaset - jfbini + 1
2385	WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
2386	CALL obs_wri_sla( cdtmp, sladata(jslaset) )
2387
2388	END DO
2389
2390	ENDIF
2391
2392	ENDIF
2393
2394	! - Sea surface temperature
2395	IF ( ln_sst ) THEN
2396
2397	! Copy data from sstdatqc to sstdata structures
2398	DO jsstset = 1, nsstsets
2399
2400	CALL obs_surf_decompress( sstdatqc(jsstset), &
2401	& sstdata(jsstset), .TRUE., numout )
2402
2403	END DO
2404
2405	jsstset = 0
2406
2407	! Write the AVISO SST data
2408
2409	IF ( ln_reysst ) THEN
2410
2411	jsstset = jsstset + 1
2412	CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
2413
2414	ENDIF
2415
2416	IF ( ln_ghrsst ) THEN
2417
2418	jsstset = jsstset + 1
2419	CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
2420
2421	ENDIF
2422
2423	IF ( ln_sstfb ) THEN
2424
2425	jfbini = jsstset + 1
2426
2427	DO jsstset = jfbini, nsstsets
2428
2429	jset = jsstset - jfbini + 1
2430	WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
2431	CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
2432
2433	END DO
2434
2435	ENDIF
2436
2437	ENDIF
2438
2439	! - Sea surface salinity
2440	IF ( ln_sss ) THEN
2441	IF(lwp) WRITE(numout,*) ' SSS currently not available'
2442	ENDIF
2443
2444	! - Sea Ice Concentration
2445	IF ( ln_seaice ) THEN
2446
2447	! Copy data from seaicedatqc to seaicedata structures
2448	DO jseaiceset = 1, nseaicesets
2449
2450	CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
2451	& seaicedata(jseaiceset), .TRUE., numout )
2452
2453	END DO
2454
2455	! Write the Sea Ice data
2456	DO jseaiceset = 1, nseaicesets
2457
2458	WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
2459	CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
2460
2461	END DO
2462
2463	ENDIF
2464
2465	! Velocity data
2466	IF( ln_vel3d ) THEN
2467
2468	! Copy data from veldatqc to velodata structures
2469	DO jveloset = 1, nvelosets
2470
2471	CALL obs_prof_decompress( veldatqc(jveloset), &
2472	& velodata(jveloset), .TRUE., numout )
2473
2474	END DO
2475
2476	! Write the profiles.
2477
2478	jveloset = 0
2479
2480	! Daily averaged data
2481
2482	IF ( ln_velavcur ) THEN
2483
2484	jveloset = jveloset + 1
2485
2486	CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
2487
2488	ENDIF
2489
2490	! High frequency data
2491
2492	IF ( ln_velhrcur ) THEN
2493
2494	jveloset = jveloset + 1
2495
2496	CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
2497
2498	ENDIF
2499
2500	! Daily averaged data
2501
2502	IF ( ln_velavadcp ) THEN
2503
2504	jveloset = jveloset + 1
2505
2506	CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
2507
2508	ENDIF
2509
2510	! High frequency data
2511
2512	IF ( ln_velhradcp ) THEN
2513
2514	jveloset = jveloset + 1
2515
2516	CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
2517
2518	ENDIF
2519
2520	! Feedback velocity data
2521
2522	IF ( ln_velfb ) THEN
2523
2524	jfbini = jveloset + 1
2525
2526	DO jveloset = jfbini, nvelosets
2527
2528	jset = jveloset - jfbini + 1
2529	WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
2530	CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
2531
2532	END DO
2533
2534	ENDIF
2535
2536	ENDIF
2537
2538	! - log10(chlorophyll)
2539	IF ( ln_logchl ) THEN
2540
2541	! Copy data from logchldatqc to logchldata structures
2542	DO jlogchlset = 1, nlogchlsets
2543
2544	CALL obs_surf_decompress( logchldatqc(jlogchlset), &
2545	& logchldata(jlogchlset), .TRUE., numout )
2546
2547	END DO
2548
2549	! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2550	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2551	DO jlogchlset = 1, nlogchlsets
2552	WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
2553	logchldata(jlogchlset)%nqc(:) = 1
2554	logchldata(jlogchlset)%robs(:,1) = obfillflt
2555	END WHERE
2556	END DO
2557
2558	! Write the logchl data
2559	DO jlogchlset = 1, nlogchlsets
2560
2561	WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
2562	CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
2563
2564	END DO
2565
2566	ENDIF
2567
2568
2569
2570
2571
2572
2573
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583	! JOZEF ADDITION
2584
2585	IF ( ln_logchlpft ) THEN
2586
2587	WRITE(numout,*)'Pfts loop started'
2588	WRITE(numout,*)'nlogchlpftsets', nlogchlpftsets
2589
2590	DO pft = 1, nn_logchlpftscc
2591
2592	! Copy data from logchldatqc to logchldata structures
2593	DO jlogchlset = 1, nlogchlpftsets
2594
2595	CALL obs_surf_decompress( logchlpftdatqc(pft), &
2596	& logchlpftdata(pft), .TRUE., numout )
2597
2598	END DO
2599
2600	! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2601	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2602	DO jlogchlset = 1, nlogchlpftsets
2603	WHERE ( logchlpftdata(pft)%rmod(:,1) == obfillflt )
2604	logchlpftdata(pft)%nqc(:) = 1
2605	logchlpftdata(pft)%robs(:,1) = obfillflt
2606	END WHERE
2607	END DO
2608
2609	! Write the logchl data
2610	DO jlogchlset = 1, nlogchlpftsets
2611
2612	WRITE(numout,*)'Code Pfts here!'
2613
2614	! WRITE(cdtmp,'(A,I2.2)')'logchlpftfb_',jlogchlset
2615	WRITE(cdtmp,'(A,I2.2,A,I2.2)')'logchlpft',pft,'fb_',jlogchlset
2616	!CALL obs_wri_logchl( cdtmp, logchlpftdata(jlogchlset,pft) )
2617	CALL obs_wri_logchl( cdtmp, logchlpftdata(pft), pft_num=jlogchlset )
2618
2619	END DO
2620
2621	END DO
2622
2623	ENDIF
2624
2625	! END
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636
2637
2638
2639
2640
2641
2642
2643
2644
2645	! - spm
2646	IF ( ln_spm ) THEN
2647
2648	! Copy data from spmdatqc to spmdata structures
2649	DO jspmset = 1, nspmsets
2650
2651	CALL obs_surf_decompress( spmdatqc(jspmset), &
2652	& spmdata(jspmset), .TRUE., numout )
2653
2654	END DO
2655
2656	! Write the spm data
2657	DO jspmset = 1, nspmsets
2658
2659	WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset
2660	CALL obs_wri_spm( cdtmp, spmdata(jspmset) )
2661
2662	END DO
2663
2664	ENDIF
2665
2666	! - fco2
2667	IF ( ln_fco2 ) THEN
2668
2669	! Copy data from fco2datqc to fco2data structures
2670	DO jfco2set = 1, nfco2sets
2671
2672	CALL obs_surf_decompress( fco2datqc(jfco2set), &
2673	& fco2data(jfco2set), .TRUE., numout )
2674
2675	END DO
2676
2677	! Mark as bad observations with no valid model counterpart due to fCO2 not being in the restart
2678	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2679	DO jfco2set = 1, nfco2sets
2680	WHERE ( fco2data(jfco2set)%rmod(:,1) == obfillflt )
2681	fco2data(jfco2set)%nqc(:) = 1
2682	fco2data(jfco2set)%robs(:,1) = obfillflt
2683	END WHERE
2684	END DO
2685
2686	! Write the fco2 data
2687	DO jfco2set = 1, nfco2sets
2688
2689	WRITE(cdtmp,'(A,I2.2)')'fco2fb_',jfco2set
2690	CALL obs_wri_fco2( cdtmp, fco2data(jfco2set) )
2691
2692	END DO
2693
2694	ENDIF
2695
2696	! - pco2
2697	IF ( ln_pco2 ) THEN
2698
2699	! Copy data from pco2datqc to pco2data structures
2700	DO jpco2set = 1, npco2sets
2701
2702	CALL obs_surf_decompress( pco2datqc(jpco2set), &
2703	& pco2data(jpco2set), .TRUE., numout )
2704
2705	END DO
2706
2707	! Mark as bad observations with no valid model counterpart due to pco2 not being in the restart
2708	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2709	DO jpco2set = 1, npco2sets
2710	WHERE ( pco2data(jpco2set)%rmod(:,1) == obfillflt )
2711	pco2data(jpco2set)%nqc(:) = 1
2712	pco2data(jpco2set)%robs(:,1) = obfillflt
2713	END WHERE
2714	END DO
2715
2716	! Write the pco2 data
2717	DO jpco2set = 1, npco2sets
2718
2719	WRITE(cdtmp,'(A,I2.2)')'pco2fb_',jpco2set
2720	CALL obs_wri_pco2( cdtmp, pco2data(jpco2set) )
2721
2722	END DO
2723
2724	ENDIF
2725
2726	END SUBROUTINE dia_obs_wri
2727
2728	SUBROUTINE dia_obs_dealloc
2729	IMPLICIT NONE
2730	!!----------------------------------------------------------------------
2731	!! * ROUTINE dia_obs_dealloc *
2732	!!
2733	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
2734	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
2735	!!
2736	!! ** Method : Clean up various arrays left behind by the obs_oper.
2737	!!
2738	!! ** Action :
2739	!!
2740	!!----------------------------------------------------------------------
2741	!! obs_grid deallocation
2742	CALL obs_grid_deallocate
2743
2744	!! diaobs deallocation
2745	IF ( nprofsets > 0 ) THEN
2746	DEALLOCATE(ld_enact, &
2747	& profdata, &
2748	& prodatqc)
2749	END IF
2750	IF ( ln_sla ) THEN
2751	DEALLOCATE(sladata, &
2752	& sladatqc)
2753	END IF
2754	IF ( ln_seaice ) THEN
2755	DEALLOCATE(sladata, &
2756	& sladatqc)
2757	END IF
2758	IF ( ln_sst ) THEN
2759	DEALLOCATE(sstdata, &
2760	& sstdatqc)
2761	END IF
2762	IF ( ln_vel3d ) THEN
2763	DEALLOCATE(ld_velav, &
2764	& velodata, &
2765	& veldatqc)
2766	END IF
2767	END SUBROUTINE dia_obs_dealloc
2768
2769	SUBROUTINE ini_date( ddobsini )
2770	!!----------------------------------------------------------------------
2771	!! * ROUTINE ini_date *
2772	!!
2773	!! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
2774	!!
2775	!! ** Method : Get initial data in double precision YYYYMMDD.HHMMSS format
2776	!!
2777	!! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format
2778	!!
2779	!! History :
2780	!! ! 06-03 (K. Mogensen) Original code
2781	!! ! 06-05 (K. Mogensen) Reformatted
2782	!! ! 06-10 (A. Weaver) Cleaning
2783	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
2784	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
2785	!!----------------------------------------------------------------------
2786	USE phycst, ONLY : & ! Physical constants
2787	& rday
2788	! USE daymod, ONLY : & ! Time variables
2789	! & nmonth_len
2790	USE dom_oce, ONLY : & ! Ocean space and time domain variables
2791	& rdt
2792
2793	IMPLICIT NONE
2794
2795	!! * Arguments
2796	REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
2797
2798	!! * Local declarations
2799	INTEGER :: iyea ! date - (year, month, day, hour, minute)
2800	INTEGER :: imon
2801	INTEGER :: iday
2802	INTEGER :: ihou
2803	INTEGER :: imin
2804	INTEGER :: imday ! Number of days in month.
2805	REAL(KIND=wp) :: zdayfrc ! Fraction of day
2806
2807	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
2808
2809	!!----------------------------------------------------------------------
2810	!! Initial date initialization (year, month, day, hour, minute)
2811	!! (This assumes that the initial date is for 00z))
2812	!!----------------------------------------------------------------------
2813	iyea = ndate0 / 10000
2814	imon = ( ndate0 - iyea * 10000 ) / 100
2815	iday = ndate0 - iyea * 10000 - imon * 100
2816	ihou = 0
2817	imin = 0
2818
2819	!!----------------------------------------------------------------------
2820	!! Compute number of days + number of hours + min since initial time
2821	!!----------------------------------------------------------------------
2822	iday = iday + ( nit000 -1 ) * rdt / rday
2823	zdayfrc = ( nit000 -1 ) * rdt / rday
2824	zdayfrc = zdayfrc - aint(zdayfrc)
2825	ihou = int( zdayfrc * 24 )
2826	imin = int( (zdayfrc * 24 - ihou) * 60 )
2827
2828	!!-----------------------------------------------------------------------
2829	!! Convert number of days (iday) into a real date
2830	!!----------------------------------------------------------------------
2831
2832	CALL calc_month_len( iyea, imonth_len )
2833
2834	DO WHILE ( iday > imonth_len(imon) )
2835	iday = iday - imonth_len(imon)
2836	imon = imon + 1
2837	IF ( imon > 12 ) THEN
2838	imon = 1
2839	iyea = iyea + 1
2840	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
2841	ENDIF
2842	END DO
2843
2844	!!----------------------------------------------------------------------
2845	!! Convert it into YYYYMMDD.HHMMSS format.
2846	!!----------------------------------------------------------------------
2847	ddobsini = iyea * 10000_dp + imon * 100_dp + &
2848	& iday + ihou * 0.01_dp + imin * 0.0001_dp
2849
2850
2851	END SUBROUTINE ini_date
2852
2853	SUBROUTINE fin_date( ddobsfin )
2854	!!----------------------------------------------------------------------
2855	!! * ROUTINE fin_date *
2856	!!
2857	!! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
2858	!!
2859	!! ** Method : Get final data in double precision YYYYMMDD.HHMMSS format
2860	!!
2861	!! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format
2862	!!
2863	!! History :
2864	!! ! 06-03 (K. Mogensen) Original code
2865	!! ! 06-05 (K. Mogensen) Reformatted
2866	!! ! 06-10 (A. Weaver) Cleaning
2867	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
2868	!!----------------------------------------------------------------------
2869	USE phycst, ONLY : & ! Physical constants
2870	& rday
2871	! USE daymod, ONLY : & ! Time variables
2872	! & nmonth_len
2873	USE dom_oce, ONLY : & ! Ocean space and time domain variables
2874	& rdt
2875
2876	IMPLICIT NONE
2877
2878	!! * Arguments
2879	REAL(KIND=dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
2880
2881	!! * Local declarations
2882	INTEGER :: iyea ! date - (year, month, day, hour, minute)
2883	INTEGER :: imon
2884	INTEGER :: iday
2885	INTEGER :: ihou
2886	INTEGER :: imin
2887	INTEGER :: imday ! Number of days in month.
2888	REAL(KIND=wp) :: zdayfrc ! Fraction of day
2889
2890	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
2891
2892	!-----------------------------------------------------------------------
2893	! Initial date initialization (year, month, day, hour, minute)
2894	! (This assumes that the initial date is for 00z)
2895	!-----------------------------------------------------------------------
2896	iyea = ndate0 / 10000
2897	imon = ( ndate0 - iyea * 10000 ) / 100
2898	iday = ndate0 - iyea * 10000 - imon * 100
2899	ihou = 0
2900	imin = 0
2901
2902	!-----------------------------------------------------------------------
2903	! Compute number of days + number of hours + min since initial time
2904	!-----------------------------------------------------------------------
2905	iday = iday + nitend * rdt / rday
2906	zdayfrc = nitend * rdt / rday
2907	zdayfrc = zdayfrc - AINT( zdayfrc )
2908	ihou = INT( zdayfrc * 24 )
2909	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
2910
2911	!-----------------------------------------------------------------------
2912	! Convert number of days (iday) into a real date
2913	!----------------------------------------------------------------------
2914
2915	CALL calc_month_len( iyea, imonth_len )
2916
2917	DO WHILE ( iday > imonth_len(imon) )
2918	iday = iday - imonth_len(imon)
2919	imon = imon + 1
2920	IF ( imon > 12 ) THEN
2921	imon = 1
2922	iyea = iyea + 1
2923	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
2924	ENDIF
2925	END DO
2926
2927	!-----------------------------------------------------------------------
2928	! Convert it into YYYYMMDD.HHMMSS format
2929	!-----------------------------------------------------------------------
2930	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
2931	& + ihou * 0.01_dp + imin * 0.0001_dp
2932
2933	END SUBROUTINE fin_date
2934
2935	END MODULE diaobs

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

Download in other formats: