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.

obs_write.F90 in branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90 @ 7992

Last change on this file since 7992 was 7992, checked in by jwhile, 7 years ago
This version of the code seems to work correctly after some bug fixes
File size: 25.2 KB

Rev	Line
[2128]	1	MODULE obs_write
	2	!!======================================================================
	3	!! * MODULE obs_write *
	4	!! Observation diagnosticss: Write observation related diagnostics
	5	!!=====================================================================
	6
	7	!!----------------------------------------------------------------------
[7992]	8	!! obs_wri_prof : Write profile observations in feedback format
	9	!! obs_wri_surf : Write surface observations in feedback format
	10	!! obs_wri_stats : Print basic statistics on the data being written out
[2128]	11	!!----------------------------------------------------------------------
	12
	13	!! * Modules used
	14	USE par_kind, ONLY : & ! Precision variables
	15	& wp
	16	USE in_out_manager ! I/O manager
	17	USE dom_oce ! Ocean space and time domain variables
	18	USE obs_types ! Observation type integer to character translation
	19	USE julian, ONLY : & ! Julian date routines
	20	& greg2jul
	21	USE obs_utils, ONLY : & ! Observation operator utility functions
	22	& chkerr
	23	USE obs_profiles_def ! Type definitions for profiles
	24	USE obs_surf_def ! Type defintions for surface observations
	25	USE obs_fbm ! Observation feedback I/O
	26	USE obs_grid ! Grid tools
	27	USE obs_conv ! Conversion between units
	28	USE obs_const
[4990]	29	USE obs_mpp ! MPP support routines for observation diagnostics
	30	USE lib_mpp ! MPP routines
[2128]	31
	32	IMPLICIT NONE
	33
	34	!! * Routine accessibility
	35	PRIVATE
[7992]	36	PUBLIC obs_wri_prof, & ! Write profile observation files
	37	& obs_wri_surf, & ! Write surface observation files
[2128]	38	& obswriinfo
	39
	40	TYPE obswriinfo
	41	INTEGER :: inum
	42	INTEGER, POINTER, DIMENSION(:) :: ipoint
	43	CHARACTER(len=ilenname), POINTER, DIMENSION(:) :: cdname
	44	CHARACTER(len=ilenlong), POINTER, DIMENSION(:,:) :: cdlong
	45	CHARACTER(len=ilenunit), POINTER, DIMENSION(:,:) :: cdunit
	46	END TYPE obswriinfo
	47
[2287]	48	!!----------------------------------------------------------------------
	49	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	50	!! $Id$
	51	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	52	!!----------------------------------------------------------------------
	53
[2128]	54	CONTAINS
	55
[7992]	56	SUBROUTINE obs_wri_prof( profdata, padd, pext )
[2128]	57	!!-----------------------------------------------------------------------
	58	!!
[7992]	59	!! * ROUTINE obs_wri_prof *
[2128]	60	!!
[7992]	61	!! ** Purpose : Write profile feedback files
[2128]	62	!!
	63	!! ** Method : NetCDF
	64	!!
	65	!! ** Action :
	66	!!
	67	!! History :
	68	!! ! 06-04 (A. Vidard) Original
	69	!! ! 06-04 (A. Vidard) Reformatted
	70	!! ! 06-10 (A. Weaver) Cleanup
	71	!! ! 07-01 (K. Mogensen) Use profile data types
	72	!! ! 07-03 (K. Mogensen) General handling of profiles
	73	!! ! 09-01 (K. Mogensen) New feedback format
[7992]	74	!! ! 15-02 (M. Martin) Combined routine for writing profiles
[2128]	75	!!-----------------------------------------------------------------------
	76
	77	!! * Arguments
	78	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
	79	TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable
	80	TYPE(obswriinfo), OPTIONAL :: pext ! Extra info
[7992]	81
[2128]	82	!! * Local declarations
	83	TYPE(obfbdata) :: fbdata
[7992]	84	CHARACTER(LEN=40) :: clfname
	85	CHARACTER(LEN=6) :: clfiletype
[2128]	86	INTEGER :: ilevel
	87	INTEGER :: jvar
	88	INTEGER :: jo
	89	INTEGER :: jk
	90	INTEGER :: ik
	91	INTEGER :: ja
	92	INTEGER :: je
[7992]	93	INTEGER :: iadd
	94	INTEGER :: iext
[2128]	95	REAL(wp) :: zpres
	96
	97	IF ( PRESENT( padd ) ) THEN
[7992]	98	iadd = padd%inum
[2128]	99	ELSE
[7992]	100	iadd = 0
[2128]	101	ENDIF
	102
	103	IF ( PRESENT( pext ) ) THEN
[7992]	104	iext = pext%inum
[2128]	105	ELSE
[7992]	106	iext = 0
[2128]	107	ENDIF
[7992]	108
[2128]	109	CALL init_obfbdata( fbdata )
	110
	111	! Find maximum level
	112	ilevel = 0
	113	DO jvar = 1, 2
	114	ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) )
	115	END DO
	116
[7992]	117	SELECT CASE ( TRIM(profdata%cvars(1)) )
	118	CASE('POTM')
	119
	120	clfiletype='profb'
	121	CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, &
	122	& 1 + iadd, 1 + iext, .TRUE. )
	123	fbdata%cname(1) = profdata%cvars(1)
	124	fbdata%cname(2) = profdata%cvars(2)
	125	fbdata%coblong(1) = 'Potential temperature'
	126	fbdata%coblong(2) = 'Practical salinity'
	127	fbdata%cobunit(1) = 'Degrees centigrade'
	128	fbdata%cobunit(2) = 'PSU'
	129	fbdata%cextname(1) = 'TEMP'
	130	fbdata%cextlong(1) = 'Insitu temperature'
	131	fbdata%cextunit(1) = 'Degrees centigrade'
	132	fbdata%caddlong(1,1) = 'Model interpolated potential temperature'
	133	fbdata%caddlong(1,2) = 'Model interpolated practical salinity'
	134	fbdata%caddunit(1,1) = 'Degrees centigrade'
	135	fbdata%caddunit(1,2) = 'PSU'
	136	fbdata%cgrid(:) = 'T'
	137	DO je = 1, iext
	138	fbdata%cextname(1+je) = pext%cdname(je)
	139	fbdata%cextlong(1+je) = pext%cdlong(je,1)
	140	fbdata%cextunit(1+je) = pext%cdunit(je,1)
[2128]	141	END DO
[7992]	142	DO ja = 1, iadd
	143	fbdata%caddname(1+ja) = padd%cdname(ja)
	144	DO jvar = 1, 2
	145	fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar)
	146	fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar)
	147	END DO
	148	END DO
[2128]	149
[7992]	150	CASE('UVEL')
	151
	152	clfiletype='velfb'
	153	CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 1, 0, .TRUE. )
	154	fbdata%cname(1) = profdata%cvars(1)
	155	fbdata%cname(2) = profdata%cvars(2)
	156	fbdata%coblong(1) = 'Zonal velocity'
	157	fbdata%coblong(2) = 'Meridional velocity'
	158	fbdata%cobunit(1) = 'm/s'
	159	fbdata%cobunit(2) = 'm/s'
	160	DO je = 1, iext
	161	fbdata%cextname(je) = pext%cdname(je)
	162	fbdata%cextlong(je) = pext%cdlong(je,1)
	163	fbdata%cextunit(je) = pext%cdunit(je,1)
	164	END DO
	165	fbdata%caddlong(1,1) = 'Model interpolated zonal velocity'
	166	fbdata%caddlong(1,2) = 'Model interpolated meridional velocity'
	167	fbdata%caddunit(1,1) = 'm/s'
	168	fbdata%caddunit(1,2) = 'm/s'
	169	fbdata%cgrid(1) = 'U'
	170	fbdata%cgrid(2) = 'V'
	171	DO ja = 1, iadd
	172	fbdata%caddname(1+ja) = padd%cdname(ja)
	173	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	174	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	175	END DO
	176
	177	END SELECT
	178
	179	fbdata%caddname(1) = 'Hx'
	180
	181	WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
	182
[2128]	183	IF(lwp) THEN
	184	WRITE(numout,*)
[7992]	185	WRITE(numout,*)'obs_wri_prof :'
[2128]	186	WRITE(numout,*)'~~~~~~~~~~~~~'
[7992]	187	WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname)
[2128]	188	ENDIF
	189
[7992]	190	! Transform obs_prof data structure into obfb data structure
[2128]	191	fbdata%cdjuldref = '19500101000000'
	192	DO jo = 1, profdata%nprof
	193	fbdata%plam(jo) = profdata%rlam(jo)
	194	fbdata%pphi(jo) = profdata%rphi(jo)
	195	WRITE(fbdata%cdtyp(jo),'(I4)') profdata%ntyp(jo)
	196	fbdata%ivqc(jo,:) = profdata%ivqc(jo,:)
	197	fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:)
[7992]	198	IF ( profdata%nqc(jo) > 255 ) THEN
	199	fbdata%ioqc(jo) = IBSET(profdata%nqc(jo),2)
[2128]	200	fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo)
[7992]	201	fbdata%ioqcf(2,jo) = profdata%nqc(jo)
[2128]	202	ELSE
	203	fbdata%ioqc(jo) = profdata%nqc(jo)
	204	fbdata%ioqcf(:,jo) = profdata%nqcf(:,jo)
	205	ENDIF
	206	fbdata%ipqc(jo) = profdata%ipqc(jo)
	207	fbdata%ipqcf(:,jo) = profdata%ipqcf(:,jo)
	208	fbdata%itqc(jo) = profdata%itqc(jo)
	209	fbdata%itqcf(:,jo) = profdata%itqcf(:,jo)
	210	fbdata%cdwmo(jo) = profdata%cwmo(jo)
	211	fbdata%kindex(jo) = profdata%npfil(jo)
	212	DO jvar = 1, profdata%nvar
	213	IF (ln_grid_global) THEN
	214	fbdata%iobsi(jo,jvar) = profdata%mi(jo,jvar)
	215	fbdata%iobsj(jo,jvar) = profdata%mj(jo,jvar)
	216	ELSE
	217	fbdata%iobsi(jo,jvar) = mig(profdata%mi(jo,jvar))
	218	fbdata%iobsj(jo,jvar) = mjg(profdata%mj(jo,jvar))
	219	ENDIF
	220	END DO
	221	CALL greg2jul( 0, &
	222	& profdata%nmin(jo), &
	223	& profdata%nhou(jo), &
	224	& profdata%nday(jo), &
	225	& profdata%nmon(jo), &
	226	& profdata%nyea(jo), &
	227	& fbdata%ptim(jo), &
	228	& krefdate = 19500101 )
	229	! Reform the profiles arrays for output
	230	DO jvar = 1, 2
	231	DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar)
	232	ik = profdata%var(jvar)%nvlidx(jk)
	233	fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)
	234	fbdata%pob(ik,jo,jvar) = profdata%var(jvar)%vobs(jk)
	235	fbdata%pdep(ik,jo) = profdata%var(jvar)%vdep(jk)
	236	fbdata%idqc(ik,jo) = profdata%var(jvar)%idqc(jk)
	237	fbdata%idqcf(:,ik,jo) = profdata%var(jvar)%idqcf(:,jk)
[7992]	238	IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN
	239	fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2)
[2128]	240	fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk)
[7992]	241	fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111')
[2128]	242	ELSE
	243	fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk)
	244	fbdata%ivlqcf(:,ik,jo,jvar) = profdata%var(jvar)%nvqcf(:,jk)
	245	ENDIF
	246	fbdata%iobsk(ik,jo,jvar) = profdata%var(jvar)%mvk(jk)
[7992]	247	DO ja = 1, iadd
[2128]	248	fbdata%padd(ik,jo,1+ja,jvar) = &
	249	& profdata%var(jvar)%vext(jk,padd%ipoint(ja))
	250	END DO
[7992]	251	DO je = 1, iext
[2128]	252	fbdata%pext(ik,jo,1+je) = &
	253	& profdata%var(jvar)%vext(jk,pext%ipoint(je))
	254	END DO
[7992]	255	IF ( ( jvar == 1 ) .AND. &
	256	& ( TRIM(profdata%cvars(1)) == 'POTM' ) ) THEN
[2128]	257	fbdata%pext(ik,jo,1) = profdata%var(jvar)%vext(jk,1)
	258	ENDIF
	259	END DO
	260	END DO
	261	END DO
	262
[7992]	263	IF ( TRIM(profdata%cvars(1)) == 'POTM' ) THEN
	264	! Convert insitu temperature to potential temperature using the model
	265	! salinity if no potential temperature
	266	DO jo = 1, fbdata%nobs
	267	IF ( fbdata%pphi(jo) < 9999.0 ) THEN
	268	DO jk = 1, fbdata%nlev
	269	IF ( ( fbdata%pob(jk,jo,1) >= 9999.0 ) .AND. &
	270	& ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
	271	& ( fbdata%padd(jk,jo,1,2) < 9999.0 ) .AND. &
	272	& ( fbdata%pext(jk,jo,1) < 9999.0 ) ) THEN
	273	zpres = dep_to_p( REAL(fbdata%pdep(jk,jo),wp), &
	274	& REAL(fbdata%pphi(jo),wp) )
	275	fbdata%pob(jk,jo,1) = potemp( &
	276	& REAL(fbdata%padd(jk,jo,1,2), wp), &
	277	& REAL(fbdata%pext(jk,jo,1), wp), &
	278	& zpres, 0.0_wp )
	279	ENDIF
	280	END DO
	281	ENDIF
	282	END DO
	283	ENDIF
	284
[2128]	285	! Write the obfbdata structure
[7992]	286	CALL write_obfbdata( clfname, fbdata )
[4990]	287
	288	! Output some basic statistics
	289	CALL obs_wri_stats( fbdata )
	290
[2128]	291	CALL dealloc_obfbdata( fbdata )
	292
[7992]	293	END SUBROUTINE obs_wri_prof
	294
	295	SUBROUTINE obs_wri_surf( surfdata, padd, pext )
[2128]	296	!!-----------------------------------------------------------------------
	297	!!
[7992]	298	!! * ROUTINE obs_wri_surf *
[2128]	299	!!
[7992]	300	!! ** Purpose : Write surface observation files
[2128]	301	!!
	302	!! ** Method : NetCDF
	303	!!
	304	!! ** Action :
	305	!!
	306	!! ! 07-03 (K. Mogensen) Original
	307	!! ! 09-01 (K. Mogensen) New feedback format.
[7992]	308	!! ! 15-02 (M. Martin) Combined surface writing routine.
[2128]	309	!!-----------------------------------------------------------------------
	310
	311	!! * Modules used
	312	IMPLICIT NONE
	313
	314	!! * Arguments
[7992]	315	TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data
[2128]	316	TYPE(obswriinfo), OPTIONAL :: padd ! Additional info for each variable
	317	TYPE(obswriinfo), OPTIONAL :: pext ! Extra info
	318
	319	!! * Local declarations
	320	TYPE(obfbdata) :: fbdata
[7992]	321	CHARACTER(LEN=40) :: clfname ! netCDF filename
	322	CHARACTER(LEN=6) :: clfiletype
	323	CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
[2128]	324	INTEGER :: jo
	325	INTEGER :: ja
	326	INTEGER :: je
[7992]	327	INTEGER :: iadd
	328	INTEGER :: iext
[2128]	329
	330	IF ( PRESENT( padd ) ) THEN
[7992]	331	iadd = padd%inum
[2128]	332	ELSE
[7992]	333	iadd = 0
[2128]	334	ENDIF
	335
	336	IF ( PRESENT( pext ) ) THEN
[7992]	337	iext = pext%inum
[2128]	338	ELSE
[7992]	339	iext = 0
[2128]	340	ENDIF
	341
	342	CALL init_obfbdata( fbdata )
	343
[7992]	344	SELECT CASE ( TRIM(surfdata%cvars(1)) )
	345	CASE('SLA')
[2128]	346
[7992]	347	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	348	& 2 + iadd, 1 + iext, .TRUE. )
[2128]	349
[7992]	350	clfiletype = 'slafb'
	351	fbdata%cname(1) = surfdata%cvars(1)
	352	fbdata%coblong(1) = 'Sea level anomaly'
	353	fbdata%cobunit(1) = 'Metres'
	354	fbdata%cextname(1) = 'MDT'
	355	fbdata%cextlong(1) = 'Mean dynamic topography'
	356	fbdata%cextunit(1) = 'Metres'
	357	DO je = 1, iext
	358	fbdata%cextname(je) = pext%cdname(je)
	359	fbdata%cextlong(je) = pext%cdlong(je,1)
	360	fbdata%cextunit(je) = pext%cdunit(je,1)
[2128]	361	END DO
[7992]	362	fbdata%caddlong(1,1) = 'Model interpolated SSH - MDT'
	363	fbdata%caddunit(1,1) = 'Metres'
	364	fbdata%caddname(2) = 'SSH'
	365	fbdata%caddlong(2,1) = 'Model Sea surface height'
	366	fbdata%caddunit(2,1) = 'Metres'
	367	fbdata%cgrid(1) = 'T'
	368	DO ja = 1, iadd
	369	fbdata%caddname(2+ja) = padd%cdname(ja)
	370	fbdata%caddlong(2+ja,1) = padd%cdlong(ja,1)
	371	fbdata%caddunit(2+ja,1) = padd%cdunit(ja,1)
[2128]	372	END DO
	373
[7992]	374	CASE('SST')
[2128]	375
[7992]	376	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	377	& 1 + iadd, iext, .TRUE. )
[4990]	378
[7992]	379	clfiletype = 'sstfb'
	380	fbdata%cname(1) = surfdata%cvars(1)
	381	fbdata%coblong(1) = 'Sea surface temperature'
	382	fbdata%cobunit(1) = 'Degree centigrade'
	383	DO je = 1, iext
	384	fbdata%cextname(je) = pext%cdname(je)
	385	fbdata%cextlong(je) = pext%cdlong(je,1)
	386	fbdata%cextunit(je) = pext%cdunit(je,1)
	387	END DO
	388	fbdata%caddlong(1,1) = 'Model interpolated SST'
	389	fbdata%caddunit(1,1) = 'Degree centigrade'
	390	fbdata%cgrid(1) = 'T'
	391	DO ja = 1, iadd
	392	fbdata%caddname(1+ja) = padd%cdname(ja)
	393	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	394	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	395	END DO
[2128]	396
[7992]	397	CASE('ICECONC')
[2128]	398
[7992]	399	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	400	& 1 + iadd, iext, .TRUE. )
[2128]	401
[7992]	402	clfiletype = 'sicfb'
	403	fbdata%cname(1) = surfdata%cvars(1)
	404	fbdata%coblong(1) = 'Sea ice'
	405	fbdata%cobunit(1) = 'Fraction'
	406	DO je = 1, iext
	407	fbdata%cextname(je) = pext%cdname(je)
	408	fbdata%cextlong(je) = pext%cdlong(je,1)
	409	fbdata%cextunit(je) = pext%cdunit(je,1)
	410	END DO
	411	fbdata%caddlong(1,1) = 'Model interpolated ICE'
	412	fbdata%caddunit(1,1) = 'Fraction'
	413	fbdata%cgrid(1) = 'T'
	414	DO ja = 1, iadd
	415	fbdata%caddname(1+ja) = padd%cdname(ja)
	416	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	417	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	418	END DO
[2128]	419
[7992]	420	CASE('SSS')
[2128]	421
[7992]	422	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	423	& 1 + iadd, iext, .TRUE. )
[2128]	424
[7992]	425	clfiletype = 'sssfb'
	426	fbdata%cname(1) = surfdata%cvars(1)
	427	fbdata%coblong(1) = 'Sea surface salinity'
	428	fbdata%cobunit(1) = 'psu'
	429	DO je = 1, iext
	430	fbdata%cextname(je) = pext%cdname(je)
	431	fbdata%cextlong(je) = pext%cdlong(je,1)
	432	fbdata%cextunit(je) = pext%cdunit(je,1)
	433	END DO
	434	fbdata%caddlong(1,1) = 'Model interpolated SSS'
	435	fbdata%caddunit(1,1) = 'psu'
	436	fbdata%cgrid(1) = 'T'
	437	DO ja = 1, iadd
	438	fbdata%caddname(1+ja) = padd%cdname(ja)
	439	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	440	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	441	END DO
[2128]	442
[7992]	443	CASE('LOGCHL')
[2128]	444
[7992]	445	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	446	& 1 + iadd, iext, .TRUE. )
[2128]	447
[7992]	448	clfiletype = 'logchlfb'
	449	fbdata%cname(1) = surfdata%cvars(1)
	450	fbdata%coblong(1) = 'logchl concentration'
	451	fbdata%cobunit(1) = 'mg/m3'
	452	DO je = 1, iext
	453	fbdata%cextname(je) = pext%cdname(je)
	454	fbdata%cextlong(je) = pext%cdlong(je,1)
	455	fbdata%cextunit(je) = pext%cdunit(je,1)
	456	END DO
	457	fbdata%caddlong(1,1) = 'Model interpolated LOGCHL'
	458	fbdata%caddunit(1,1) = 'mg/m3'
	459	fbdata%cgrid(1) = 'T'
	460	DO ja = 1, iadd
	461	fbdata%caddname(1+ja) = padd%cdname(ja)
	462	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	463	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	464	END DO
[2128]	465
[7992]	466	CASE('SPM')
[2128]	467
[7992]	468	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	469	& 1 + iadd, iext, .TRUE. )
[2128]	470
[7992]	471	clfiletype = 'spmfb'
	472	fbdata%cname(1) = surfdata%cvars(1)
	473	fbdata%coblong(1) = 'spm'
	474	fbdata%cobunit(1) = 'g/m3'
	475	DO je = 1, iext
	476	fbdata%cextname(je) = pext%cdname(je)
	477	fbdata%cextlong(je) = pext%cdlong(je,1)
	478	fbdata%cextunit(je) = pext%cdunit(je,1)
[2128]	479	END DO
[7992]	480	fbdata%caddlong(1,1) = 'Model interpolated spm'
	481	fbdata%caddunit(1,1) = 'g/m3'
	482	fbdata%cgrid(1) = 'T'
	483	DO ja = 1, iadd
	484	fbdata%caddname(1+ja) = padd%cdname(ja)
	485	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	486	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
[2128]	487	END DO
	488
[7992]	489	CASE('FCO2')
[2128]	490
[7992]	491	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	492	& 1 + iadd, iext, .TRUE. )
[2128]	493
[7992]	494	clfiletype = 'fco2fb'
	495	fbdata%cname(1) = surfdata%cvars(1)
	496	fbdata%coblong(1) = 'fco2'
	497	fbdata%cobunit(1) = 'uatm'
	498	DO je = 1, iext
	499	fbdata%cextname(je) = pext%cdname(je)
	500	fbdata%cextlong(je) = pext%cdlong(je,1)
	501	fbdata%cextunit(je) = pext%cdunit(je,1)
	502	END DO
	503	fbdata%caddlong(1,1) = 'Model interpolated fco2'
	504	fbdata%caddunit(1,1) = 'uatm'
	505	fbdata%cgrid(1) = 'T'
	506	DO ja = 1, iadd
	507	fbdata%caddname(1+ja) = padd%cdname(ja)
	508	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	509	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	510	END DO
[2128]	511
[7992]	512	CASE('PCO2')
[4990]	513
[7992]	514	CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, &
	515	& 1 + iadd, iext, .TRUE. )
[2128]	516
[7992]	517	clfiletype = 'pco2fb'
	518	fbdata%cname(1) = surfdata%cvars(1)
	519	fbdata%coblong(1) = 'pco2'
	520	fbdata%cobunit(1) = 'uatm'
	521	DO je = 1, iext
	522	fbdata%cextname(je) = pext%cdname(je)
	523	fbdata%cextlong(je) = pext%cdlong(je,1)
	524	fbdata%cextunit(je) = pext%cdunit(je,1)
	525	END DO
	526	fbdata%caddlong(1,1) = 'Model interpolated pco2'
	527	fbdata%caddunit(1,1) = 'uatm'
	528	fbdata%cgrid(1) = 'T'
	529	DO ja = 1, iadd
	530	fbdata%caddname(1+ja) = padd%cdname(ja)
	531	fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1)
	532	fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1)
	533	END DO
[2128]	534
[7992]	535	CASE DEFAULT
[2128]	536
[7992]	537	CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' )
[2128]	538
[7992]	539	END SELECT
[2128]	540
	541	fbdata%caddname(1) = 'Hx'
	542
[7992]	543	WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
[2128]	544
	545	IF(lwp) THEN
	546	WRITE(numout,*)
[7992]	547	WRITE(numout,*)'obs_wri_surf :'
	548	WRITE(numout,*)'~~~~~~~~~~~~~'
	549	WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname)
[2128]	550	ENDIF
	551
[7992]	552	! Transform surf data structure into obfbdata structure
[2128]	553	fbdata%cdjuldref = '19500101000000'
[7992]	554	DO jo = 1, surfdata%nsurf
	555	fbdata%plam(jo) = surfdata%rlam(jo)
	556	fbdata%pphi(jo) = surfdata%rphi(jo)
	557	WRITE(fbdata%cdtyp(jo),'(I4)') surfdata%ntyp(jo)
[2128]	558	fbdata%ivqc(jo,:) = 0
	559	fbdata%ivqcf(:,jo,:) = 0
[7992]	560	IF ( surfdata%nqc(jo) > 255 ) THEN
[2128]	561	fbdata%ioqc(jo) = 4
	562	fbdata%ioqcf(1,jo) = 0
[7992]	563	fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
[2128]	564	ELSE
[7992]	565	fbdata%ioqc(jo) = surfdata%nqc(jo)
[2128]	566	fbdata%ioqcf(:,jo) = 0
	567	ENDIF
	568	fbdata%ipqc(jo) = 0
	569	fbdata%ipqcf(:,jo) = 0
	570	fbdata%itqc(jo) = 0
	571	fbdata%itqcf(:,jo) = 0
[7992]	572	fbdata%cdwmo(jo) = surfdata%cwmo(jo)
	573	fbdata%kindex(jo) = surfdata%nsfil(jo)
[2128]	574	IF (ln_grid_global) THEN
[7992]	575	fbdata%iobsi(jo,1) = surfdata%mi(jo)
	576	fbdata%iobsj(jo,1) = surfdata%mj(jo)
[2128]	577	ELSE
[7992]	578	fbdata%iobsi(jo,1) = mig(surfdata%mi(jo))
	579	fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo))
[2128]	580	ENDIF
	581	CALL greg2jul( 0, &
[7992]	582	& surfdata%nmin(jo), &
	583	& surfdata%nhou(jo), &
	584	& surfdata%nday(jo), &
	585	& surfdata%nmon(jo), &
	586	& surfdata%nyea(jo), &
[2128]	587	& fbdata%ptim(jo), &
	588	& krefdate = 19500101 )
[7992]	589	fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)
	590	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)
	591	fbdata%pob(1,jo,1) = surfdata%robs(jo,1)
[2128]	592	fbdata%pdep(1,jo) = 0.0
	593	fbdata%idqc(1,jo) = 0
	594	fbdata%idqcf(:,1,jo) = 0
[7992]	595	IF ( surfdata%nqc(jo) > 255 ) THEN
	596	fbdata%ivqc(jo,1) = 4
	597	fbdata%ivlqc(1,jo,1) = 4
[2128]	598	fbdata%ivlqcf(1,1,jo,1) = 0
[7992]	599	fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111')
[2128]	600	ELSE
[7992]	601	fbdata%ivqc(jo,1) = surfdata%nqc(jo)
	602	fbdata%ivlqc(1,jo,1) = surfdata%nqc(jo)
[2128]	603	fbdata%ivlqcf(:,1,jo,1) = 0
	604	ENDIF
	605	fbdata%iobsk(1,jo,1) = 0
[7992]	606	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)
	607	DO ja = 1, iadd
	608	fbdata%padd(1,jo,2+ja,1) = &
	609	& surfdata%rext(jo,padd%ipoint(ja))
[2128]	610	END DO
[7992]	611	DO je = 1, iext
	612	fbdata%pext(1,jo,1+je) = &
	613	& surfdata%rext(jo,pext%ipoint(je))
[2128]	614	END DO
	615	END DO
	616
	617	! Write the obfbdata structure
[7992]	618	CALL write_obfbdata( clfname, fbdata )
[2128]	619
[4990]	620	! Output some basic statistics
	621	CALL obs_wri_stats( fbdata )
	622
[2128]	623	CALL dealloc_obfbdata( fbdata )
	624
[7992]	625	END SUBROUTINE obs_wri_surf
[2128]	626
[4990]	627	SUBROUTINE obs_wri_stats( fbdata )
	628	!!-----------------------------------------------------------------------
	629	!!
	630	!! * ROUTINE obs_wri_stats *
	631	!!
	632	!! ** Purpose : Output some basic statistics of the data being written out
	633	!!
	634	!! ** Method :
	635	!!
	636	!! ** Action :
	637	!!
	638	!! ! 2014-08 (D. Lea) Initial version
	639	!!-----------------------------------------------------------------------
	640
	641	!! * Arguments
	642	TYPE(obfbdata) :: fbdata
	643
	644	!! * Local declarations
	645	INTEGER :: jvar
	646	INTEGER :: jo
	647	INTEGER :: jk
[7992]	648	INTEGER :: inumgoodobs
	649	INTEGER :: inumgoodobsmpp
[4990]	650	REAL(wp) :: zsumx
	651	REAL(wp) :: zsumx2
	652	REAL(wp) :: zomb
[7992]	653
[4990]	654
	655	IF (lwp) THEN
	656	WRITE(numout,*) ''
	657	WRITE(numout,*) 'obs_wri_stats :'
[7992]	658	WRITE(numout,*) '~~~~~~~~~~~~~~~'
[4990]	659	ENDIF
	660
	661	DO jvar = 1, fbdata%nvar
	662	zsumx=0.0_wp
	663	zsumx2=0.0_wp
[7992]	664	inumgoodobs=0
[4990]	665	DO jo = 1, fbdata%nobs
	666	DO jk = 1, fbdata%nlev
	667	IF ( ( fbdata%pob(jk,jo,jvar) < 9999.0 ) .AND. &
	668	& ( fbdata%pdep(jk,jo) < 9999.0 ) .AND. &
	669	& ( fbdata%padd(jk,jo,1,jvar) < 9999.0 ) ) THEN
[7992]	670
	671	zomb=fbdata%pob(jk, jo, jvar)-fbdata%padd(jk, jo, 1, jvar)
[4990]	672	zsumx=zsumx+zomb
	673	zsumx2=zsumx2+zomb**2
[7992]	674	inumgoodobs=inumgoodobs+1
	675	ENDIF
[4990]	676	ENDDO
	677	ENDDO
	678
[7992]	679	CALL obs_mpp_sum_integer( inumgoodobs, inumgoodobsmpp )
[4990]	680	CALL mpp_sum(zsumx)
	681	CALL mpp_sum(zsumx2)
	682
	683	IF (lwp) THEN
[7992]	684	WRITE(numout,*) 'Type: ',fbdata%cname(jvar),' Total number of good observations: ',inumgoodobsmpp
	685	WRITE(numout,*) 'Overall mean obs minus model of the good observations: ',zsumx/inumgoodobsmpp
	686	WRITE(numout,*) 'Overall RMS obs minus model of the good observations: ',sqrt( zsumx2/inumgoodobsmpp )
	687	WRITE(numout,*) ''
[4990]	688	ENDIF
[7992]	689
[4990]	690	ENDDO
	691
	692	END SUBROUTINE obs_wri_stats
	693
[2128]	694	END MODULE obs_write

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

Download in other formats: