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

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 9354

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

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

Download in other formats: