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

Context Navigation

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 7837

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

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

Download in other formats: