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 NEMO/trunk/src/OCE/OBS – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/OBS/obs_oper.F90 @ 13802

Last change on this file since 13802 was 13295, checked in by acc, 4 years ago
Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision
Property svn:keywords set to `Id`
File size: 37.5 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	!!----------------------------------------------------------------------
[6140]	9	!! obs_prof_opt : Compute the model counterpart of profile data
	10	!! obs_surf_opt : Compute the model counterpart of surface data
[2128]	11	!!----------------------------------------------------------------------
[9490]	12	USE obs_inter_sup ! Interpolation support
	13	USE obs_inter_h2d, ONLY : obs_int_h2d, obs_int_h2d_init ! Horizontal interpolation to the obs pt
	14	USE obs_averg_h2d, ONLY : obs_avg_h2d, obs_avg_h2d_init, obs_max_fpsize ! Horizontal averaging to the obs footprint
	15	USE obs_inter_z1d, ONLY : obs_int_z1d, obs_int_z1d_spl ! Vertical interpolation to the obs pt
	16	USE obs_const , ONLY : obfillflt ! Obs fill value
	17	USE dom_oce, ONLY : glamt, glamf, gphit, gphif ! lat/lon of ocean grid-points
	18	USE lib_mpp, ONLY : ctl_warn, ctl_stop ! Warning and stopping routines
	19	USE sbcdcy, ONLY : sbc_dcy, nday_qsr ! For calculation of where it is night-time
	20	USE obs_grid, ONLY : obs_level_search
	21	!
	22	USE par_kind , ONLY : wp ! Precision variables
	23	USE in_out_manager ! I/O manager
[2128]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
[9490]	28	PUBLIC obs_prof_opt !: Compute the model counterpart of profile obs
	29	PUBLIC obs_surf_opt !: Compute the model counterpart of surface obs
[2128]	30
[9490]	31	INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 !: Max number of daily avgd obs types
[2128]	32
[12377]	33	!! * Substitutions
	34	# include "do_loop_substitute.h90"
[2287]	35	!!----------------------------------------------------------------------
[9570]	36	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2287]	37	!! $Id$
[10068]	38	!! Software governed by the CeCILL license (see ./LICENSE)
[2287]	39	!!----------------------------------------------------------------------
[2128]	40	CONTAINS
	41
[6140]	42	SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
	43	& kit000, kdaystp, &
[9023]	44	& pvar1, pvar2, pgdept, pgdepw, &
	45	& pmask1, pmask2, &
[6140]	46	& plam1, plam2, pphi1, pphi2, &
	47	& k1dint, k2dint, kdailyavtypes )
[2128]	48	!!-----------------------------------------------------------------------
	49	!! * ROUTINE obs_pro_opt *
	50	!!
	51	!! ** Purpose : Compute the model counterpart of profiles
	52	!! data by interpolating from the model grid to the
	53	!! observation point.
	54	!!
	55	!! ** Method : Linearly interpolate to each observation point using
	56	!! the model values at the corners of the surrounding grid box.
	57	!!
	58	!! First, a vertical profile of horizontally interpolated model
[6140]	59	!! now values is computed at the obs (lon, lat) point.
[2128]	60	!! Several horizontal interpolation schemes are available:
	61	!! - distance-weighted (great circle) (k2dint = 0)
	62	!! - distance-weighted (small angle) (k2dint = 1)
	63	!! - bilinear (geographical grid) (k2dint = 2)
	64	!! - bilinear (quadrilateral grid) (k2dint = 3)
	65	!! - polynomial (quadrilateral grid) (k2dint = 4)
	66	!!
[6140]	67	!! Next, the vertical profile is interpolated to the
[2128]	68	!! data depth points. Two vertical interpolation schemes are
	69	!! available:
	70	!! - linear (k1dint = 0)
	71	!! - Cubic spline (k1dint = 1)
	72	!!
	73	!! For the cubic spline the 2nd derivative of the interpolating
	74	!! polynomial is computed before entering the vertical interpolation
	75	!! routine.
	76	!!
[6140]	77	!! If the logical is switched on, the model equivalent is
[2128]	78	!! a daily mean model temperature field. So, we first compute
	79	!! the mean, then interpolate only at the end of the day.
	80	!!
[6140]	81	!! Note: in situ temperature observations must be converted
[2128]	82	!! to potential temperature (the model variable) prior to
	83	!! assimilation.
	84	!!
	85	!! ** Action :
	86	!!
	87	!! History :
	88	!! ! 97-11 (A. Weaver, S. Ricci, N. Daget)
	89	!! ! 06-03 (G. Smith) NEMOVAR migration
	90	!! ! 06-10 (A. Weaver) Cleanup
	91	!! ! 07-01 (K. Mogensen) Merge of temperature and salinity
	92	!! ! 07-03 (K. Mogensen) General handling of profiles
[6140]	93	!! ! 15-02 (M. Martin) Combined routine for all profile types
[9023]	94	!! ! 17-02 (M. Martin) Include generalised vertical coordinate changes
[2128]	95	!!-----------------------------------------------------------------------
	96	USE obs_profiles_def ! Definition of storage space for profile obs.
	97
	98	IMPLICIT NONE
	99
[9490]	100	TYPE(obs_prof), INTENT(inout) :: prodatqc ! Subset of profile data passing QC
	101	INTEGER , INTENT(in ) :: kt ! Time step
	102	INTEGER , INTENT(in ) :: kpi, kpj, kpk ! Model grid parameters
	103	INTEGER , INTENT(in ) :: kit000 ! Number of the first time step (kit000-1 = restart time)
	104	INTEGER , INTENT(in ) :: k1dint ! Vertical interpolation type (see header)
	105	INTEGER , INTENT(in ) :: k2dint ! Horizontal interpolation type (see header)
	106	INTEGER , INTENT(in ) :: kdaystp ! Number of time steps per day
	107	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pvar1 , pvar2 ! Model field 1 and 2
	108	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pmask1, pmask2 ! Land-sea mask 1 and 2
	109	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj) :: plam1 , plam2 ! Model longitude 1 and 2
	110	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj) :: pphi1 , pphi2 ! Model latitudes 1 and 2
	111	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pgdept, pgdepw ! depth of T and W levels
	112	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: kdailyavtypes ! Types for daily averages
[6140]	113
[2128]	114	!! * Local declarations
	115	INTEGER :: ji
	116	INTEGER :: jj
	117	INTEGER :: jk
	118	INTEGER :: jobs
	119	INTEGER :: inrc
	120	INTEGER :: ipro
	121	INTEGER :: idayend
	122	INTEGER :: ista
	123	INTEGER :: iend
	124	INTEGER :: iobs
[9023]	125	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
	126	INTEGER :: inum_obs
[2128]	127	INTEGER, DIMENSION(imaxavtypes) :: &
	128	& idailyavtypes
[6140]	129	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
	130	& igrdi1, &
	131	& igrdi2, &
	132	& igrdj1, &
	133	& igrdj2
[9023]	134	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
	135
[2128]	136	REAL(KIND=wp) :: zlam
	137	REAL(KIND=wp) :: zphi
	138	REAL(KIND=wp) :: zdaystp
	139	REAL(KIND=wp), DIMENSION(kpk) :: &
[6140]	140	& zobsmask1, &
	141	& zobsmask2, &
[2128]	142	& zobsk, &
	143	& zobs2k
[9023]	144	REAL(KIND=wp), DIMENSION(2,2,1) :: &
[6140]	145	& zweig1, &
[9023]	146	& zweig2, &
	147	& zweig
[2128]	148	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
[6140]	149	& zmask1, &
	150	& zmask2, &
[9023]	151	& zint1, &
	152	& zint2, &
	153	& zinm1, &
	154	& zinm2, &
	155	& zgdept, &
	156	& zgdepw
[2128]	157	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
[6140]	158	& zglam1, &
	159	& zglam2, &
	160	& zgphi1, &
	161	& zgphi2
[9023]	162	REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2
	163	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
	164
[6140]	165	LOGICAL :: ld_dailyav
[2128]	166
	167	!------------------------------------------------------------------------
	168	! Local initialization
	169	!------------------------------------------------------------------------
[6140]	170	! Record and data counters
[2128]	171	inrc = kt - kit000 + 2
	172	ipro = prodatqc%npstp(inrc)
[6140]	173
[2128]	174	! Daily average types
[6140]	175	ld_dailyav = .FALSE.
[2128]	176	IF ( PRESENT(kdailyavtypes) ) THEN
	177	idailyavtypes(:) = kdailyavtypes(:)
[6140]	178	IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
[2128]	179	ELSE
	180	idailyavtypes(:) = -1
	181	ENDIF
	182
[6140]	183	! Daily means are calculated for values over timesteps:
	184	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
[2128]	185	idayend = MOD( kt - kit000 + 1, kdaystp )
	186
[6140]	187	IF ( ld_dailyav ) THEN
	188
	189	! Initialize daily mean for first timestep of the day
	190	IF ( idayend == 1 .OR. kt == 0 ) THEN
[13295]	191	DO_3D( 1, 1, 1, 1, 1, jpk )
[12377]	192	prodatqc%vdmean(ji,jj,jk,1) = 0.0
	193	prodatqc%vdmean(ji,jj,jk,2) = 0.0
	194	END_3D
[6140]	195	ENDIF
[2128]	196
[13295]	197	DO_3D( 1, 1, 1, 1, 1, jpk )
[12377]	198	! Increment field 1 for computing daily mean
	199	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
	200	& + pvar1(ji,jj,jk)
	201	! Increment field 2 for computing daily mean
	202	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
	203	& + pvar2(ji,jj,jk)
	204	END_3D
[6140]	205
	206	! Compute the daily mean at the end of day
	207	zdaystp = 1.0 / REAL( kdaystp )
	208	IF ( idayend == 0 ) THEN
	209	IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
	210	CALL FLUSH(numout)
[13295]	211	DO_3D( 1, 1, 1, 1, 1, jpk )
[12377]	212	prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
	213	& * zdaystp
	214	prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
	215	& * zdaystp
	216	END_3D
[6140]	217	ENDIF
	218
[2128]	219	ENDIF
	220
	221	! Get the data for interpolation
	222	ALLOCATE( &
[6140]	223	& igrdi1(2,2,ipro), &
	224	& igrdi2(2,2,ipro), &
	225	& igrdj1(2,2,ipro), &
	226	& igrdj2(2,2,ipro), &
	227	& zglam1(2,2,ipro), &
	228	& zglam2(2,2,ipro), &
	229	& zgphi1(2,2,ipro), &
	230	& zgphi2(2,2,ipro), &
	231	& zmask1(2,2,kpk,ipro), &
	232	& zmask2(2,2,kpk,ipro), &
[9023]	233	& zint1(2,2,kpk,ipro), &
	234	& zint2(2,2,kpk,ipro), &
	235	& zgdept(2,2,kpk,ipro), &
	236	& zgdepw(2,2,kpk,ipro) &
[2128]	237	& )
	238
	239	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	240	iobs = jobs - prodatqc%nprofup
[6140]	241	igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
	242	igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
	243	igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
	244	igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
	245	igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
	246	igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
	247	igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
	248	igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
	249	igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
	250	igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
	251	igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
	252	igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
	253	igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
	254	igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
	255	igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
	256	igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
[2128]	257	END DO
	258
[9023]	259	! Initialise depth arrays
	260	zgdept(:,:,:,:) = 0.0
	261	zgdepw(:,:,:,:) = 0.0
	262
[6140]	263	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
	264	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
	265	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
	266	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 )
	267
	268	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
	269	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
	270	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
	271	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 )
[2128]	272
[9023]	273	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )
	274	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )
	275
[2128]	276	! At the end of the day also get interpolated means
[6140]	277	IF ( ld_dailyav .AND. idayend == 0 ) THEN
[2128]	278
	279	ALLOCATE( &
[6140]	280	& zinm1(2,2,kpk,ipro), &
	281	& zinm2(2,2,kpk,ipro) &
[2128]	282	& )
	283
[6140]	284	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
	285	& prodatqc%vdmean(:,:,:,1), zinm1 )
	286	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
	287	& prodatqc%vdmean(:,:,:,2), zinm2 )
[2128]	288
	289	ENDIF
	290
[9023]	291	! Return if no observations to process
	292	! Has to be done after comm commands to ensure processors
	293	! stay in sync
	294	IF ( ipro == 0 ) RETURN
	295
[2128]	296	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	297
	298	iobs = jobs - prodatqc%nprofup
	299
	300	IF ( kt /= prodatqc%mstp(jobs) ) THEN
[6140]	301
[2128]	302	IF(lwp) THEN
	303	WRITE(numout,*)
	304	WRITE(numout,*) ' E R R O R : Observation', &
	305	& ' time step is not consistent with the', &
	306	& ' model time step'
	307	WRITE(numout,*) ' ========='
	308	WRITE(numout,*)
	309	WRITE(numout,*) ' Record = ', jobs, &
	310	& ' kt = ', kt, &
	311	& ' mstp = ', prodatqc%mstp(jobs), &
	312	& ' ntyp = ', prodatqc%ntyp(jobs)
	313	ENDIF
	314	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
	315	ENDIF
[6140]	316
[2128]	317	zlam = prodatqc%rlam(jobs)
	318	zphi = prodatqc%rphi(jobs)
[6140]	319
[9023]	320	! Horizontal weights
	321	! Masked values are calculated later.
[6140]	322	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
[2128]	323
[9023]	324	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
[6140]	325	& zglam1(:,:,iobs), zgphi1(:,:,iobs), &
[9023]	326	& zmask1(:,:,1,iobs), zweig1, zmsk_1 )
[2128]	327
	328	ENDIF
	329
[6140]	330	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
	331
[9023]	332	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
[6140]	333	& zglam2(:,:,iobs), zgphi2(:,:,iobs), &
[12377]	334	& zmask2(:,:,1,iobs), zweig2, zmsk_2)
[6140]	335
	336	ENDIF
	337
[2128]	338	IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
	339
	340	zobsk(:) = obfillflt
	341
[6140]	342	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
[2128]	343
	344	IF ( idayend == 0 ) THEN
[6140]	345	! Daily averaged data
[2128]	346
[9023]	347	! vertically interpolate all 4 corners
	348	ista = prodatqc%npvsta(jobs,1)
	349	iend = prodatqc%npvend(jobs,1)
	350	inum_obs = iend - ista + 1
	351	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
[6140]	352
[9023]	353	DO iin=1,2
	354	DO ijn=1,2
[6140]	355
[9023]	356	IF ( k1dint == 1 ) THEN
	357	CALL obs_int_z1d_spl( kpk, &
	358	& zinm1(iin,ijn,:,iobs), &
	359	& zobs2k, zgdept(iin,ijn,:,iobs), &
	360	& zmask1(iin,ijn,:,iobs))
	361	ENDIF
	362
	363	CALL obs_level_search(kpk, &
	364	& zgdept(iin,ijn,:,iobs), &
	365	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
	366	& iv_indic)
[2128]	367
[9023]	368	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
	369	& prodatqc%var(1)%vdep(ista:iend), &
	370	& zinm1(iin,ijn,:,iobs), &
	371	& zobs2k, interp_corner(iin,ijn,:), &
	372	& zgdept(iin,ijn,:,iobs), &
	373	& zmask1(iin,ijn,:,iobs))
	374
	375	ENDDO
	376	ENDDO
[2128]	377
[9023]	378	ENDIF !idayend
[6140]	379
[9023]	380	ELSE
[6140]	381
[9023]	382	! Point data
	383
	384	! vertically interpolate all 4 corners
	385	ista = prodatqc%npvsta(jobs,1)
	386	iend = prodatqc%npvend(jobs,1)
	387	inum_obs = iend - ista + 1
	388	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
	389	DO iin=1,2
	390	DO ijn=1,2
	391
	392	IF ( k1dint == 1 ) THEN
	393	CALL obs_int_z1d_spl( kpk, &
	394	& zint1(iin,ijn,:,iobs),&
	395	& zobs2k, zgdept(iin,ijn,:,iobs), &
	396	& zmask1(iin,ijn,:,iobs))
	397
	398	ENDIF
	399
	400	CALL obs_level_search(kpk, &
	401	& zgdept(iin,ijn,:,iobs),&
	402	& inum_obs, prodatqc%var(1)%vdep(ista:iend), &
	403	& iv_indic)
[2128]	404
[9023]	405	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
	406	& prodatqc%var(1)%vdep(ista:iend), &
	407	& zint1(iin,ijn,:,iobs), &
	408	& zobs2k,interp_corner(iin,ijn,:), &
	409	& zgdept(iin,ijn,:,iobs), &
	410	& zmask1(iin,ijn,:,iobs) )
	411
	412	ENDDO
	413	ENDDO
	414
	415	ENDIF
[2128]	416
[9023]	417	!-------------------------------------------------------------
	418	! Compute the horizontal interpolation for every profile level
	419	!-------------------------------------------------------------
	420
	421	DO ikn=1,inum_obs
	422	iend=ista+ikn-1
	423
	424	zweig(:,:,1) = 0._wp
	425
	426	! This code forces the horizontal weights to be
	427	! zero IF the observation is below the bottom of the
	428	! corners of the interpolation nodes, Or if it is in
	429	! the mask. This is important for observations near
	430	! steep bathymetry
	431	DO iin=1,2
	432	DO ijn=1,2
	433
	434	depth_loop1: DO ik=kpk,2,-1
	435	IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN
	436
	437	zweig(iin,ijn,1) = &
	438	& zweig1(iin,ijn,1) * &
	439	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
	440	& - prodatqc%var(1)%vdep(iend)),0._wp)
	441
	442	EXIT depth_loop1
	443
	444	ENDIF
	445
	446	ENDDO depth_loop1
	447
	448	ENDDO
	449	ENDDO
	450
	451	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
	452	& prodatqc%var(1)%vmod(iend:iend) )
	453
	454	! Set QC flag for any observations found below the bottom
	455	! needed as the check here is more strict than that in obs_prep
	456	IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
	457
	458	ENDDO
	459
	460	DEALLOCATE(interp_corner,iv_indic)
	461
	462	ENDIF
	463
	464	! For the second variable
[2128]	465	IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
	466
	467	zobsk(:) = obfillflt
	468
	469	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
	470
	471	IF ( idayend == 0 ) THEN
[6140]	472	! Daily averaged data
[2128]	473
[9023]	474	! vertically interpolate all 4 corners
	475	ista = prodatqc%npvsta(jobs,2)
	476	iend = prodatqc%npvend(jobs,2)
	477	inum_obs = iend - ista + 1
	478	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
[2128]	479
[9023]	480	DO iin=1,2
	481	DO ijn=1,2
[6140]	482
[9023]	483	IF ( k1dint == 1 ) THEN
	484	CALL obs_int_z1d_spl( kpk, &
	485	& zinm2(iin,ijn,:,iobs), &
	486	& zobs2k, zgdept(iin,ijn,:,iobs), &
	487	& zmask2(iin,ijn,:,iobs))
	488	ENDIF
	489
	490	CALL obs_level_search(kpk, &
	491	& zgdept(iin,ijn,:,iobs), &
	492	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
	493	& iv_indic)
[2128]	494
[9023]	495	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
	496	& prodatqc%var(2)%vdep(ista:iend), &
	497	& zinm2(iin,ijn,:,iobs), &
	498	& zobs2k, interp_corner(iin,ijn,:), &
	499	& zgdept(iin,ijn,:,iobs), &
	500	& zmask2(iin,ijn,:,iobs))
	501
	502	ENDDO
	503	ENDDO
[2128]	504
[9023]	505	ENDIF !idayend
[2128]	506
[9023]	507	ELSE
[6140]	508
[9023]	509	! Point data
	510
	511	! vertically interpolate all 4 corners
	512	ista = prodatqc%npvsta(jobs,2)
	513	iend = prodatqc%npvend(jobs,2)
	514	inum_obs = iend - ista + 1
	515	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
	516	DO iin=1,2
	517	DO ijn=1,2
	518
	519	IF ( k1dint == 1 ) THEN
	520	CALL obs_int_z1d_spl( kpk, &
	521	& zint2(iin,ijn,:,iobs),&
	522	& zobs2k, zgdept(iin,ijn,:,iobs), &
	523	& zmask2(iin,ijn,:,iobs))
	524
	525	ENDIF
	526
	527	CALL obs_level_search(kpk, &
	528	& zgdept(iin,ijn,:,iobs),&
	529	& inum_obs, prodatqc%var(2)%vdep(ista:iend), &
	530	& iv_indic)
[6140]	531
[9023]	532	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
	533	& prodatqc%var(2)%vdep(ista:iend), &
	534	& zint2(iin,ijn,:,iobs), &
	535	& zobs2k,interp_corner(iin,ijn,:), &
	536	& zgdept(iin,ijn,:,iobs), &
	537	& zmask2(iin,ijn,:,iobs) )
	538
	539	ENDDO
	540	ENDDO
	541
	542	ENDIF
[2128]	543
[9023]	544	!-------------------------------------------------------------
	545	! Compute the horizontal interpolation for every profile level
	546	!-------------------------------------------------------------
	547
	548	DO ikn=1,inum_obs
	549	iend=ista+ikn-1
	550
	551	zweig(:,:,1) = 0._wp
	552
	553	! This code forces the horizontal weights to be
	554	! zero IF the observation is below the bottom of the
	555	! corners of the interpolation nodes, Or if it is in
	556	! the mask. This is important for observations near
	557	! steep bathymetry
	558	DO iin=1,2
	559	DO ijn=1,2
	560
	561	depth_loop2: DO ik=kpk,2,-1
	562	IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN
	563
	564	zweig(iin,ijn,1) = &
	565	& zweig2(iin,ijn,1) * &
	566	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
	567	& - prodatqc%var(2)%vdep(iend)),0._wp)
	568
	569	EXIT depth_loop2
[2128]	570
[9023]	571	ENDIF
[6140]	572
[9023]	573	ENDDO depth_loop2
	574
	575	ENDDO
	576	ENDDO
	577
	578	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
	579	& prodatqc%var(2)%vmod(iend:iend) )
	580
	581	! Set QC flag for any observations found below the bottom
	582	! needed as the check here is more strict than that in obs_prep
	583	IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
	584
	585	ENDDO
	586
	587	DEALLOCATE(interp_corner,iv_indic)
	588
	589	ENDIF
	590
	591	ENDDO
	592
[2128]	593	! Deallocate the data for interpolation
	594	DEALLOCATE( &
[6140]	595	& igrdi1, &
	596	& igrdi2, &
	597	& igrdj1, &
	598	& igrdj2, &
	599	& zglam1, &
	600	& zglam2, &
	601	& zgphi1, &
	602	& zgphi2, &
	603	& zmask1, &
	604	& zmask2, &
	605	& zint1, &
[9023]	606	& zint2, &
	607	& zgdept, &
	608	& zgdepw &
[2128]	609	& )
[6140]	610
[2128]	611	! At the end of the day also get interpolated means
[6140]	612	IF ( ld_dailyav .AND. idayend == 0 ) THEN
[2128]	613	DEALLOCATE( &
[6140]	614	& zinm1, &
	615	& zinm2 &
[2128]	616	& )
	617	ENDIF
	618
	619	prodatqc%nprofup = prodatqc%nprofup + ipro
[6140]	620
	621	END SUBROUTINE obs_prof_opt
	622
[9023]	623	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
	624	& kit000, kdaystp, psurf, psurfmask, &
	625	& k2dint, ldnightav, plamscl, pphiscl, &
	626	& lindegrees )
[2128]	627
	628	!!-----------------------------------------------------------------------
	629	!!
[6140]	630	!! * ROUTINE obs_surf_opt *
[2128]	631	!!
[6140]	632	!! ** Purpose : Compute the model counterpart of surface
[2128]	633	!! data by interpolating from the model grid to the
	634	!! observation point.
	635	!!
	636	!! ** Method : Linearly interpolate to each observation point using
	637	!! the model values at the corners of the surrounding grid box.
	638	!!
[6140]	639	!! The new model value is first computed at the obs (lon, lat) point.
[2128]	640	!!
	641	!! Several horizontal interpolation schemes are available:
	642	!! - distance-weighted (great circle) (k2dint = 0)
	643	!! - distance-weighted (small angle) (k2dint = 1)
	644	!! - bilinear (geographical grid) (k2dint = 2)
	645	!! - bilinear (quadrilateral grid) (k2dint = 3)
	646	!! - polynomial (quadrilateral grid) (k2dint = 4)
	647	!!
[9023]	648	!! Two horizontal averaging schemes are also available:
	649	!! - weighted radial footprint (k2dint = 5)
	650	!! - weighted rectangular footprint (k2dint = 6)
[2128]	651	!!
[9023]	652	!!
[2128]	653	!! ** Action :
	654	!!
	655	!! History :
[6140]	656	!! ! 07-03 (A. Weaver)
	657	!! ! 15-02 (M. Martin) Combined routine for surface types
[9023]	658	!! ! 17-03 (M. Martin) Added horizontal averaging options
[2128]	659	!!-----------------------------------------------------------------------
	660	USE obs_surf_def ! Definition of storage space for surface observations
	661
	662	IMPLICIT NONE
	663
	664	TYPE(obs_surf), INTENT(INOUT) :: &
[6140]	665	& surfdataqc ! Subset of surface data passing QC
[2128]	666	INTEGER, INTENT(IN) :: kt ! Time step
	667	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	668	INTEGER, INTENT(IN) :: kpj
	669	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	670	! (kit000-1 = restart time)
[6140]	671	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
[2128]	672	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
[6140]	673	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
	674	& psurf, & ! Model surface field
	675	& psurfmask ! Land-sea mask
	676	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
[9023]	677	REAL(KIND=wp), INTENT(IN) :: &
	678	& plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions
	679	& pphiscl ! This is the full width (rather than half-width)
	680	LOGICAL, INTENT(IN) :: &
	681	& lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
[3651]	682
[2128]	683	!! * Local declarations
	684	INTEGER :: ji
	685	INTEGER :: jj
	686	INTEGER :: jobs
	687	INTEGER :: inrc
[6140]	688	INTEGER :: isurf
[2128]	689	INTEGER :: iobs
[9023]	690	INTEGER :: imaxifp, imaxjfp
	691	INTEGER :: imodi, imodj
[3651]	692	INTEGER :: idayend
[6140]	693	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
[9023]	694	& igrdi, &
	695	& igrdj, &
	696	& igrdip1, &
	697	& igrdjp1
[3651]	698	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
[6140]	699	& icount_night, &
[3651]	700	& imask_night
[6140]	701	REAL(wp) :: zlam
	702	REAL(wp) :: zphi
	703	REAL(wp), DIMENSION(1) :: zext, zobsmask
	704	REAL(wp) :: zdaystp
[2128]	705	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
[9023]	706	& zweig, &
[6140]	707	& zmask, &
	708	& zsurf, &
	709	& zsurfm, &
[9023]	710	& zsurftmp, &
[6140]	711	& zglam, &
[9023]	712	& zgphi, &
	713	& zglamf, &
	714	& zgphif
	715
[6140]	716	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
	717	& zintmp, &
	718	& zouttmp, &
	719	& zmeanday ! to compute model sst in region of 24h daylight (pole)
[2128]	720
[6140]	721	!------------------------------------------------------------------------
[2128]	722	! Local initialization
[6140]	723	!------------------------------------------------------------------------
	724	! Record and data counters
[2128]	725	inrc = kt - kit000 + 2
[6140]	726	isurf = surfdataqc%nsstp(inrc)
[2128]	727
[9023]	728	! Work out the maximum footprint size for the
	729	! interpolation/averaging in model grid-points - has to be even.
	730
	731	CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
	732
	733
[6140]	734	IF ( ldnightav ) THEN
[3651]	735
	736	! Initialize array for night mean
[6140]	737	IF ( kt == 0 ) THEN
	738	ALLOCATE ( icount_night(kpi,kpj) )
	739	ALLOCATE ( imask_night(kpi,kpj) )
	740	ALLOCATE ( zintmp(kpi,kpj) )
	741	ALLOCATE ( zouttmp(kpi,kpj) )
	742	ALLOCATE ( zmeanday(kpi,kpj) )
	743	nday_qsr = -1 ! initialisation flag for nbc_dcy
	744	ENDIF
[3651]	745
[6140]	746	! Night-time means are calculated for night-time values over timesteps:
	747	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
	748	idayend = MOD( kt - kit000 + 1, kdaystp )
[3651]	749
[6140]	750	! Initialize night-time mean for first timestep of the day
	751	IF ( idayend == 1 .OR. kt == 0 ) THEN
[13295]	752	DO_2D( 1, 1, 1, 1 )
[12377]	753	surfdataqc%vdmean(ji,jj) = 0.0
	754	zmeanday(ji,jj) = 0.0
	755	icount_night(ji,jj) = 0
	756	END_2D
[6140]	757	ENDIF
[3651]	758
[6140]	759	zintmp(:,:) = 0.0
	760	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
	761	imask_night(:,:) = INT( zouttmp(:,:) )
[3651]	762
[13295]	763	DO_2D( 1, 1, 1, 1 )
[12377]	764	! Increment the temperature field for computing night mean and counter
	765	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
	766	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
	767	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
	768	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
	769	END_2D
[3651]	770
[6140]	771	! Compute the night-time mean at the end of the day
	772	zdaystp = 1.0 / REAL( kdaystp )
	773	IF ( idayend == 0 ) THEN
	774	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
[13295]	775	DO_2D( 1, 1, 1, 1 )
[12377]	776	! Test if "no night" point
	777	IF ( icount_night(ji,jj) > 0 ) THEN
	778	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
	779	& / REAL( icount_night(ji,jj) )
	780	ELSE
	781	!At locations where there is no night (e.g. poles),
	782	! calculate daily mean instead of night-time mean.
	783	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
	784	ENDIF
	785	END_2D
[6140]	786	ENDIF
	787
[3651]	788	ENDIF
	789
[2128]	790	! Get the data for interpolation
[6140]	791
[2128]	792	ALLOCATE( &
[9023]	793	& zweig(imaxifp,imaxjfp,1), &
	794	& igrdi(imaxifp,imaxjfp,isurf), &
	795	& igrdj(imaxifp,imaxjfp,isurf), &
	796	& zglam(imaxifp,imaxjfp,isurf), &
	797	& zgphi(imaxifp,imaxjfp,isurf), &
	798	& zmask(imaxifp,imaxjfp,isurf), &
	799	& zsurf(imaxifp,imaxjfp,isurf), &
	800	& zsurftmp(imaxifp,imaxjfp,isurf), &
	801	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
	802	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
	803	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
	804	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
[2128]	805	& )
[6140]	806
	807	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
	808	iobs = jobs - surfdataqc%nsurfup
[9023]	809	DO ji = 0, imaxifp
	810	imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
[9490]	811	!
[9023]	812	!Deal with wrap around in longitude
	813	IF ( imodi < 1 ) imodi = imodi + jpiglo
	814	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
[9490]	815	!
[9023]	816	DO jj = 0, imaxjfp
	817	imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
	818	!If model values are out of the domain to the north/south then
	819	!set them to be the edge of the domain
	820	IF ( imodj < 1 ) imodj = 1
	821	IF ( imodj > jpjglo ) imodj = jpjglo
[9490]	822	!
[9023]	823	igrdip1(ji+1,jj+1,iobs) = imodi
	824	igrdjp1(ji+1,jj+1,iobs) = imodj
[9490]	825	!
[9023]	826	IF ( ji >= 1 .AND. jj >= 1 ) THEN
	827	igrdi(ji,jj,iobs) = imodi
	828	igrdj(ji,jj,iobs) = imodj
	829	ENDIF
[9490]	830	!
[9023]	831	END DO
	832	END DO
[2128]	833	END DO
[6140]	834
[9023]	835	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[2128]	836	& igrdi, igrdj, glamt, zglam )
[9023]	837	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[2128]	838	& igrdi, igrdj, gphit, zgphi )
[9023]	839	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[6140]	840	& igrdi, igrdj, psurfmask, zmask )
[9023]	841	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[6140]	842	& igrdi, igrdj, psurf, zsurf )
[9023]	843	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
	844	& igrdip1, igrdjp1, glamf, zglamf )
	845	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
	846	& igrdip1, igrdjp1, gphif, zgphif )
[3651]	847
	848	! At the end of the day get interpolated means
[9023]	849	IF ( idayend == 0 .AND. ldnightav ) THEN
[3651]	850
[9023]	851	ALLOCATE( &
	852	& zsurfm(imaxifp,imaxjfp,isurf) &
	853	& )
[3651]	854
[9023]	855	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
	856	& surfdataqc%vdmean(:,:), zsurfm )
[3651]	857
	858	ENDIF
	859
[2128]	860	! Loop over observations
[6140]	861	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
[2128]	862
[6140]	863	iobs = jobs - surfdataqc%nsurfup
	864
	865	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
	866
[2128]	867	IF(lwp) THEN
	868	WRITE(numout,*)
	869	WRITE(numout,*) ' E R R O R : Observation', &
	870	& ' time step is not consistent with the', &
	871	& ' model time step'
	872	WRITE(numout,*) ' ========='
	873	WRITE(numout,*)
	874	WRITE(numout,*) ' Record = ', jobs, &
	875	& ' kt = ', kt, &
[6140]	876	& ' mstp = ', surfdataqc%mstp(jobs), &
	877	& ' ntyp = ', surfdataqc%ntyp(jobs)
[2128]	878	ENDIF
[6140]	879	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
	880
[2128]	881	ENDIF
[6140]	882
	883	zlam = surfdataqc%rlam(jobs)
	884	zphi = surfdataqc%rphi(jobs)
	885
	886	IF ( ldnightav .AND. idayend == 0 ) THEN
	887	! Night-time averaged data
[9023]	888	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
[6140]	889	ELSE
[9023]	890	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
[6140]	891	ENDIF
[3651]	892
[9023]	893	IF ( k2dint <= 4 ) THEN
	894
	895	! Get weights to interpolate the model value to the observation point
	896	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
	897	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	898	& zmask(:,:,iobs), zweig, zobsmask )
	899
	900	! Interpolate the model value to the observation point
	901	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
	902
	903	ELSE
	904
	905	! Get weights to average the model SLA to the observation footprint
	906	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
	907	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	908	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
	909	& zmask(:,:,iobs), plamscl, pphiscl, &
[12377]	910	& lindegrees, zweig )
[9023]	911
	912	! Average the model SST to the observation footprint
	913	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
	914	& zweig, zsurftmp(:,:,iobs), zext )
	915
	916	ENDIF
	917
[6140]	918	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
	919	! ... Remove the MDT from the SSH at the observation point to get the SLA
	920	surfdataqc%rext(jobs,1) = zext(1)
	921	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
[3651]	922	ELSE
[6140]	923	surfdataqc%rmod(jobs,1) = zext(1)
	924	ENDIF
[9023]	925
	926	IF ( zext(1) == obfillflt ) THEN
	927	! If the observation value is a fill value, set QC flag to bad
	928	surfdataqc%nqc(jobs) = 4
	929	ENDIF
[3651]	930
[6140]	931	END DO
[3651]	932
[2128]	933	! Deallocate the data for interpolation
	934	DEALLOCATE( &
[9023]	935	& zweig, &
[2128]	936	& igrdi, &
	937	& igrdj, &
	938	& zglam, &
	939	& zgphi, &
	940	& zmask, &
[9023]	941	& zsurf, &
	942	& zsurftmp, &
	943	& zglamf, &
	944	& zgphif, &
	945	& igrdip1,&
	946	& igrdjp1 &
[2128]	947	& )
[3651]	948
[6140]	949	! At the end of the day also deallocate night-time mean array
[9023]	950	IF ( idayend == 0 .AND. ldnightav ) THEN
	951	DEALLOCATE( &
	952	& zsurfm &
	953	& )
[3651]	954	ENDIF
[9490]	955	!
[6140]	956	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
[9490]	957	!
[6140]	958	END SUBROUTINE obs_surf_opt
[2128]	959
[9490]	960	!!======================================================================
[2128]	961	END MODULE obs_oper

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

Download in other formats: