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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 2715

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

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

Download in other formats: