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

Context Navigation

source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 11202

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

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

Download in other formats: