Context Navigation

obs_oper.F90 @ 14548

Last change on this file since 14548 was 14056, checked in by ayoung, 3 years ago
Adding branch for ticket #2567 to trunk.
Property svn:keywords set to `Id`
File size: 30.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
[14056]	42	SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
	43	& kit000, kdaystp, kvar, &
	44	& pvar, pgdept, pgdepw, &
	45	& pmask, &
	46	& plam, pphi, &
[6140]	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
[14056]	107	INTEGER , INTENT(in ) :: kvar ! Number of variables in prodatqc
	108	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pvar ! Model field
	109	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pmask ! Land-sea mask
	110	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj) :: plam ! Model longitude
	111	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj) :: pphi ! Model latitudes
[9490]	112	REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pgdept, pgdepw ! depth of T and W levels
	113	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: kdailyavtypes ! Types for daily averages
[6140]	114
[2128]	115	!! * Local declarations
	116	INTEGER :: ji
	117	INTEGER :: jj
	118	INTEGER :: jk
	119	INTEGER :: jobs
	120	INTEGER :: inrc
	121	INTEGER :: ipro
	122	INTEGER :: idayend
	123	INTEGER :: ista
	124	INTEGER :: iend
	125	INTEGER :: iobs
[9023]	126	INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes
	127	INTEGER :: inum_obs
[2128]	128	INTEGER, DIMENSION(imaxavtypes) :: &
	129	& idailyavtypes
[6140]	130	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
[14056]	131	& igrdi, &
	132	& igrdj
[9023]	133	INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
	134
[2128]	135	REAL(KIND=wp) :: zlam
	136	REAL(KIND=wp) :: zphi
	137	REAL(KIND=wp) :: zdaystp
	138	REAL(KIND=wp), DIMENSION(kpk) :: &
[14056]	139	& zobsk, &
[2128]	140	& zobs2k
[9023]	141	REAL(KIND=wp), DIMENSION(2,2,1) :: &
[6140]	142	& zweig1, &
[9023]	143	& zweig
[2128]	144	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
[14056]	145	& zmask, &
	146	& zint, &
	147	& zinm, &
[9023]	148	& zgdept, &
	149	& zgdepw
[2128]	150	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
[14056]	151	& zglam, &
	152	& zgphi
	153	REAL(KIND=wp), DIMENSION(1) :: zmsk
[9023]	154	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
	155
[6140]	156	LOGICAL :: ld_dailyav
[2128]	157
	158	!------------------------------------------------------------------------
	159	! Local initialization
	160	!------------------------------------------------------------------------
[6140]	161	! Record and data counters
[2128]	162	inrc = kt - kit000 + 2
	163	ipro = prodatqc%npstp(inrc)
[6140]	164
[2128]	165	! Daily average types
[6140]	166	ld_dailyav = .FALSE.
[2128]	167	IF ( PRESENT(kdailyavtypes) ) THEN
	168	idailyavtypes(:) = kdailyavtypes(:)
[6140]	169	IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
[2128]	170	ELSE
	171	idailyavtypes(:) = -1
	172	ENDIF
	173
[6140]	174	! Daily means are calculated for values over timesteps:
	175	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
[2128]	176	idayend = MOD( kt - kit000 + 1, kdaystp )
	177
[6140]	178	IF ( ld_dailyav ) THEN
	179
	180	! Initialize daily mean for first timestep of the day
	181	IF ( idayend == 1 .OR. kt == 0 ) THEN
[13295]	182	DO_3D( 1, 1, 1, 1, 1, jpk )
[14056]	183	prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
[12377]	184	END_3D
[6140]	185	ENDIF
[2128]	186
[13295]	187	DO_3D( 1, 1, 1, 1, 1, jpk )
[12377]	188	! Increment field 1 for computing daily mean
[14056]	189	prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
	190	& + pvar(ji,jj,jk)
[12377]	191	END_3D
[6140]	192
	193	! Compute the daily mean at the end of day
	194	zdaystp = 1.0 / REAL( kdaystp )
	195	IF ( idayend == 0 ) THEN
	196	IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
	197	CALL FLUSH(numout)
[13295]	198	DO_3D( 1, 1, 1, 1, 1, jpk )
[14056]	199	prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
	200	& * zdaystp
[12377]	201	END_3D
[6140]	202	ENDIF
	203
[2128]	204	ENDIF
	205
	206	! Get the data for interpolation
	207	ALLOCATE( &
[14056]	208	& igrdi(2,2,ipro), &
	209	& igrdj(2,2,ipro), &
	210	& zglam(2,2,ipro), &
	211	& zgphi(2,2,ipro), &
	212	& zmask(2,2,kpk,ipro), &
	213	& zint(2,2,kpk,ipro), &
	214	& zgdept(2,2,kpk,ipro), &
	215	& zgdepw(2,2,kpk,ipro) &
[2128]	216	& )
	217
	218	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	219	iobs = jobs - prodatqc%nprofup
[14056]	220	igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
	221	igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
	222	igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
	223	igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
	224	igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
	225	igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
	226	igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
	227	igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
[2128]	228	END DO
	229
[9023]	230	! Initialise depth arrays
	231	zgdept(:,:,:,:) = 0.0
	232	zgdepw(:,:,:,:) = 0.0
	233
[14056]	234	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
	235	CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
	236	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask )
	237	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint )
[2128]	238
[14056]	239	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )
	240	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )
[9023]	241
[2128]	242	! At the end of the day also get interpolated means
[6140]	243	IF ( ld_dailyav .AND. idayend == 0 ) THEN
[2128]	244
[14056]	245	ALLOCATE( zinm(2,2,kpk,ipro) )
[2128]	246
[14056]	247	CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
	248	& prodatqc%vdmean(:,:,:,kvar), zinm )
[2128]	249
	250	ENDIF
	251
[9023]	252	! Return if no observations to process
	253	! Has to be done after comm commands to ensure processors
	254	! stay in sync
	255	IF ( ipro == 0 ) RETURN
	256
[2128]	257	DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
	258
	259	iobs = jobs - prodatqc%nprofup
	260
	261	IF ( kt /= prodatqc%mstp(jobs) ) THEN
[6140]	262
[2128]	263	IF(lwp) THEN
	264	WRITE(numout,*)
	265	WRITE(numout,*) ' E R R O R : Observation', &
	266	& ' time step is not consistent with the', &
	267	& ' model time step'
	268	WRITE(numout,*) ' ========='
	269	WRITE(numout,*)
	270	WRITE(numout,*) ' Record = ', jobs, &
	271	& ' kt = ', kt, &
	272	& ' mstp = ', prodatqc%mstp(jobs), &
	273	& ' ntyp = ', prodatqc%ntyp(jobs)
	274	ENDIF
	275	CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
	276	ENDIF
[6140]	277
[2128]	278	zlam = prodatqc%rlam(jobs)
	279	zphi = prodatqc%rphi(jobs)
[6140]	280
[9023]	281	! Horizontal weights
	282	! Masked values are calculated later.
[14056]	283	IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
[2128]	284
[9023]	285	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
[14056]	286	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	287	& zmask(:,:,1,iobs), zweig1, zmsk )
[2128]	288
	289	ENDIF
	290
[14056]	291	IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
[6140]	292
[2128]	293	zobsk(:) = obfillflt
	294
[6140]	295	IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
[2128]	296
	297	IF ( idayend == 0 ) THEN
[6140]	298	! Daily averaged data
[2128]	299
[9023]	300	! vertically interpolate all 4 corners
[14056]	301	ista = prodatqc%npvsta(jobs,kvar)
	302	iend = prodatqc%npvend(jobs,kvar)
[9023]	303	inum_obs = iend - ista + 1
	304	ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))
[6140]	305
[9023]	306	DO iin=1,2
	307	DO ijn=1,2
[6140]	308
[9023]	309	IF ( k1dint == 1 ) THEN
	310	CALL obs_int_z1d_spl( kpk, &
[14056]	311	& zinm(iin,ijn,:,iobs), &
[9023]	312	& zobs2k, zgdept(iin,ijn,:,iobs), &
[14056]	313	& zmask(iin,ijn,:,iobs))
[9023]	314	ENDIF
	315
	316	CALL obs_level_search(kpk, &
	317	& zgdept(iin,ijn,:,iobs), &
[14056]	318	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
[9023]	319	& iv_indic)
[2128]	320
[9023]	321	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
[14056]	322	& prodatqc%var(kvar)%vdep(ista:iend), &
	323	& zinm(iin,ijn,:,iobs), &
[9023]	324	& zobs2k, interp_corner(iin,ijn,:), &
	325	& zgdept(iin,ijn,:,iobs), &
[14056]	326	& zmask(iin,ijn,:,iobs))
[9023]	327
	328	ENDDO
	329	ENDDO
[2128]	330
[9023]	331	ENDIF !idayend
[6140]	332
[9023]	333	ELSE
[6140]	334
[9023]	335	! Point data
	336
	337	! vertically interpolate all 4 corners
[14056]	338	ista = prodatqc%npvsta(jobs,kvar)
	339	iend = prodatqc%npvend(jobs,kvar)
[9023]	340	inum_obs = iend - ista + 1
	341	ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))
	342	DO iin=1,2
	343	DO ijn=1,2
	344
	345	IF ( k1dint == 1 ) THEN
	346	CALL obs_int_z1d_spl( kpk, &
[14056]	347	& zint(iin,ijn,:,iobs),&
[9023]	348	& zobs2k, zgdept(iin,ijn,:,iobs), &
[14056]	349	& zmask(iin,ijn,:,iobs))
[9023]	350
	351	ENDIF
	352
	353	CALL obs_level_search(kpk, &
	354	& zgdept(iin,ijn,:,iobs),&
[14056]	355	& inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &
[9023]	356	& iv_indic)
[2128]	357
[9023]	358	CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &
[14056]	359	& prodatqc%var(kvar)%vdep(ista:iend), &
	360	& zint(iin,ijn,:,iobs), &
[9023]	361	& zobs2k,interp_corner(iin,ijn,:), &
	362	& zgdept(iin,ijn,:,iobs), &
[14056]	363	& zmask(iin,ijn,:,iobs) )
[9023]	364
	365	ENDDO
	366	ENDDO
	367
	368	ENDIF
[2128]	369
[9023]	370	!-------------------------------------------------------------
	371	! Compute the horizontal interpolation for every profile level
	372	!-------------------------------------------------------------
	373
	374	DO ikn=1,inum_obs
	375	iend=ista+ikn-1
	376
	377	zweig(:,:,1) = 0._wp
	378
	379	! This code forces the horizontal weights to be
	380	! zero IF the observation is below the bottom of the
	381	! corners of the interpolation nodes, Or if it is in
	382	! the mask. This is important for observations near
	383	! steep bathymetry
	384	DO iin=1,2
	385	DO ijn=1,2
	386
[14056]	387	depth_loop: DO ik=kpk,2,-1
	388	IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN
[9023]	389
	390	zweig(iin,ijn,1) = &
	391	& zweig1(iin,ijn,1) * &
	392	& MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &
[14056]	393	& - prodatqc%var(kvar)%vdep(iend)),0._wp)
[9023]	394
[14056]	395	EXIT depth_loop
[9023]	396
	397	ENDIF
	398
[14056]	399	ENDDO depth_loop
[9023]	400
	401	ENDDO
	402	ENDDO
	403
	404	CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &
[14056]	405	& prodatqc%var(kvar)%vmod(iend:iend) )
[9023]	406
	407	! Set QC flag for any observations found below the bottom
	408	! needed as the check here is more strict than that in obs_prep
[14056]	409	IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
[9023]	410
	411	ENDDO
	412
	413	DEALLOCATE(interp_corner,iv_indic)
	414
[14056]	415	ENDIF
[9023]	416
	417	ENDDO
	418
[2128]	419	! Deallocate the data for interpolation
[14056]	420	DEALLOCATE( &
	421	& igrdi, &
	422	& igrdj, &
	423	& zglam, &
	424	& zgphi, &
	425	& zmask, &
	426	& zint, &
[9023]	427	& zgdept, &
	428	& zgdepw &
[2128]	429	& )
[6140]	430
[2128]	431	! At the end of the day also get interpolated means
[6140]	432	IF ( ld_dailyav .AND. idayend == 0 ) THEN
[14056]	433	DEALLOCATE( zinm )
[2128]	434	ENDIF
	435
[14056]	436	IF ( kvar == prodatqc%nvar ) THEN
	437	prodatqc%nprofup = prodatqc%nprofup + ipro
	438	ENDIF
[6140]	439
	440	END SUBROUTINE obs_prof_opt
	441
[9023]	442	SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, &
	443	& kit000, kdaystp, psurf, psurfmask, &
	444	& k2dint, ldnightav, plamscl, pphiscl, &
	445	& lindegrees )
[2128]	446
	447	!!-----------------------------------------------------------------------
	448	!!
[6140]	449	!! * ROUTINE obs_surf_opt *
[2128]	450	!!
[6140]	451	!! ** Purpose : Compute the model counterpart of surface
[2128]	452	!! data by interpolating from the model grid to the
	453	!! observation point.
	454	!!
	455	!! ** Method : Linearly interpolate to each observation point using
	456	!! the model values at the corners of the surrounding grid box.
	457	!!
[6140]	458	!! The new model value is first computed at the obs (lon, lat) point.
[2128]	459	!!
	460	!! Several horizontal interpolation schemes are available:
	461	!! - distance-weighted (great circle) (k2dint = 0)
	462	!! - distance-weighted (small angle) (k2dint = 1)
	463	!! - bilinear (geographical grid) (k2dint = 2)
	464	!! - bilinear (quadrilateral grid) (k2dint = 3)
	465	!! - polynomial (quadrilateral grid) (k2dint = 4)
	466	!!
[9023]	467	!! Two horizontal averaging schemes are also available:
	468	!! - weighted radial footprint (k2dint = 5)
	469	!! - weighted rectangular footprint (k2dint = 6)
[2128]	470	!!
[9023]	471	!!
[2128]	472	!! ** Action :
	473	!!
	474	!! History :
[6140]	475	!! ! 07-03 (A. Weaver)
	476	!! ! 15-02 (M. Martin) Combined routine for surface types
[9023]	477	!! ! 17-03 (M. Martin) Added horizontal averaging options
[2128]	478	!!-----------------------------------------------------------------------
	479	USE obs_surf_def ! Definition of storage space for surface observations
	480
	481	IMPLICIT NONE
	482
	483	TYPE(obs_surf), INTENT(INOUT) :: &
[6140]	484	& surfdataqc ! Subset of surface data passing QC
[2128]	485	INTEGER, INTENT(IN) :: kt ! Time step
	486	INTEGER, INTENT(IN) :: kpi ! Model grid parameters
	487	INTEGER, INTENT(IN) :: kpj
	488	INTEGER, INTENT(IN) :: kit000 ! Number of the first time step
	489	! (kit000-1 = restart time)
[6140]	490	INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
[2128]	491	INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header)
[6140]	492	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
	493	& psurf, & ! Model surface field
	494	& psurfmask ! Land-sea mask
	495	LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
[9023]	496	REAL(KIND=wp), INTENT(IN) :: &
	497	& plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions
	498	& pphiscl ! This is the full width (rather than half-width)
	499	LOGICAL, INTENT(IN) :: &
	500	& lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
[3651]	501
[2128]	502	!! * Local declarations
	503	INTEGER :: ji
	504	INTEGER :: jj
	505	INTEGER :: jobs
	506	INTEGER :: inrc
[6140]	507	INTEGER :: isurf
[2128]	508	INTEGER :: iobs
[9023]	509	INTEGER :: imaxifp, imaxjfp
	510	INTEGER :: imodi, imodj
[3651]	511	INTEGER :: idayend
[6140]	512	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
[9023]	513	& igrdi, &
	514	& igrdj, &
	515	& igrdip1, &
	516	& igrdjp1
[3651]	517	INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
[6140]	518	& icount_night, &
[3651]	519	& imask_night
[6140]	520	REAL(wp) :: zlam
	521	REAL(wp) :: zphi
	522	REAL(wp), DIMENSION(1) :: zext, zobsmask
	523	REAL(wp) :: zdaystp
[2128]	524	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
[9023]	525	& zweig, &
[6140]	526	& zmask, &
	527	& zsurf, &
	528	& zsurfm, &
[9023]	529	& zsurftmp, &
[6140]	530	& zglam, &
[9023]	531	& zgphi, &
	532	& zglamf, &
	533	& zgphif
	534
[6140]	535	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
	536	& zintmp, &
	537	& zouttmp, &
	538	& zmeanday ! to compute model sst in region of 24h daylight (pole)
[2128]	539
[6140]	540	!------------------------------------------------------------------------
[2128]	541	! Local initialization
[6140]	542	!------------------------------------------------------------------------
	543	! Record and data counters
[2128]	544	inrc = kt - kit000 + 2
[6140]	545	isurf = surfdataqc%nsstp(inrc)
[2128]	546
[9023]	547	! Work out the maximum footprint size for the
	548	! interpolation/averaging in model grid-points - has to be even.
	549
	550	CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
	551
	552
[6140]	553	IF ( ldnightav ) THEN
[3651]	554
	555	! Initialize array for night mean
[6140]	556	IF ( kt == 0 ) THEN
	557	ALLOCATE ( icount_night(kpi,kpj) )
	558	ALLOCATE ( imask_night(kpi,kpj) )
	559	ALLOCATE ( zintmp(kpi,kpj) )
	560	ALLOCATE ( zouttmp(kpi,kpj) )
	561	ALLOCATE ( zmeanday(kpi,kpj) )
	562	nday_qsr = -1 ! initialisation flag for nbc_dcy
	563	ENDIF
[3651]	564
[6140]	565	! Night-time means are calculated for night-time values over timesteps:
	566	! [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
	567	idayend = MOD( kt - kit000 + 1, kdaystp )
[3651]	568
[6140]	569	! Initialize night-time mean for first timestep of the day
	570	IF ( idayend == 1 .OR. kt == 0 ) THEN
[13295]	571	DO_2D( 1, 1, 1, 1 )
[12377]	572	surfdataqc%vdmean(ji,jj) = 0.0
	573	zmeanday(ji,jj) = 0.0
	574	icount_night(ji,jj) = 0
	575	END_2D
[6140]	576	ENDIF
[3651]	577
[6140]	578	zintmp(:,:) = 0.0
	579	zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
	580	imask_night(:,:) = INT( zouttmp(:,:) )
[3651]	581
[13295]	582	DO_2D( 1, 1, 1, 1 )
[12377]	583	! Increment the temperature field for computing night mean and counter
	584	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
	585	& + psurf(ji,jj) * REAL( imask_night(ji,jj) )
	586	zmeanday(ji,jj) = zmeanday(ji,jj) + psurf(ji,jj)
	587	icount_night(ji,jj) = icount_night(ji,jj) + imask_night(ji,jj)
	588	END_2D
[3651]	589
[6140]	590	! Compute the night-time mean at the end of the day
	591	zdaystp = 1.0 / REAL( kdaystp )
	592	IF ( idayend == 0 ) THEN
	593	IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
[13295]	594	DO_2D( 1, 1, 1, 1 )
[12377]	595	! Test if "no night" point
	596	IF ( icount_night(ji,jj) > 0 ) THEN
	597	surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
	598	& / REAL( icount_night(ji,jj) )
	599	ELSE
	600	!At locations where there is no night (e.g. poles),
	601	! calculate daily mean instead of night-time mean.
	602	surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
	603	ENDIF
	604	END_2D
[6140]	605	ENDIF
	606
[3651]	607	ENDIF
	608
[2128]	609	! Get the data for interpolation
[6140]	610
[2128]	611	ALLOCATE( &
[9023]	612	& zweig(imaxifp,imaxjfp,1), &
	613	& igrdi(imaxifp,imaxjfp,isurf), &
	614	& igrdj(imaxifp,imaxjfp,isurf), &
	615	& zglam(imaxifp,imaxjfp,isurf), &
	616	& zgphi(imaxifp,imaxjfp,isurf), &
	617	& zmask(imaxifp,imaxjfp,isurf), &
	618	& zsurf(imaxifp,imaxjfp,isurf), &
	619	& zsurftmp(imaxifp,imaxjfp,isurf), &
	620	& zglamf(imaxifp+1,imaxjfp+1,isurf), &
	621	& zgphif(imaxifp+1,imaxjfp+1,isurf), &
	622	& igrdip1(imaxifp+1,imaxjfp+1,isurf), &
	623	& igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
[2128]	624	& )
[6140]	625
	626	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
	627	iobs = jobs - surfdataqc%nsurfup
[9023]	628	DO ji = 0, imaxifp
	629	imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
[9490]	630	!
[9023]	631	!Deal with wrap around in longitude
	632	IF ( imodi < 1 ) imodi = imodi + jpiglo
	633	IF ( imodi > jpiglo ) imodi = imodi - jpiglo
[9490]	634	!
[9023]	635	DO jj = 0, imaxjfp
	636	imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
	637	!If model values are out of the domain to the north/south then
	638	!set them to be the edge of the domain
	639	IF ( imodj < 1 ) imodj = 1
	640	IF ( imodj > jpjglo ) imodj = jpjglo
[9490]	641	!
[9023]	642	igrdip1(ji+1,jj+1,iobs) = imodi
	643	igrdjp1(ji+1,jj+1,iobs) = imodj
[9490]	644	!
[9023]	645	IF ( ji >= 1 .AND. jj >= 1 ) THEN
	646	igrdi(ji,jj,iobs) = imodi
	647	igrdj(ji,jj,iobs) = imodj
	648	ENDIF
[9490]	649	!
[9023]	650	END DO
	651	END DO
[2128]	652	END DO
[6140]	653
[9023]	654	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[2128]	655	& igrdi, igrdj, glamt, zglam )
[9023]	656	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[2128]	657	& igrdi, igrdj, gphit, zgphi )
[9023]	658	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[6140]	659	& igrdi, igrdj, psurfmask, zmask )
[9023]	660	CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
[6140]	661	& igrdi, igrdj, psurf, zsurf )
[9023]	662	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
	663	& igrdip1, igrdjp1, glamf, zglamf )
	664	CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
	665	& igrdip1, igrdjp1, gphif, zgphif )
[3651]	666
	667	! At the end of the day get interpolated means
[9023]	668	IF ( idayend == 0 .AND. ldnightav ) THEN
[3651]	669
[9023]	670	ALLOCATE( &
	671	& zsurfm(imaxifp,imaxjfp,isurf) &
	672	& )
[3651]	673
[9023]	674	CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
	675	& surfdataqc%vdmean(:,:), zsurfm )
[3651]	676
	677	ENDIF
	678
[2128]	679	! Loop over observations
[6140]	680	DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
[2128]	681
[6140]	682	iobs = jobs - surfdataqc%nsurfup
	683
	684	IF ( kt /= surfdataqc%mstp(jobs) ) THEN
	685
[2128]	686	IF(lwp) THEN
	687	WRITE(numout,*)
	688	WRITE(numout,*) ' E R R O R : Observation', &
	689	& ' time step is not consistent with the', &
	690	& ' model time step'
	691	WRITE(numout,*) ' ========='
	692	WRITE(numout,*)
	693	WRITE(numout,*) ' Record = ', jobs, &
	694	& ' kt = ', kt, &
[6140]	695	& ' mstp = ', surfdataqc%mstp(jobs), &
	696	& ' ntyp = ', surfdataqc%ntyp(jobs)
[2128]	697	ENDIF
[6140]	698	CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
	699
[2128]	700	ENDIF
[6140]	701
	702	zlam = surfdataqc%rlam(jobs)
	703	zphi = surfdataqc%rphi(jobs)
	704
	705	IF ( ldnightav .AND. idayend == 0 ) THEN
	706	! Night-time averaged data
[9023]	707	zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
[6140]	708	ELSE
[9023]	709	zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
[6140]	710	ENDIF
[3651]	711
[9023]	712	IF ( k2dint <= 4 ) THEN
	713
	714	! Get weights to interpolate the model value to the observation point
	715	CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &
	716	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	717	& zmask(:,:,iobs), zweig, zobsmask )
	718
	719	! Interpolate the model value to the observation point
	720	CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
	721
	722	ELSE
	723
	724	! Get weights to average the model SLA to the observation footprint
	725	CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, &
	726	& zglam(:,:,iobs), zgphi(:,:,iobs), &
	727	& zglamf(:,:,iobs), zgphif(:,:,iobs), &
	728	& zmask(:,:,iobs), plamscl, pphiscl, &
[12377]	729	& lindegrees, zweig )
[9023]	730
	731	! Average the model SST to the observation footprint
	732	CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
	733	& zweig, zsurftmp(:,:,iobs), zext )
	734
	735	ENDIF
	736
[6140]	737	IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
	738	! ... Remove the MDT from the SSH at the observation point to get the SLA
	739	surfdataqc%rext(jobs,1) = zext(1)
	740	surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
[3651]	741	ELSE
[6140]	742	surfdataqc%rmod(jobs,1) = zext(1)
	743	ENDIF
[9023]	744
	745	IF ( zext(1) == obfillflt ) THEN
	746	! If the observation value is a fill value, set QC flag to bad
	747	surfdataqc%nqc(jobs) = 4
	748	ENDIF
[3651]	749
[6140]	750	END DO
[3651]	751
[2128]	752	! Deallocate the data for interpolation
	753	DEALLOCATE( &
[9023]	754	& zweig, &
[2128]	755	& igrdi, &
	756	& igrdj, &
	757	& zglam, &
	758	& zgphi, &
	759	& zmask, &
[9023]	760	& zsurf, &
	761	& zsurftmp, &
	762	& zglamf, &
	763	& zgphif, &
	764	& igrdip1,&
	765	& igrdjp1 &
[2128]	766	& )
[3651]	767
[6140]	768	! At the end of the day also deallocate night-time mean array
[9023]	769	IF ( idayend == 0 .AND. ldnightav ) THEN
	770	DEALLOCATE( &
	771	& zsurfm &
	772	& )
[3651]	773	ENDIF
[9490]	774	!
[6140]	775	surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
[9490]	776	!
[6140]	777	END SUBROUTINE obs_surf_opt
[2128]	778
[9490]	779	!!======================================================================
[2128]	780	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 @ 14548

Download in other formats: