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_oper.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 48.8 KB

Rev	Line
[2128]	1	MODULE obs_oper
	2	!!======================================================================
	3	!! * MODULE obs_oper *
	4	!! Observation diagnostics: Observation operators for various observation
	5	!! types
	6	!!======================================================================
	7
	8	!!----------------------------------------------------------------------
	9	!! obs_pro_opt : Compute the model counterpart of temperature and
	10	!! salinity observations from profiles
	11	!! obs_sla_opt : Compute the model counterpart of sea level anomaly
	12	!! observations
	13	!! obs_sst_opt : Compute the model counterpart of sea surface temperature
	14	!! observations
	15	!! obs_sss_opt : Compute the model counterpart of sea surface salinity
	16	!! observations
	17	!! obs_seaice_opt : Compute the model counterpart of sea ice concentration
	18	!! observations
	19	!!
	20	!! obs_vel_opt : Compute the model counterpart of zonal and meridional
	21	!! components of velocity from observations.
	22	!!----------------------------------------------------------------------
	23
	24	!! * Modules used
	25	USE par_kind, ONLY : & ! Precision variables
	26	& wp
	27	USE in_out_manager ! I/O manager
	28	USE obs_inter_sup ! Interpolation support
	29	USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the observation pt
	30	& obs_int_h2d, &
	31	& obs_int_h2d_init
	32	USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the observation pt
	33	& obs_int_z1d, &
	34	& obs_int_z1d_spl
	35	USE obs_const, ONLY : &
	36	& obfillflt ! Fillvalue
	37	USE dom_oce, ONLY : &
	38	& glamt, glamu, glamv, &
	39	& gphit, gphiu, gphiv
[2715]	40	USE lib_mpp, ONLY : &
	41	& ctl_warn, ctl_stop
[2128]	42
	43	IMPLICIT NONE
	44
	45	!! * Routine accessibility
	46	PRIVATE
	47
	48	PUBLIC obs_pro_opt, & ! Compute the model counterpart of profile observations
	49	& obs_sla_opt, & ! Compute the model counterpart of SLA observations
	50	& obs_sst_opt, & ! Compute the model counterpart of SST observations
	51	& obs_sss_opt, & ! Compute the model counterpart of SSS observations
	52	& obs_seaice_opt, &
	53	& obs_vel_opt ! Compute the model counterpart of velocity profile data
	54
	55	INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
	56
[3211]	57	!! * Control permutation of array indices
	58	!! None required from dom_oce_ftrans.h90
	59
[2287]	60	!!----------------------------------------------------------------------
	61	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	62	!! $Id$
	63	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	64	!!----------------------------------------------------------------------
	65
[2128]	66	CONTAINS
	67
	68	SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
	69	& ptn, psn, pgdept, ptmask, k1dint, k2dint, &
	70	& kdailyavtypes )
	71	!!-----------------------------------------------------------------------
	72	!!
	73	!! * ROUTINE obs_pro_opt *
	74	!!
	75	!! ** Purpose : Compute the model counterpart of profiles
	76	!! data by interpolating from the model grid to the
	77	!! observation point.
	78	!!
	79	!! ** Method : Linearly interpolate to each observation point using
	80	!! the model values at the corners of the surrounding grid box.
	81	!!
	82	!! First, a vertical profile of horizontally interpolated model
	83	!! now temperatures is computed at the obs (lon, lat) point.
	84	!! Several horizontal interpolation schemes are available:
	85	!! - distance-weighted (great circle) (k2dint = 0)
	86	!! - distance-weighted (small angle) (k2dint = 1)
	87	!! - bilinear (geographical grid) (k2dint = 2)
	88	!! - bilinear (quadrilateral grid) (k2dint = 3)
	89	!! - polynomial (quadrilateral grid) (k2dint = 4)
	90	!!
	91	!! Next, the vertical temperature profile is interpolated to the
	92	!! data depth points. Two vertical interpolation schemes are
	93	!! available:
	94	!! - linear (k1dint = 0)
	95	!! - Cubic spline (k1dint = 1)
	96	!!
	97	!! For the cubic spline the 2nd derivative of the interpolating
	98	!! polynomial is computed before entering the vertical interpolation
	99	!! routine.
	100	!!
	101	!! For ENACT moored buoy data (e.g., TAO), the model equivalent is
	102	!! a daily mean model temperature field. So, we first compute
	103	!! the mean, then interpolate only at the end of the day.
	104	!!
	105	!! Note: the in situ temperature observations must be converted
	106	!! to potential temperature (the model variable) prior to
	107	!! assimilation.
	108	!!??????????????????????????????????????????????????????????????
	109	!! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
	110	!!??????????????????????????????????????????????????????????????
	111	!!
	112	!! ** Action :
	113	!!
	114	!! History :
	115	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
	116	!! ! 06-03 (G. Smith) NEMOVAR migration
	117	!! ! 06-10 (A. Weaver) Cleanup
	118	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
	119	!! ! 07-03 (K. Mogensen) General handling of profiles
	120	!!-----------------------------------------------------------------------
	121
	122	!! * Modules used
	123	USE obs_profiles_def ! Definition of storage space for profile obs.
	124
	125	IMPLICIT NONE
	126
	127	!! * Arguments
	128	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
	129	INTEGER, INTENT(IN) :: kt ! Time step
	130	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	131	INTEGER, INTENT(IN) :: kpj
	132	INTEGER, INTENT(IN) :: kpk
	133	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	134	! (kit000-1 = restart time)
	135	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
	136	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
	137	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
[3211]	138	!! DCSE_NEMO : this style defeats ftrans
	139	! REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
	140	! & ptn, & ! Model temperature field
	141	! & psn, & ! Model salinity field
	142	! & ptmask ! Land-sea mask
	143
	144	!FTRANS ptn psn ptmask :I :I :z
	145	REAL(KIND=wp), INTENT(IN) :: ptn(kpi,kpj,kpk) ! Model temperature field
	146	REAL(KIND=wp), INTENT(IN) :: psn(kpi,kpj,kpk) ! Model salinity field
	147	REAL(KIND=wp), INTENT(IN) :: ptmask(kpi,kpj,kpk) ! Land-sea mask
	148
[2128]	149	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
	150	& pgdept ! Model array of depth levels
	151	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
	152	& kdailyavtypes! Types for daily averages
	153	!! * Local declarations
	154	INTEGER :: ji
	155	INTEGER :: jj
	156	INTEGER :: jk
	157	INTEGER :: jobs
	158	INTEGER :: inrc
	159	INTEGER :: ipro
	160	INTEGER :: idayend
	161	INTEGER :: ista
	162	INTEGER :: iend
	163	INTEGER :: iobs
	164	INTEGER, DIMENSION(imaxavtypes) :: &
	165	& idailyavtypes
	166	REAL(KIND=wp) :: zlam
	167	REAL(KIND=wp) :: zphi
	168	REAL(KIND=wp) :: zdaystp
	169	REAL(KIND=wp), DIMENSION(kpk) :: &
	170	& zobsmask, &
	171	& zobsk, &
	172	& zobs2k
	173	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
	174	& zweig
	175	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
	176	& zmask, &
	177	& zintt, &
	178	& zints, &
	179	& zinmt, &
	180	& zinms
	181	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	182	& zglam, &
	183	& zgphi
	184	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
	185	& igrdi, &
	186	& igrdj
	187
	188	!------------------------------------------------------------------------
	189	! Local initialization
	190	!------------------------------------------------------------------------
	191	! ... Record and data counters
	192	inrc = kt - kit000 + 2
	193	ipro = prodatqc%npstp(inrc)
	194
	195	! Daily average types
	196	IF ( PRESENT(kdailyavtypes) ) THEN
	197	idailyavtypes(:) = kdailyavtypes(:)
	198	ELSE
	199	idailyavtypes(:) = -1
	200	ENDIF
	201
	202	! Initialize daily mean for first timestep
	203	idayend = MOD( kt - kit000 + 1, kdaystp )
	204
	205	! Added kt == 0 test to catch restart case
	206	IF ( idayend == 1 .OR. kt == 0) THEN
	207	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
	208	DO jk = 1, jpk
	209	DO jj = 1, jpj
	210	DO ji = 1, jpi
	211	prodatqc%vdmean(ji,jj,jk,1) = 0.0
	212	prodatqc%vdmean(ji,jj,jk,2) = 0.0
	213	END DO
	214	END DO
	215	END DO
	216	ENDIF
	217
	218	DO jk = 1, jpk
	219	DO jj = 1, jpj
	220	DO ji = 1, jpi
	221	! Increment the temperature field for computing daily mean
	222	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
	223	& + ptn(ji,jj,jk)
	224	! Increment the salinity field for computing daily mean
	225	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
	226	& + psn(ji,jj,jk)
	227	END DO
	228	END DO
	229	END DO
	230
	231	! Compute the daily mean at the end of day
	232	zdaystp = 1.0 / REAL( kdaystp )
	233	IF ( idayend == 0 ) THEN
	234	DO jk = 1, jpk
	235	DO jj = 1, jpj
	236	DO ji = 1, jpi
	237	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
	238	& * zdaystp
	239	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
	240	& * zdaystp
	241	END DO
	242	END DO
	243	END DO
	244	ENDIF
	245
	246	! Get the data for interpolation
	247	ALLOCATE( &
	248	& igrdi(2,2,ipro), &
	249	& igrdj(2,2,ipro), &
	250	& zglam(2,2,ipro), &
	251	& zgphi(2,2,ipro), &
	252	& zmask(2,2,kpk,ipro), &
	253	& zintt(2,2,kpk,ipro), &
	254	& zints(2,2,kpk,ipro) &
	255	& )
	256
	257	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	258	iobs = jobs - prodatqc%nprofup
	259	igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
	260	igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
	261	igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
	262	igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
	263	igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
	264	igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
	265	igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
	266	igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
	267	END DO
	268
	269	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
	270	CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
	271	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
	272	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn, zintt )
	273	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn, zints )
	274
	275	! At the end of the day also get interpolated means
	276	IF ( idayend == 0 ) THEN
	277
	278	ALLOCATE( &
	279	& zinmt(2,2,kpk,ipro), &
	280	& zinms(2,2,kpk,ipro) &
	281	& )
	282
	283	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
	284	& prodatqc%vdmean(:,:,:,1), zinmt )
	285	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
	286	& prodatqc%vdmean(:,:,:,2), zinms )
	287
	288	ENDIF
	289
	290	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	291
	292	iobs = jobs - prodatqc%nprofup
	293
	294	IF ( kt /= prodatqc%mstp(jobs) ) THEN
	295
	296	IF(lwp) THEN
	297	WRITE(numout,*)
	298	WRITE(numout,*) ' E R R O R : Observation', &
	299	& ' time step is not consistent with the', &
	300	& ' model time step'
	301	WRITE(numout,*) ' ========='
	302	WRITE(numout,*)
	303	WRITE(numout,*) ' Record = ', jobs, &
	304	& ' kt = ', kt, &
	305	& ' mstp = ', prodatqc%mstp(jobs), &
	306	& ' ntyp = ', prodatqc%ntyp(jobs)
	307	ENDIF
	308	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
	309	ENDIF
	310
	311	zlam = prodatqc%rlam(jobs)
	312	zphi = prodatqc%rphi(jobs)
	313
	314	! Horizontal weights and vertical mask
	315
	316	IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
	317	& ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
	318
	319	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
	320	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	321	& zmask(:,:,:,iobs), zweig, zobsmask )
	322
	323	ENDIF
	324
	325	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
	326
	327	zobsk(:) = obfillflt
	328
	329	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
	330
	331	IF ( idayend == 0 ) THEN
	332
	333	! Daily averaged moored buoy (MRB) data
	334
	335	CALL obs_int_h2d( kpk, kpk, &
	336	& zweig, zinmt(:,:,:,iobs), zobsk )
	337
	338
	339	ELSE
	340
	341	CALL ctl_stop( ' A nonzero' // &
	342	& ' number of profile T BUOY data should' // &
	343	& ' only occur at the end of a given day' )
	344
	345	ENDIF
	346
	347	ELSE
	348
	349	! Point data
	350
	351	CALL obs_int_h2d( kpk, kpk, &
	352	& zweig, zintt(:,:,:,iobs), zobsk )
	353
	354	ENDIF
	355
	356	!-------------------------------------------------------------
	357	! Compute vertical second-derivative of the interpolating
	358	! polynomial at obs points
	359	!-------------------------------------------------------------
	360
	361	IF ( k1dint == 1 ) THEN
	362	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
	363	& pgdept, zobsmask )
	364	ENDIF
	365
	366	!-----------------------------------------------------------------
	367	! Vertical interpolation to the observation point
	368	!-----------------------------------------------------------------
	369	ista = prodatqc%npvsta(jobs,1)
	370	iend = prodatqc%npvend(jobs,1)
	371	CALL obs_int_z1d( kpk, &
	372	& prodatqc%var(1)%mvk(ista:iend), &
	373	& k1dint, iend - ista + 1, &
	374	& prodatqc%var(1)%vdep(ista:iend), &
	375	& zobsk, zobs2k, &
	376	& prodatqc%var(1)%vmod(ista:iend), &
	377	& pgdept, zobsmask )
	378
	379	ENDIF
	380
	381	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
	382
	383	zobsk(:) = obfillflt
	384
	385	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
	386
	387	IF ( idayend == 0 ) THEN
	388
	389	! Daily averaged moored buoy (MRB) data
	390
	391	CALL obs_int_h2d( kpk, kpk, &
	392	& zweig, zinms(:,:,:,iobs), zobsk )
	393
	394	ELSE
	395
	396	CALL ctl_stop( ' A nonzero' // &
	397	& ' number of profile S BUOY data should' // &
	398	& ' only occur at the end of a given day' )
	399
	400	ENDIF
	401
	402	ELSE
	403
	404	! Point data
	405
	406	CALL obs_int_h2d( kpk, kpk, &
	407	& zweig, zints(:,:,:,iobs), zobsk )
	408
	409	ENDIF
	410
	411
	412	!-------------------------------------------------------------
	413	! Compute vertical second-derivative of the interpolating
	414	! polynomial at obs points
	415	!-------------------------------------------------------------
	416
	417	IF ( k1dint == 1 ) THEN
	418	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
	419	& pgdept, zobsmask )
	420	ENDIF
	421
	422	!----------------------------------------------------------------
	423	! Vertical interpolation to the observation point
	424	!----------------------------------------------------------------
	425	ista = prodatqc%npvsta(jobs,2)
	426	iend = prodatqc%npvend(jobs,2)
	427	CALL obs_int_z1d( kpk, &
	428	& prodatqc%var(2)%mvk(ista:iend),&
	429	& k1dint, iend - ista + 1, &
	430	& prodatqc%var(2)%vdep(ista:iend),&
	431	& zobsk, zobs2k, &
	432	& prodatqc%var(2)%vmod(ista:iend),&
	433	& pgdept, zobsmask )
	434
	435	ENDIF
	436
	437	END DO
	438
	439	! Deallocate the data for interpolation
	440	DEALLOCATE( &
	441	& igrdi, &
	442	& igrdj, &
	443	& zglam, &
	444	& zgphi, &
	445	& zmask, &
	446	& zintt, &
	447	& zints &
	448	& )
	449	! At the end of the day also get interpolated means
	450	IF ( idayend == 0 ) THEN
	451	DEALLOCATE( &
	452	& zinmt, &
	453	& zinms &
	454	& )
	455	ENDIF
	456
	457	prodatqc%nprofup = prodatqc%nprofup + ipro
	458
	459	END SUBROUTINE obs_pro_opt
	460
	461	SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
	462	& psshn, psshmask, k2dint )
	463	!!-----------------------------------------------------------------------
	464	!!
	465	!! * ROUTINE obs_sla_opt *
	466	!!
	467	!! ** Purpose : Compute the model counterpart of sea level anomaly
	468	!! data by interpolating from the model grid to the
	469	!! observation point.
	470	!!
	471	!! ** Method : Linearly interpolate to each observation point using
	472	!! the model values at the corners of the surrounding grid box.
	473	!!
	474	!! The now model SSH is first computed at the obs (lon, lat) point.
	475	!!
	476	!! Several horizontal interpolation schemes are available:
	477	!! - distance-weighted (great circle) (k2dint = 0)
	478	!! - distance-weighted (small angle) (k2dint = 1)
	479	!! - bilinear (geographical grid) (k2dint = 2)
	480	!! - bilinear (quadrilateral grid) (k2dint = 3)
	481	!! - polynomial (quadrilateral grid) (k2dint = 4)
	482	!!
	483	!! The sea level anomaly at the observation points is then computed
	484	!! by removing a mean dynamic topography (defined at the obs. point).
	485	!!
	486	!! ** Action :
	487	!!
	488	!! History :
	489	!! ! 07-03 (A. Weaver)
	490	!!-----------------------------------------------------------------------
	491
	492	!! * Modules used
	493	USE obs_surf_def ! Definition of storage space for surface observations
	494
	495	IMPLICIT NONE
	496
	497	!! * Arguments
	498	TYPE(obs_surf), INTENT(INOUT) :: sladatqc ! Subset of surface data not failing screening
	499	INTEGER, INTENT(IN) :: kt ! Time step
	500	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	501	INTEGER, INTENT(IN) :: kpj
	502	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	503	! (kit000-1 = restart time)
	504	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
	505	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
	506	& psshn, & ! Model SSH field
	507	& psshmask ! Land-sea mask
	508
	509	!! * Local declarations
	510	INTEGER :: ji
	511	INTEGER :: jj
	512	INTEGER :: jobs
	513	INTEGER :: inrc
	514	INTEGER :: isla
	515	INTEGER :: iobs
	516	REAL(KIND=wp) :: zlam
	517	REAL(KIND=wp) :: zphi
	518	REAL(KIND=wp) :: zext(1), zobsmask(1)
	519	REAL(kind=wp), DIMENSION(2,2,1) :: &
	520	& zweig
	521	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	522	& zmask, &
	523	& zsshl, &
	524	& zglam, &
	525	& zgphi
	526	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
	527	& igrdi, &
	528	& igrdj
	529
	530	!------------------------------------------------------------------------
	531	! Local initialization
	532	!------------------------------------------------------------------------
	533	! ... Record and data counters
	534	inrc = kt - kit000 + 2
	535	isla = sladatqc%nsstp(inrc)
	536
	537	! Get the data for interpolation
	538
	539	ALLOCATE( &
	540	& igrdi(2,2,isla), &
	541	& igrdj(2,2,isla), &
	542	& zglam(2,2,isla), &
	543	& zgphi(2,2,isla), &
	544	& zmask(2,2,isla), &
	545	& zsshl(2,2,isla) &
	546	& )
	547
	548	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
	549	iobs = jobs - sladatqc%nsurfup
	550	igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
	551	igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
	552	igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
	553	igrdj(1,2,iobs) = sladatqc%mj(jobs)
	554	igrdi(2,1,iobs) = sladatqc%mi(jobs)
	555	igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
	556	igrdi(2,2,iobs) = sladatqc%mi(jobs)
	557	igrdj(2,2,iobs) = sladatqc%mj(jobs)
	558	END DO
	559
	560	CALL obs_int_comm_2d( 2, 2, isla, &
	561	& igrdi, igrdj, glamt, zglam )
	562	CALL obs_int_comm_2d( 2, 2, isla, &
	563	& igrdi, igrdj, gphit, zgphi )
	564	CALL obs_int_comm_2d( 2, 2, isla, &
	565	& igrdi, igrdj, psshmask, zmask )
	566	CALL obs_int_comm_2d( 2, 2, isla, &
	567	& igrdi, igrdj, psshn, zsshl )
	568
	569	! Loop over observations
	570
	571	DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
	572
	573	iobs = jobs - sladatqc%nsurfup
	574
	575	IF ( kt /= sladatqc%mstp(jobs) ) THEN
	576
	577	IF(lwp) THEN
	578	WRITE(numout,*)
	579	WRITE(numout,*) ' E R R O R : Observation', &
	580	& ' time step is not consistent with the', &
	581	& ' model time step'
	582	WRITE(numout,*) ' ========='
	583	WRITE(numout,*)
	584	WRITE(numout,*) ' Record = ', jobs, &
	585	& ' kt = ', kt, &
	586	& ' mstp = ', sladatqc%mstp(jobs), &
	587	& ' ntyp = ', sladatqc%ntyp(jobs)
	588	ENDIF
	589	CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
	590
	591	ENDIF
	592
	593	zlam = sladatqc%rlam(jobs)
	594	zphi = sladatqc%rphi(jobs)
	595
	596	! Get weights to interpolate the model SSH to the observation point
	597	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
	598	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	599	& zmask(:,:,iobs), zweig, zobsmask )
	600
	601
	602	! Interpolate the model SSH to the observation point
	603	CALL obs_int_h2d( 1, 1, &
	604	& zweig, zsshl(:,:,iobs), zext )
	605
	606	sladatqc%rext(jobs,1) = zext(1)
	607	! ... Remove the MDT at the observation point
	608	sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
	609
	610	END DO
	611
	612	! Deallocate the data for interpolation
	613	DEALLOCATE( &
	614	& igrdi, &
	615	& igrdj, &
	616	& zglam, &
	617	& zgphi, &
	618	& zmask, &
	619	& zsshl &
	620	& )
	621
	622	sladatqc%nsurfup = sladatqc%nsurfup + isla
	623
	624	END SUBROUTINE obs_sla_opt
	625
	626	SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, &
	627	& psstn, psstmask, k2dint )
	628
	629	!!-----------------------------------------------------------------------
	630	!!
	631	!! * ROUTINE obs_sst_opt *
	632	!!
	633	!! ** Purpose : Compute the model counterpart of surface temperature
	634	!! data by interpolating from the model grid to the
	635	!! observation point.
	636	!!
	637	!! ** Method : Linearly interpolate to each observation point using
	638	!! the model values at the corners of the surrounding grid box.
	639	!!
	640	!! The now model SST is first computed at the obs (lon, lat) point.
	641	!!
	642	!! Several horizontal interpolation schemes are available:
	643	!! - distance-weighted (great circle) (k2dint = 0)
	644	!! - distance-weighted (small angle) (k2dint = 1)
	645	!! - bilinear (geographical grid) (k2dint = 2)
	646	!! - bilinear (quadrilateral grid) (k2dint = 3)
	647	!! - polynomial (quadrilateral grid) (k2dint = 4)
	648	!!
	649	!!
	650	!! ** Action :
	651	!!
	652	!! History :
	653	!! ! 07-07 (S. Ricci ) : Original
	654	!!
	655	!!-----------------------------------------------------------------------
	656
	657	!! * Modules used
	658	USE obs_surf_def ! Definition of storage space for surface observations
	659
	660	IMPLICIT NONE
	661
	662	!! * Arguments
	663	TYPE(obs_surf), INTENT(INOUT) :: &
	664	& sstdatqc ! Subset of surface data not failing screening
	665	INTEGER, INTENT(IN) :: kt ! Time step
	666	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	667	INTEGER, INTENT(IN) :: kpj
	668	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	669	! (kit000-1 = restart time)
	670	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
	671	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
	672	& psstn, & ! Model SST field
	673	& psstmask ! Land-sea mask
	674
	675	!! * Local declarations
	676	INTEGER :: ji
	677	INTEGER :: jj
	678	INTEGER :: jobs
	679	INTEGER :: inrc
	680	INTEGER :: isst
	681	INTEGER :: iobs
	682	REAL(KIND=wp) :: zlam
	683	REAL(KIND=wp) :: zphi
	684	REAL(KIND=wp) :: zext(1), zobsmask(1)
	685	REAL(kind=wp), DIMENSION(2,2,1) :: &
	686	& zweig
	687	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	688	& zmask, &
	689	& zsstl, &
	690	& zglam, &
	691	& zgphi
	692	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
	693	& igrdi, &
	694	& igrdj
	695
	696	!-----------------------------------------------------------------------
	697	! Local initialization
	698	!-----------------------------------------------------------------------
	699	! ... Record and data counters
	700	inrc = kt - kit000 + 2
	701	isst = sstdatqc%nsstp(inrc)
	702
	703	! Get the data for interpolation
	704
	705	ALLOCATE( &
	706	& igrdi(2,2,isst), &
	707	& igrdj(2,2,isst), &
	708	& zglam(2,2,isst), &
	709	& zgphi(2,2,isst), &
	710	& zmask(2,2,isst), &
	711	& zsstl(2,2,isst) &
	712	& )
	713
	714	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
	715	iobs = jobs - sstdatqc%nsurfup
	716	igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
	717	igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
	718	igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
	719	igrdj(1,2,iobs) = sstdatqc%mj(jobs)
	720	igrdi(2,1,iobs) = sstdatqc%mi(jobs)
	721	igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
	722	igrdi(2,2,iobs) = sstdatqc%mi(jobs)
	723	igrdj(2,2,iobs) = sstdatqc%mj(jobs)
	724	END DO
	725
	726	CALL obs_int_comm_2d( 2, 2, isst, &
	727	& igrdi, igrdj, glamt, zglam )
	728	CALL obs_int_comm_2d( 2, 2, isst, &
	729	& igrdi, igrdj, gphit, zgphi )
	730	CALL obs_int_comm_2d( 2, 2, isst, &
	731	& igrdi, igrdj, psstmask, zmask )
	732	CALL obs_int_comm_2d( 2, 2, isst, &
	733	& igrdi, igrdj, psstn, zsstl )
	734
	735	! Loop over observations
	736
	737	DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
	738
	739	iobs = jobs - sstdatqc%nsurfup
	740
	741	IF ( kt /= sstdatqc%mstp(jobs) ) THEN
	742
	743	IF(lwp) THEN
	744	WRITE(numout,*)
	745	WRITE(numout,*) ' E R R O R : Observation', &
	746	& ' time step is not consistent with the', &
	747	& ' model time step'
	748	WRITE(numout,*) ' ========='
	749	WRITE(numout,*)
	750	WRITE(numout,*) ' Record = ', jobs, &
	751	& ' kt = ', kt, &
	752	& ' mstp = ', sstdatqc%mstp(jobs), &
	753	& ' ntyp = ', sstdatqc%ntyp(jobs)
	754	ENDIF
	755	CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
	756
	757	ENDIF
	758
	759	zlam = sstdatqc%rlam(jobs)
	760	zphi = sstdatqc%rphi(jobs)
	761
	762	! Get weights to interpolate the model SST to the observation point
	763	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
	764	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	765	& zmask(:,:,iobs), zweig, zobsmask )
	766
	767	! Interpolate the model SST to the observation point
	768	CALL obs_int_h2d( 1, 1, &
	769	& zweig, zsstl(:,:,iobs), zext )
	770
	771	sstdatqc%rmod(jobs,1) = zext(1)
	772
	773	END DO
	774
	775	! Deallocate the data for interpolation
	776	DEALLOCATE( &
	777	& igrdi, &
	778	& igrdj, &
	779	& zglam, &
	780	& zgphi, &
	781	& zmask, &
	782	& zsstl &
	783	& )
	784
	785	sstdatqc%nsurfup = sstdatqc%nsurfup + isst
	786
	787	END SUBROUTINE obs_sst_opt
	788
	789	SUBROUTINE obs_sss_opt
	790	!!-----------------------------------------------------------------------
	791	!!
	792	!! * ROUTINE obs_sss_opt *
	793	!!
	794	!! ** Purpose : Compute the model counterpart of sea surface salinity
	795	!! data by interpolating from the model grid to the
	796	!! observation point.
	797	!!
	798	!! ** Method :
	799	!!
	800	!! ** Action :
	801	!!
	802	!! History :
	803	!! ! ??-??
	804	!!-----------------------------------------------------------------------
	805
	806	IMPLICIT NONE
	807
	808	END SUBROUTINE obs_sss_opt
	809
	810	SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
	811	& pseaicen, pseaicemask, k2dint )
	812
	813	!!-----------------------------------------------------------------------
	814	!!
	815	!! * ROUTINE obs_seaice_opt *
	816	!!
	817	!! ** Purpose : Compute the model counterpart of surface temperature
	818	!! data by interpolating from the model grid to the
	819	!! observation point.
	820	!!
	821	!! ** Method : Linearly interpolate to each observation point using
	822	!! the model values at the corners of the surrounding grid box.
	823	!!
	824	!! The now model sea ice is first computed at the obs (lon, lat) point.
	825	!!
	826	!! Several horizontal interpolation schemes are available:
	827	!! - distance-weighted (great circle) (k2dint = 0)
	828	!! - distance-weighted (small angle) (k2dint = 1)
	829	!! - bilinear (geographical grid) (k2dint = 2)
	830	!! - bilinear (quadrilateral grid) (k2dint = 3)
	831	!! - polynomial (quadrilateral grid) (k2dint = 4)
	832	!!
	833	!!
	834	!! ** Action :
	835	!!
	836	!! History :
	837	!! ! 07-07 (S. Ricci ) : Original
	838	!!
	839	!!-----------------------------------------------------------------------
	840
	841	!! * Modules used
	842	USE obs_surf_def ! Definition of storage space for surface observations
	843
	844	IMPLICIT NONE
	845
	846	!! * Arguments
	847	TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc ! Subset of surface data not failing screening
	848	INTEGER, INTENT(IN) :: kt ! Time step
	849	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	850	INTEGER, INTENT(IN) :: kpj
	851	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	852	! (kit000-1 = restart time)
	853	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
	854	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
	855	& pseaicen, & ! Model sea ice field
	856	& pseaicemask ! Land-sea mask
	857
	858	!! * Local declarations
	859	INTEGER :: ji
	860	INTEGER :: jj
	861	INTEGER :: jobs
	862	INTEGER :: inrc
	863	INTEGER :: iseaice
	864	INTEGER :: iobs
	865
	866	REAL(KIND=wp) :: zlam
	867	REAL(KIND=wp) :: zphi
	868	REAL(KIND=wp) :: zext(1), zobsmask(1)
	869	REAL(kind=wp), DIMENSION(2,2,1) :: &
	870	& zweig
	871	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	872	& zmask, &
	873	& zseaicel, &
	874	& zglam, &
	875	& zgphi
	876	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
	877	& igrdi, &
	878	& igrdj
	879
	880	!------------------------------------------------------------------------
	881	! Local initialization
	882	!------------------------------------------------------------------------
	883	! ... Record and data counters
	884	inrc = kt - kit000 + 2
	885	iseaice = seaicedatqc%nsstp(inrc)
	886
	887	! Get the data for interpolation
	888
	889	ALLOCATE( &
	890	& igrdi(2,2,iseaice), &
	891	& igrdj(2,2,iseaice), &
	892	& zglam(2,2,iseaice), &
	893	& zgphi(2,2,iseaice), &
	894	& zmask(2,2,iseaice), &
	895	& zseaicel(2,2,iseaice) &
	896	& )
	897
	898	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
	899	iobs = jobs - seaicedatqc%nsurfup
	900	igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
	901	igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
	902	igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
	903	igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
	904	igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
	905	igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
	906	igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
	907	igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
	908	END DO
	909
	910	CALL obs_int_comm_2d( 2, 2, iseaice, &
	911	& igrdi, igrdj, glamt, zglam )
	912	CALL obs_int_comm_2d( 2, 2, iseaice, &
	913	& igrdi, igrdj, gphit, zgphi )
	914	CALL obs_int_comm_2d( 2, 2, iseaice, &
	915	& igrdi, igrdj, pseaicemask, zmask )
	916	CALL obs_int_comm_2d( 2, 2, iseaice, &
	917	& igrdi, igrdj, pseaicen, zseaicel )
	918
	919	DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
	920
	921	iobs = jobs - seaicedatqc%nsurfup
	922
	923	IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
	924
	925	IF(lwp) THEN
	926	WRITE(numout,*)
	927	WRITE(numout,*) ' E R R O R : Observation', &
	928	& ' time step is not consistent with the', &
	929	& ' model time step'
	930	WRITE(numout,*) ' ========='
	931	WRITE(numout,*)
	932	WRITE(numout,*) ' Record = ', jobs, &
	933	& ' kt = ', kt, &
	934	& ' mstp = ', seaicedatqc%mstp(jobs), &
	935	& ' ntyp = ', seaicedatqc%ntyp(jobs)
	936	ENDIF
	937	CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
	938
	939	ENDIF
	940
	941	zlam = seaicedatqc%rlam(jobs)
	942	zphi = seaicedatqc%rphi(jobs)
	943
	944	! Get weights to interpolate the model sea ice to the observation point
	945	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
	946	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	947	& zmask(:,:,iobs), zweig, zobsmask )
	948
	949	! ... Interpolate the model sea ice to the observation point
	950	CALL obs_int_h2d( 1, 1, &
	951	& zweig, zseaicel(:,:,iobs), zext )
	952
	953	seaicedatqc%rmod(jobs,1) = zext(1)
	954
	955	END DO
	956
	957	! Deallocate the data for interpolation
	958	DEALLOCATE( &
	959	& igrdi, &
	960	& igrdj, &
	961	& zglam, &
	962	& zgphi, &
	963	& zmask, &
	964	& zseaicel &
	965	& )
	966
	967	seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
	968
	969	END SUBROUTINE obs_seaice_opt
	970
	971	SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
	972	& pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
	973	& ld_dailyav )
	974	!!-----------------------------------------------------------------------
	975	!!
	976	!! * ROUTINE obs_vel_opt *
	977	!!
	978	!! ** Purpose : Compute the model counterpart of velocity profile
	979	!! data by interpolating from the model grid to the
	980	!! observation point.
	981	!!
	982	!! ** Method : Linearly interpolate zonal and meridional components of velocity
	983	!! to each observation point using the model values at the corners of
	984	!! the surrounding grid box. The model velocity components are on a
	985	!! staggered C- grid.
	986	!!
	987	!! For velocity data from the TAO array, the model equivalent is
	988	!! a daily mean velocity field. So, we first compute
	989	!! the mean, then interpolate only at the end of the day.
	990	!!
	991	!! ** Action :
	992	!!
	993	!! History :
	994	!! ! 07-03 (K. Mogensen) : Temperature and Salinity profiles
	995	!! ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
	996	!!-----------------------------------------------------------------------
	997
	998	!! * Modules used
	999	USE obs_profiles_def ! Definition of storage space for profile obs.
	1000
	1001	IMPLICIT NONE
	1002
	1003	!! * Arguments
	1004	TYPE(obs_prof), INTENT(INOUT) :: &
	1005	& prodatqc ! Subset of profile data not failing screening
	1006	INTEGER, INTENT(IN) :: kt ! Time step
	1007	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	1008	INTEGER, INTENT(IN) :: kpj
	1009	INTEGER, INTENT(IN) :: kpk
	1010	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	1011	! (kit000-1 = restart time)
	1012	INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header)
	1013	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
	1014	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
	1015	REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
	1016	& pun, & ! Model zonal component of velocity
	1017	& pvn, & ! Model meridional component of velocity
	1018	& pumask, & ! Land-sea mask
	1019	& pvmask ! Land-sea mask
	1020	REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
	1021	& pgdept ! Model array of depth levels
	1022	LOGICAL, INTENT(IN) :: ld_dailyav
	1023
	1024	!! * Local declarations
	1025	INTEGER :: ji
	1026	INTEGER :: jj
	1027	INTEGER :: jk
	1028	INTEGER :: jobs
	1029	INTEGER :: inrc
	1030	INTEGER :: ipro
	1031	INTEGER :: idayend
	1032	INTEGER :: ista
	1033	INTEGER :: iend
	1034	INTEGER :: iobs
	1035	INTEGER, DIMENSION(imaxavtypes) :: &
	1036	& idailyavtypes
	1037	REAL(KIND=wp) :: zlam
	1038	REAL(KIND=wp) :: zphi
	1039	REAL(KIND=wp) :: zdaystp
	1040	REAL(KIND=wp), DIMENSION(kpk) :: &
[2576]	1041	& zobsmasku, &
	1042	& zobsmaskv, &
	1043	& zobsmask, &
	1044	& zobsk, &
[2128]	1045	& zobs2k
	1046	REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
	1047	& zweigu,zweigv
	1048	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
	1049	& zumask, zvmask, &
	1050	& zintu, &
	1051	& zintv, &
	1052	& zinmu, &
	1053	& zinmv
	1054	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	1055	& zglamu, zglamv, &
	1056	& zgphiu, zgphiv
	1057	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
	1058	& igrdiu, &
	1059	& igrdju, &
	1060	& igrdiv, &
	1061	& igrdjv
	1062
	1063	!------------------------------------------------------------------------
	1064	! Local initialization
	1065	!------------------------------------------------------------------------
	1066	! ... Record and data counters
	1067	inrc = kt - kit000 + 2
	1068	ipro = prodatqc%npstp(inrc)
	1069
	1070	! Initialize daily mean for first timestep
	1071	idayend = MOD( kt - kit000 + 1, kdaystp )
	1072
	1073	! Added kt == 0 test to catch restart case
	1074	IF ( idayend == 1 .OR. kt == 0) THEN
	1075	IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
	1076	prodatqc%vdmean(:,:,:,1) = 0.0
	1077	prodatqc%vdmean(:,:,:,2) = 0.0
	1078	ENDIF
	1079
	1080	! Increment the zonal velocity field for computing daily mean
	1081	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
	1082	! Increment the meridional velocity field for computing daily mean
	1083	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
	1084
	1085	! Compute the daily mean at the end of day
	1086	zdaystp = 1.0 / REAL( kdaystp )
	1087	IF ( idayend == 0 ) THEN
	1088	prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
	1089	prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
	1090	ENDIF
	1091
	1092	! Get the data for interpolation
	1093	ALLOCATE( &
	1094	& igrdiu(2,2,ipro), &
	1095	& igrdju(2,2,ipro), &
	1096	& igrdiv(2,2,ipro), &
	1097	& igrdjv(2,2,ipro), &
	1098	& zglamu(2,2,ipro), zglamv(2,2,ipro), &
	1099	& zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
	1100	& zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
	1101	& zintu(2,2,kpk,ipro), &
	1102	& zintv(2,2,kpk,ipro) &
	1103	& )
	1104
	1105	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	1106	iobs = jobs - prodatqc%nprofup
	1107	igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
	1108	igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
	1109	igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
	1110	igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
	1111	igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
	1112	igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
	1113	igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
	1114	igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
	1115	igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
	1116	igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
	1117	igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
	1118	igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
	1119	igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
	1120	igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
	1121	igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
	1122	igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
	1123	END DO
	1124
	1125	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
	1126	CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
	1127	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
	1128	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
	1129
	1130	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
	1131	CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
	1132	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
	1133	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
	1134
	1135	! At the end of the day also get interpolated means
	1136	IF ( idayend == 0 ) THEN
	1137
	1138	ALLOCATE( &
	1139	& zinmu(2,2,kpk,ipro), &
	1140	& zinmv(2,2,kpk,ipro) &
	1141	& )
	1142
	1143	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
	1144	& prodatqc%vdmean(:,:,:,1), zinmu )
	1145	CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
	1146	& prodatqc%vdmean(:,:,:,2), zinmv )
	1147
	1148	ENDIF
	1149
	1150	! loop over observations
	1151
	1152	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	1153
	1154	iobs = jobs - prodatqc%nprofup
	1155
	1156	IF ( kt /= prodatqc%mstp(jobs) ) THEN
	1157
	1158	IF(lwp) THEN
	1159	WRITE(numout,*)
	1160	WRITE(numout,*) ' E R R O R : Observation', &
	1161	& ' time step is not consistent with the', &
	1162	& ' model time step'
	1163	WRITE(numout,*) ' ========='
	1164	WRITE(numout,*)
	1165	WRITE(numout,*) ' Record = ', jobs, &
	1166	& ' kt = ', kt, &
	1167	& ' mstp = ', prodatqc%mstp(jobs), &
	1168	& ' ntyp = ', prodatqc%ntyp(jobs)
	1169	ENDIF
	1170	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
	1171	ENDIF
	1172
	1173	zlam = prodatqc%rlam(jobs)
	1174	zphi = prodatqc%rphi(jobs)
	1175
[2576]	1176	! Initialize observation masks
	1177
	1178	zobsmasku(:) = 0.0
	1179	zobsmaskv(:) = 0.0
	1180
[2128]	1181	! Horizontal weights and vertical mask
	1182
	1183	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
	1184
	1185	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
	1186	& zglamu(:,:,iobs), zgphiu(:,:,iobs), &
[2576]	1187	& zumask(:,:,:,iobs), zweigu, zobsmasku )
[2128]	1188
	1189	ENDIF
	1190
	1191
	1192	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
	1193
	1194	CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, &
	1195	& zglamv(:,:,iobs), zgphiv(:,:,iobs), &
[2576]	1196	& zvmask(:,:,:,iobs), zweigv, zobsmasku )
[2128]	1197
	1198	ENDIF
	1199
[2576]	1200	! Ensure that the vertical mask on u and v are consistent.
[2128]	1201
[2576]	1202	zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
	1203
[2128]	1204	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
	1205
	1206	zobsk(:) = obfillflt
	1207
	1208	IF ( ld_dailyav ) THEN
	1209
	1210	IF ( idayend == 0 ) THEN
	1211
	1212	! Daily averaged data
	1213
	1214	CALL obs_int_h2d( kpk, kpk, &
	1215	& zweigu, zinmu(:,:,:,iobs), zobsk )
	1216
	1217
	1218	ELSE
	1219
	1220	CALL ctl_stop( ' A nonzero' // &
	1221	& ' number of U profile data should' // &
	1222	& ' only occur at the end of a given day' )
	1223
	1224	ENDIF
	1225
	1226	ELSE
	1227
	1228	! Point data
	1229
	1230	CALL obs_int_h2d( kpk, kpk, &
	1231	& zweigu, zintu(:,:,:,iobs), zobsk )
	1232
	1233	ENDIF
	1234
	1235	!-------------------------------------------------------------
	1236	! Compute vertical second-derivative of the interpolating
	1237	! polynomial at obs points
	1238	!-------------------------------------------------------------
	1239
	1240	IF ( k1dint == 1 ) THEN
	1241	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
	1242	& pgdept, zobsmask )
	1243	ENDIF
	1244
	1245	!-----------------------------------------------------------------
	1246	! Vertical interpolation to the observation point
	1247	!-----------------------------------------------------------------
	1248	ista = prodatqc%npvsta(jobs,1)
	1249	iend = prodatqc%npvend(jobs,1)
	1250	CALL obs_int_z1d( kpk, &
	1251	& prodatqc%var(1)%mvk(ista:iend), &
	1252	& k1dint, iend - ista + 1, &
	1253	& prodatqc%var(1)%vdep(ista:iend), &
	1254	& zobsk, zobs2k, &
	1255	& prodatqc%var(1)%vmod(ista:iend), &
	1256	& pgdept, zobsmask )
	1257
	1258	ENDIF
	1259
	1260	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
	1261
	1262	zobsk(:) = obfillflt
	1263
	1264	IF ( ld_dailyav ) THEN
	1265
	1266	IF ( idayend == 0 ) THEN
	1267
	1268	! Daily averaged data
	1269
	1270	CALL obs_int_h2d( kpk, kpk, &
	1271	& zweigv, zinmv(:,:,:,iobs), zobsk )
	1272
	1273	ELSE
	1274
	1275	CALL ctl_stop( ' A nonzero' // &
	1276	& ' number of V profile data should' // &
	1277	& ' only occur at the end of a given day' )
	1278
	1279	ENDIF
	1280
	1281	ELSE
	1282
	1283	! Point data
	1284
	1285	CALL obs_int_h2d( kpk, kpk, &
	1286	& zweigv, zintv(:,:,:,iobs), zobsk )
	1287
	1288	ENDIF
	1289
	1290
	1291	!-------------------------------------------------------------
	1292	! Compute vertical second-derivative of the interpolating
	1293	! polynomial at obs points
	1294	!-------------------------------------------------------------
	1295
	1296	IF ( k1dint == 1 ) THEN
	1297	CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
	1298	& pgdept, zobsmask )
	1299	ENDIF
	1300
	1301	!----------------------------------------------------------------
	1302	! Vertical interpolation to the observation point
	1303	!----------------------------------------------------------------
	1304	ista = prodatqc%npvsta(jobs,2)
	1305	iend = prodatqc%npvend(jobs,2)
	1306	CALL obs_int_z1d( kpk, &
	1307	& prodatqc%var(2)%mvk(ista:iend),&
	1308	& k1dint, iend - ista + 1, &
	1309	& prodatqc%var(2)%vdep(ista:iend),&
	1310	& zobsk, zobs2k, &
	1311	& prodatqc%var(2)%vmod(ista:iend),&
	1312	& pgdept, zobsmask )
	1313
	1314	ENDIF
	1315
	1316	END DO
	1317
	1318	! Deallocate the data for interpolation
	1319	DEALLOCATE( &
	1320	& igrdiu, &
	1321	& igrdju, &
	1322	& igrdiv, &
	1323	& igrdjv, &
	1324	& zglamu, zglamv, &
	1325	& zgphiu, zgphiv, &
	1326	& zumask, zvmask, &
	1327	& zintu, &
	1328	& zintv &
	1329	& )
	1330	! At the end of the day also get interpolated means
	1331	IF ( idayend == 0 ) THEN
	1332	DEALLOCATE( &
	1333	& zinmu, &
	1334	& zinmv &
	1335	& )
	1336	ENDIF
	1337
	1338	prodatqc%nprofup = prodatqc%nprofup + ipro
	1339
	1340	END SUBROUTINE obs_vel_opt
	1341
	1342	END MODULE obs_oper
	1343

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

Download in other formats: