Context Navigation

obs_oper.F90 @ 10712

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

Download in other formats: