Context Navigation

obs_oper.F90 @ 12779

Last change on this file since 12779 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 37.5 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: