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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_medusa/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 8466

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

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

Download in other formats: