742 lines
34 KiB
Fortran
742 lines
34 KiB
Fortran
module module_NWP_match_to_sat
|
|
|
|
USE omp_lib
|
|
use netcdf
|
|
use module_satalite
|
|
use module_read_NWP
|
|
implicit none
|
|
!include 'netcdf.inc'
|
|
|
|
public :: interpolation, allocate_profiles,&
|
|
deallocate_profiles
|
|
private :: select_obs_region
|
|
integer(kind=4),public,allocatable :: qcmark(:,:)
|
|
real(kind=4),public,allocatable :: T(:,:,:),q(:,:,:),cfrac(:,:,:),ozone(:,:,:)
|
|
real(kind=4),public,allocatable :: cloud(:,:,:,:)
|
|
real(kind=4),public,allocatable :: ts(:,:),T2m(:,:), q2m(:,:), u10m(:,:), v10m(:,:)
|
|
real(kind=4),public,allocatable :: PSFC(:,:),seaice(:,:)
|
|
real(kind=4),allocatable :: CO2(:),CH4(:),N2O(:),CO(:)
|
|
integer,public,parameter :: geos=0,polar=1
|
|
|
|
contains
|
|
!=================================================================================================
|
|
!=================================================================================================
|
|
subroutine allocate_profiles(allocstatus)
|
|
integer :: allocstatus
|
|
allocate(T (nscanp, nscanl, nlevel))
|
|
allocate(q (nscanp, nscanl, nlevel))
|
|
allocate(psfc (nscanp, nscanl))
|
|
allocate(seaice (nscanp, nscanl))
|
|
allocate(T2m (nscanp, nscanl))
|
|
allocate(ts (nscanp, nscanl))
|
|
allocate(q2m (nscanp, nscanl))
|
|
allocate(u10m (nscanp, nscanl))
|
|
allocate(cloud (nscanp, nscanl, 6, nlevel-1))
|
|
allocate(cfrac (nscanp, nscanl, nlevel-1))
|
|
allocate(v10m (nscanp, nscanl))
|
|
allocate(qcmark (nscanp, nscanl))
|
|
allocate(ozone (nscanp, nscanl, nlevel))
|
|
allocate(CO2 (nlevel))
|
|
allocate(CH4 (nlevel))
|
|
allocate(N2O (nlevel))
|
|
allocate(CO (nlevel))
|
|
allocstatus = allocated
|
|
return
|
|
end subroutine allocate_profiles
|
|
|
|
!=================================================================================================
|
|
!=================================================================================================
|
|
subroutine deallocate_profiles(allocstatus)
|
|
integer :: allocstatus
|
|
deallocate(T )
|
|
deallocate(q )
|
|
deallocate(psfc )
|
|
deallocate(seaice )
|
|
deallocate(T2m )
|
|
deallocate(ts )
|
|
deallocate(q2m )
|
|
deallocate(u10m )
|
|
deallocate(cloud )
|
|
deallocate(cfrac )
|
|
deallocate(v10m )
|
|
deallocate(qcmark )
|
|
deallocate(ozone )
|
|
deallocate(CO2 )
|
|
deallocate(CH4 )
|
|
deallocate(N2O )
|
|
deallocate(CO )
|
|
allocstatus = deallocated
|
|
return
|
|
end subroutine deallocate_profiles
|
|
!=================================================================================================
|
|
!=================================================================================================
|
|
subroutine interpolation(time1,time2,time3,period,NWP_resolution,sattype)
|
|
|
|
integer(kind=8) :: time1, time2, time3
|
|
real(kind=8) :: second1, second2, second3
|
|
integer :: sattype
|
|
real(kind=4) :: NWP_resolution
|
|
real(kind=4) :: period
|
|
integer :: i,j,k,nlstart,nlend
|
|
integer,parameter :: judge_N=100000
|
|
integer(kind=4),allocatable,dimension(:,:) :: index_lat,index_lon
|
|
integer(kind=4),allocatable,dimension(:) :: indexx0,indexy0,indexx1,indexy1
|
|
real(kind=4),allocatable,dimension(:,:) :: s1,s2,s3,s4,sumgeo
|
|
real(kind=4),allocatable,dimension(:,:) :: fgeo1,fgeo2,fgeo3,fgeo4
|
|
real(kind=4),allocatable,dimension(:,:) :: fractime1,fractime2
|
|
real(kind=4),allocatable,dimension(:,:) :: dtime1,dtime2,dtime3
|
|
real(kind=4),allocatable,dimension(:) :: lonflat,latflat,lonflatt
|
|
|
|
real(kind=4),allocatable :: Tt1(:,:,:),qt1(:,:,:),cfract1(:,:,:),ozonet1(:,:,:)
|
|
real(kind=4),allocatable :: cloudt1(:,:,:,:)
|
|
real(kind=4),allocatable :: tst1(:,:),T2mt1(:,:), q2mt1(:,:), u10mt1(:,:), v10mt1(:,:)
|
|
real(kind=4),allocatable :: PSFCt1(:,:),seaicet1(:,:)
|
|
|
|
real(kind=4),allocatable :: Tt2(:,:,:),qt2(:,:,:),cfract2(:,:,:),ozonet2(:,:,:)
|
|
real(kind=4),allocatable :: cloudt2(:,:,:,:)
|
|
real(kind=4),allocatable :: tst2(:,:),T2mt2(:,:), q2mt2(:,:), u10mt2(:,:), v10mt2(:,:)
|
|
real(kind=4),allocatable :: PSFCt2(:,:),seaicet2(:,:)
|
|
|
|
real(kind=4),allocatable :: Tt3(:,:,:),qt3(:,:,:),cfract3(:,:,:),ozonet3(:,:,:)
|
|
real(kind=4),allocatable :: cloudt3(:,:,:,:)
|
|
real(kind=4),allocatable :: tst3(:,:),T2mt3(:,:), q2mt3(:,:), u10mt3(:,:), v10mt3(:,:)
|
|
real(kind=4),allocatable :: PSFCt3(:,:),seaicet3(:,:)
|
|
|
|
real(kind=4),dimension(101) :: pressure101,CO2101,CH4101,N2O101,CO101,fact
|
|
real(kind=4) :: debug
|
|
|
|
|
|
allocate(index_lat(nscanp,nscanl))
|
|
allocate(index_lon(nscanp,nscanl))
|
|
allocate(dtime1(nscanp,nscanl))
|
|
allocate(dtime2(nscanp,nscanl))
|
|
allocate(dtime3(nscanp,nscanl))
|
|
allocate(fractime1(nscanp,nscanl))
|
|
allocate(fractime2(nscanp,nscanl))
|
|
allocate(fgeo1(nscanp,nscanl))
|
|
allocate(fgeo2(nscanp,nscanl))
|
|
allocate(fgeo3(nscanp,nscanl))
|
|
allocate(fgeo4(nscanp,nscanl))
|
|
allocate(s1 (nscanp,nscanl))
|
|
allocate(s2 (nscanp,nscanl))
|
|
allocate(s3 (nscanp,nscanl))
|
|
allocate(s4 (nscanp,nscanl))
|
|
allocate(sumgeo(nscanp,nscanl))
|
|
|
|
allocate(Tt1 (nscanp, nscanl, nlevel))
|
|
allocate(ozonet1 (nscanp, nscanl, nlevel))
|
|
allocate(qt1 (nscanp, nscanl, nlevel))
|
|
allocate(psfct1 (nscanp, nscanl))
|
|
allocate(seaicet1 (nscanp, nscanl))
|
|
allocate(T2mt1 (nscanp, nscanl))
|
|
allocate(tst1 (nscanp, nscanl))
|
|
allocate(q2mt1 (nscanp, nscanl))
|
|
allocate(u10mt1 (nscanp, nscanl))
|
|
allocate(cloudt1 (nscanp, nscanl, 6, nlevel-1))
|
|
allocate(cfract1 (nscanp, nscanl, nlevel-1))
|
|
allocate(v10mt1 (nscanp, nscanl))
|
|
|
|
allocate(Tt2 (nscanp, nscanl, nlevel))
|
|
allocate(ozonet2 (nscanp, nscanl, nlevel))
|
|
allocate(qt2 (nscanp, nscanl, nlevel))
|
|
allocate(psfct2 (nscanp, nscanl))
|
|
allocate(seaicet2 (nscanp, nscanl))
|
|
allocate(T2mt2 (nscanp, nscanl))
|
|
allocate(tst2 (nscanp, nscanl))
|
|
allocate(q2mt2 (nscanp, nscanl))
|
|
allocate(u10mt2 (nscanp, nscanl))
|
|
allocate(cloudt2 (nscanp, nscanl, 6, nlevel-1))
|
|
allocate(cfract2 (nscanp, nscanl, nlevel-1))
|
|
allocate(v10mt2 (nscanp, nscanl))
|
|
|
|
if(sattype==polar)then
|
|
allocate(Tt3 (nscanp, nscanl, nlevel))
|
|
allocate(ozonet3 (nscanp, nscanl, nlevel))
|
|
allocate(qt3 (nscanp, nscanl, nlevel))
|
|
allocate(psfct3 (nscanp, nscanl))
|
|
allocate(seaicet3 (nscanp, nscanl))
|
|
allocate(T2mt3 (nscanp, nscanl))
|
|
allocate(tst3 (nscanp, nscanl))
|
|
allocate(q2mt3 (nscanp, nscanl))
|
|
allocate(u10mt3 (nscanp, nscanl))
|
|
allocate(cloudt3 (nscanp, nscanl, 6, nlevel-1))
|
|
allocate(cfract3 (nscanp, nscanl, nlevel-1))
|
|
allocate(v10mt3 (nscanp, nscanl))
|
|
endif
|
|
|
|
Tt1 = -9999.9
|
|
ozonet1 = -9999.9
|
|
qt1 = -9999.9
|
|
psfct1 = -9999.9
|
|
seaicet1 = -9999.9
|
|
T2mt1 = -9999.9
|
|
tst1 = -9999.9
|
|
q2mt1 = -9999.9
|
|
u10mt1 = -9999.9
|
|
cloudt1 = -9999.9
|
|
cfract1 = -9999.9
|
|
v10mt1 = -9999.9
|
|
|
|
Tt2 = -9999.9
|
|
ozonet2 = -9999.9
|
|
qt2 = -9999.9
|
|
psfct2 = -9999.9
|
|
seaicet2 = -9999.9
|
|
T2mt2 = -9999.9
|
|
tst2 = -9999.9
|
|
q2mt2 = -9999.9
|
|
u10mt2 = -9999.9
|
|
cloudt2 = -9999.9
|
|
cfract2 = -9999.9
|
|
v10mt2 = -9999.9
|
|
|
|
if(sattype==polar)then
|
|
Tt3 = -9999.9
|
|
ozonet3 = -9999.9
|
|
qt3 = -9999.9
|
|
psfct3 = -9999.9
|
|
seaicet3 = -9999.9
|
|
T2mt3 = -9999.9
|
|
tst3 = -9999.9
|
|
q2mt3 = -9999.9
|
|
u10mt3 = -9999.9
|
|
cloudt3 = -9999.9
|
|
cfract3 = -9999.9
|
|
v10mt3 = -9999.9
|
|
endif
|
|
|
|
pressure101(:) = (/ &
|
|
.005, .016, .038, .077, &
|
|
.137, .224, .345, .506, .714, .975, 1.297, &
|
|
1.687, 2.153, 2.701, 3.340, 4.077, 4.920, 5.878, &
|
|
6.957, 8.165, 9.512, 11.004, 12.649, 14.456, 16.432, &
|
|
18.585, 20.922, 23.453, 26.183, 29.121, 32.274, 35.651, &
|
|
39.257, 43.100, 47.188, 51.528, 56.126, 60.989, 66.125, &
|
|
71.540, 77.240, 83.231, 89.520, 96.114, 103.017, 110.237, &
|
|
117.777, 125.646, 133.846, 142.385, 151.266, 160.496, 170.078, &
|
|
180.018, 190.320, 200.989, 212.028, 223.441, 235.234, 247.408, &
|
|
259.969, 272.919, 286.262, 300.000, 314.137, 328.675, 343.618, &
|
|
358.966, 374.724, 390.893, 407.474, 424.470, 441.882, 459.712, &
|
|
477.961, 496.630, 515.720, 535.232, 555.167, 575.525, 596.306, &
|
|
617.511, 639.140, 661.192, 683.667, 706.565, 729.886, 753.628, &
|
|
777.790, 802.371, 827.371, 852.788, 878.620, 904.866, 931.524, &
|
|
958.591, 986.067, 1013.948, 1042.232, 1070.917, 1100.000 /)
|
|
CO2101 = (/ 0.374361,0.374362,0.374365,0.374369,0.374375,0.374385,0.374397,0.374413,0.374443,0.374491, &
|
|
0.374551,0.374609,0.374649,0.374681,0.374717,0.374762,0.374824,0.374899,0.374978,0.375054, &
|
|
0.375116,0.375165,0.375208,0.375240,0.375275,0.375313,0.375370,0.375435,0.375532,0.375639, &
|
|
0.375778,0.375926,0.376100,0.376286,0.376482,0.376689,0.376930,0.377232,0.377551,0.377915, &
|
|
0.378302,0.378660,0.378980,0.379313,0.379535,0.379768,0.379977,0.380169,0.380362,0.380544, &
|
|
0.380734,0.380893,0.381048,0.381202,0.381359,0.381515,0.381673,0.381833,0.381989,0.382141, &
|
|
0.382286,0.382413,0.382535,0.382619,0.382694,0.382752,0.382803,0.382843,0.382880,0.382912, &
|
|
0.382946,0.382980,0.383022,0.383066,0.383108,0.383149,0.383180,0.383208,0.383228,0.383246, &
|
|
0.383258,0.383272,0.383284,0.383299,0.383318,0.383339,0.383355,0.383375,0.383414,0.383476, &
|
|
0.383546,0.383617,0.383696,0.383778,0.383853,0.383916,0.383967,0.384003,0.384027,0.384036, &
|
|
0.384038/)
|
|
|
|
do i = 1, 101
|
|
fact(i) = 1.0
|
|
end do
|
|
CO2101 = CO2101*1000*fact
|
|
|
|
CH4101 = (/ 0.150E+00, 0.148E+00, 0.151E+00, 0.162E+00, 0.189E+00, 0.225E+00, 0.260E+00, &
|
|
0.293E+00, 0.331E+00, 0.381E+00, 0.452E+00, 0.531E+00, 0.615E+00, 0.708E+00, &
|
|
0.792E+00, 0.861E+00, 0.921E+00, 0.966E+00, 0.100E+01, 0.103E+01, 0.106E+01, &
|
|
0.110E+01, 0.115E+01, 0.119E+01, 0.123E+01, 0.128E+01, 0.132E+01, 0.136E+01, &
|
|
0.140E+01, 0.142E+01, 0.144E+01, 0.147E+01, 0.150E+01, 0.152E+01, 0.154E+01, &
|
|
0.156E+01, 0.158E+01, 0.160E+01, 0.161E+01, 0.162E+01, 0.163E+01, 0.164E+01, &
|
|
0.165E+01, 0.168E+01, 0.169E+01, 0.171E+01, 0.172E+01, 0.173E+01, 0.174E+01, &
|
|
0.174E+01, 0.175E+01, 0.176E+01, 0.176E+01, 0.176E+01, 0.177E+01, 0.177E+01, &
|
|
0.177E+01, 0.177E+01, 0.178E+01, 0.178E+01, 0.178E+01, 0.178E+01, 0.178E+01, &
|
|
0.178E+01, 0.178E+01, 0.178E+01, 0.178E+01, 0.179E+01, 0.179E+01, 0.179E+01, &
|
|
0.179E+01, 0.179E+01, 0.179E+01, 0.179E+01, 0.179E+01, 0.179E+01, 0.179E+01, &
|
|
0.179E+01, 0.179E+01, 0.179E+01, 0.179E+01, 0.179E+01, 0.180E+01, 0.180E+01, &
|
|
0.180E+01, 0.180E+01, 0.180E+01, 0.180E+01, 0.180E+01, 0.180E+01, 0.180E+01, &
|
|
0.180E+01, 0.181E+01, 0.181E+01, 0.181E+01, 0.181E+01, 0.181E+01, 0.181E+01, &
|
|
0.181E+01, 0.181E+01, 0.181E+01 /)
|
|
CH4101 = CH4101*0.96
|
|
|
|
N2O101 = (/ 0.814E-02, 0.604E-02, 0.452E-02, 0.355E-02, 0.292E-02, 0.209E-02, 0.189E-02, &
|
|
0.220E-02, 0.409E-02, 0.671E-02, 0.972E-02, 0.125E-01, 0.151E-01, 0.200E-01, &
|
|
0.251E-01, 0.315E-01, 0.382E-01, 0.459E-01, 0.537E-01, 0.637E-01, 0.733E-01, &
|
|
0.849E-01, 0.972E-01, 0.108E+00, 0.122E+00, 0.135E+00, 0.148E+00, 0.159E+00, &
|
|
0.169E+00, 0.179E+00, 0.188E+00, 0.195E+00, 0.203E+00, 0.210E+00, 0.218E+00, &
|
|
0.225E+00, 0.232E+00, 0.239E+00, 0.246E+00, 0.253E+00, 0.260E+00, 0.267E+00, &
|
|
0.274E+00, 0.281E+00, 0.287E+00, 0.293E+00, 0.297E+00, 0.302E+00, 0.305E+00, &
|
|
0.308E+00, 0.311E+00, 0.313E+00, 0.315E+00, 0.316E+00, 0.317E+00, 0.319E+00, &
|
|
0.319E+00, 0.320E+00, 0.320E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, &
|
|
0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, &
|
|
0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, &
|
|
0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, &
|
|
0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, &
|
|
0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, 0.321E+00, &
|
|
0.321E+00, 0.321E+00, 0.321E+00 /)
|
|
N2O101 = N2O101*0.992
|
|
|
|
CO101 = (/ 0.844E+00, 0.826E+00, 0.791E+00, 0.735E+00, 0.658E+00, 0.564E+00, 0.442E+00, &
|
|
0.320E+00, 0.245E+00, 0.197E+00, 0.155E+00, 0.116E+00, 0.885E-01, 0.710E-01, &
|
|
0.602E-01, 0.513E-01, 0.439E-01, 0.412E-01, 0.396E-01, 0.377E-01, 0.350E-01, &
|
|
0.327E-01, 0.306E-01, 0.282E-01, 0.264E-01, 0.249E-01, 0.242E-01, 0.237E-01, &
|
|
0.236E-01, 0.236E-01, 0.237E-01, 0.239E-01, 0.244E-01, 0.251E-01, 0.263E-01, &
|
|
0.280E-01, 0.302E-01, 0.334E-01, 0.373E-01, 0.414E-01, 0.462E-01, 0.516E-01, &
|
|
0.574E-01, 0.644E-01, 0.679E-01, 0.725E-01, 0.763E-01, 0.797E-01, 0.833E-01, &
|
|
0.859E-01, 0.888E-01, 0.913E-01, 0.940E-01, 0.971E-01, 0.100E+00, 0.104E+00, &
|
|
0.108E+00, 0.113E+00, 0.117E+00, 0.121E+00, 0.125E+00, 0.129E+00, 0.131E+00, &
|
|
0.134E+00, 0.136E+00, 0.137E+00, 0.138E+00, 0.139E+00, 0.140E+00, 0.141E+00, &
|
|
0.141E+00, 0.141E+00, 0.142E+00, 0.142E+00, 0.143E+00, 0.143E+00, 0.144E+00, &
|
|
0.144E+00, 0.144E+00, 0.145E+00, 0.145E+00, 0.145E+00, 0.146E+00, 0.147E+00, &
|
|
0.148E+00, 0.149E+00, 0.150E+00, 0.151E+00, 0.151E+00, 0.152E+00, 0.154E+00, &
|
|
0.156E+00, 0.159E+00, 0.161E+00, 0.163E+00, 0.166E+00, 0.167E+00, 0.169E+00, &
|
|
0.170E+00, 0.172E+00, 0.174E+00 /)
|
|
CO101 = CO101*0.88
|
|
|
|
call INTINP_INT2PRTTOV(101,nlevel,pressure101,P,CO2101, CO2 )
|
|
call INTINP_INT2PRTTOV(101,nlevel,pressure101,P,CH4101, CH4 )
|
|
call INTINP_INT2PRTTOV(101,nlevel,pressure101,P,N2O101, N2O )
|
|
call INTINP_INT2PRTTOV(101,nlevel,pressure101,P,CO101 , CO )
|
|
|
|
where(lat>90 .or. lat<-90)
|
|
lat=LATP(1,1)-NWP_resolution*(judge_N-1)
|
|
endwhere
|
|
where(lon>360 .or. lon<-180)
|
|
lon=LONP(1,1)+NWP_resolution*(judge_N-1)
|
|
endwhere
|
|
|
|
print*,'-----------After profile interpolation!-----------------'
|
|
|
|
!index_lat = aint((LATP(1,1)-lat)/NWP_resolution)+1
|
|
!index_lon = aint((lon-LONP(1,1))/NWP_resolution)+1
|
|
|
|
index_lat = abs(aint((LATP(1,1)-lat)/NWP_resolution))+1
|
|
index_lon = abs(aint((lon-LONP(1,1))/NWP_resolution))+1
|
|
|
|
qcmark(:,:) = 0
|
|
|
|
debug=0
|
|
do i=1,nscanl
|
|
do j=1,nscanp
|
|
if (index_lat(j,i)>n_lat) then
|
|
qcmark(j,i)=1
|
|
lat(j,i)=90
|
|
end if
|
|
enddo
|
|
enddo
|
|
!where(index_lat>n_lat)
|
|
! qcmark=1
|
|
! lat=90.
|
|
!endwhere
|
|
|
|
do i=1,nscanl
|
|
do j=1,nscanp
|
|
if (index_lon(j,i)>n_lon) then
|
|
qcmark(j,i)=1
|
|
lon(j,i)=360
|
|
end if
|
|
enddo
|
|
enddo
|
|
|
|
!where(index_lon>n_lon)
|
|
! qcmark=1
|
|
! lon=360.
|
|
!endwhere
|
|
|
|
where(time<0)
|
|
qcmark=1
|
|
endwhere
|
|
|
|
where(index_lon>=n_lon)
|
|
index_lon = n_lon-1
|
|
endwhere
|
|
where(index_lat>=n_lat)
|
|
index_lat = n_lat-1
|
|
endwhere
|
|
|
|
print*,'time1=',time1,'time2=',time2,'time3=',time3
|
|
call CALCULATE_TIME2MSECOND(dble(time1),second1)
|
|
call CALCULATE_TIME2MSECOND(dble(time2),second2)
|
|
call CALCULATE_TIME2MSECOND(dble(time3),second3)
|
|
print*,'msecond1=',second1,'msecond2=',second2,'msecond3=',second3
|
|
print*,'time=',time(1,1:5)
|
|
|
|
dtime1 = sngl(time(:,:)-second1)
|
|
dtime2 = sngl(time(:,:)-second2)
|
|
if(sattype==polar)then
|
|
dtime3 = sngl(time(:,:)-second3)
|
|
endif
|
|
|
|
! print*,'dtime1=',dtime1(1,1:5),dtime1(nscanp,1)
|
|
! print*,'dtime2=',dtime2(1,1:5),dtime2(nscanp,1)
|
|
! print*,'dtime3=',dtime3(1,1:5),dtime3(nscanp,1)
|
|
|
|
|
|
! print*,minval(LATP),maxval(LATP),minval(LONP),maxval(LONP) !wrf
|
|
! print*,minval(lat),maxval(lat),minval(lon),maxval(lon) !sat
|
|
|
|
|
|
if(time2==0 .and. time3==0)then
|
|
print*,'single time NWP (not need time interp)'
|
|
where(lat<minval(LATP) .or. lat>maxval(LATP) .or. lon<minval(LONP) .or. lon>maxval(LONP))
|
|
qcmark=1
|
|
endwhere
|
|
do i=1,nscanl
|
|
|
|
call OMP_SET_NUM_THREADS(16)
|
|
|
|
!$OMP PARALLEL DO default(shared)
|
|
do j=1,nscanp
|
|
if(qcmark(j, i)==0)then
|
|
s1(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i), index_lat(j,i) ))**2 + (lon(j,i)-lonp(index_lon(j,i),index_lat(j,i)))**2)
|
|
s2(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i), index_lat(j,i)+1))**2 + (lon(j,i)-lonp(index_lon(j,i),index_lat(j,i)+1))**2)
|
|
s3(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i)+1,index_lat(j,i)+1))**2 + (lon(j,i)-lonp(index_lon(j,i)+1,index_lat(j,i)+1))**2)
|
|
s4(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i)+1,index_lat(j,i) ))**2 + (lon(j,i)-lonp(index_lon(j,i)+1,index_lat(j,i)))**2)
|
|
sumgeo(j,i) = 1/s1(j,i)+1/s2(j,i)+1/s3(j,i)+1/s4(j,i)
|
|
fgeo1(j,i) = 1/s1(j,i)/sumgeo(j,i)
|
|
fgeo2(j,i) = 1/s2(j,i)/sumgeo(j,i)
|
|
fgeo3(j,i) = 1/s3(j,i)/sumgeo(j,i)
|
|
fgeo4(j,i) = 1/s4(j,i)/sumgeo(j,i)
|
|
call select_obs_region(index_lon(j,i),index_lat(j,i),index_lon(j,i)+1,index_lat(j,i)+1,fgeo1(j,i),fgeo2(j,i),fgeo3(j,i),fgeo4(j,i), &
|
|
T1 , q1 , cfrac1 , ozone1 , cloud1 , ts1 , T2m1 , q2m1 , u10m1 , v10m1 , PSFC1 , seaice1, &
|
|
Tt1(j,i,:) , qt1(j,i,:) , cfract1(j,i,:) , ozonet1(j,i,:) , cloudt1(j,i,:,:) , tst1(j,i) , T2mt1(j,i) , q2mt1(j,i) , &
|
|
u10mt1(j,i) , v10mt1(j,i) , PSFCt1(j,i) , seaicet1(j,i) , j , i)
|
|
else
|
|
Tt1(j,i,:) = -9999.9
|
|
qt1(j,i,:) = -9999.9
|
|
cfract1(j,i,:) = -9999.9
|
|
ozonet1(j,i,:) = -9999.9
|
|
cloudt1(j,i,:,:) = -9999.9
|
|
tst1(j,i) = -9999.9
|
|
T2mt1(j,i) = -9999.9
|
|
q2mt1(j,i) = -9999.9
|
|
u10mt1(j,i) = -9999.9
|
|
v10mt1(j,i) = -9999.9
|
|
PSFCt1(j,i) = -9999.9
|
|
seaicet1(j,i) = -9999.9
|
|
endif
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
enddo
|
|
if(time2==0 .and. time3==0)then
|
|
T = Tt1
|
|
q = qt1
|
|
cfrac = cfract1
|
|
ozone = ozonet1
|
|
cloud = cloudt1
|
|
ts = tst1
|
|
T2m = T2mt1
|
|
q2m = q2mt1
|
|
u10m = u10mt1
|
|
v10m = v10mt1
|
|
PSFC = PSFCt1
|
|
seaice = seaicet1
|
|
endif
|
|
|
|
else
|
|
print*,'Three times NWP interp, sumqcmark=',sum(qcmark)
|
|
do i=1,nscanl
|
|
|
|
call OMP_SET_NUM_THREADS(16)
|
|
!$OMP PARALLEL DO default(shared)
|
|
do j=1,nscanp
|
|
if(qcmark(j, i)==0)then
|
|
!print*,'+++++ thread=',OMP_GET_THREAD_NUM(),' i,j=',i,j
|
|
s1(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i), index_lat(j,i) ))**2 + (lon(j,i)-lonp(index_lon(j,i),index_lat(j,i)))**2)
|
|
s2(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i), index_lat(j,i)+1))**2 + (lon(j,i)-lonp(index_lon(j,i),index_lat(j,i)+1))**2)
|
|
s3(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i)+1,index_lat(j,i)+1))**2 + (lon(j,i)-lonp(index_lon(j,i)+1,index_lat(j,i)+1))**2)
|
|
s4(j,i) = sqrt((lat(j,i)-latp(index_lon(j,i)+1,index_lat(j,i) ))**2 + (lon(j,i)-lonp(index_lon(j,i)+1,index_lat(j,i)))**2)
|
|
sumgeo(j,i) = 1/s1(j,i)+1/s2(j,i)+1/s3(j,i)+1/s4(j,i)
|
|
fgeo1(j,i) = 1/s1(j,i)/sumgeo(j,i)
|
|
fgeo2(j,i) = 1/s2(j,i)/sumgeo(j,i)
|
|
fgeo3(j,i) = 1/s3(j,i)/sumgeo(j,i)
|
|
fgeo4(j,i) = 1/s4(j,i)/sumgeo(j,i)
|
|
|
|
call select_obs_region(index_lon(j,i),index_lat(j,i),index_lon(j,i)+1,index_lat(j,i)+1,fgeo1(j,i),fgeo2(j,i),fgeo3(j,i),fgeo4(j,i), &
|
|
T1 , q1 , cfrac1 , ozone1 , cloud1 , ts1 , T2m1 , q2m1 , u10m1 , v10m1 , PSFC1 , seaice1, &
|
|
Tt1(j,i,:) , qt1(j,i,:) , cfract1(j,i,:) , ozonet1(j,i,:) , cloudt1(j,i,:,:) , tst1(j,i) , T2mt1(j,i) , q2mt1(j,i) , &
|
|
u10mt1(j,i) , v10mt1(j,i) , PSFCt1(j,i) , seaicet1(j,i) , j , i)
|
|
call select_obs_region(index_lon(j,i),index_lat(j,i),index_lon(j,i)+1,index_lat(j,i)+1,fgeo1(j,i),fgeo2(j,i),fgeo3(j,i),fgeo4(j,i), &
|
|
T2 , q2 , cfrac2 , ozone2 , cloud2 , ts2 , T2m2 , q2m2 , u10m2 , v10m2 , PSFC2 , seaice2, &
|
|
Tt2(j,i,:) , qt2(j,i,:) , cfract2(j,i,:) , ozonet2(j,i,:) , cloudt2(j,i,:,:) , tst2(j,i) , T2mt2(j,i) , q2mt2(j,i) , u10mt2(j,i) , &
|
|
v10mt2(j,i) , PSFCt2(j,i) , seaicet2(j,i) , j , i)
|
|
if(sattype==polar)then
|
|
call select_obs_region(index_lon(j,i),index_lat(j,i),index_lon(j,i)+1,index_lat(j,i)+1,fgeo1(j,i),fgeo2(j,i),fgeo3(j,i),fgeo4(j,i), &
|
|
T3 , q3 , cfrac3 , ozone3 , cloud3 , ts3 , T2m3 , q2m3 , u10m3 , v10m3 , PSFC3 , seaice3, &
|
|
Tt3(j,i,:) , qt3(j,i,:) , cfract3(j,i,:) , ozonet3(j,i,:) , cloudt3(j,i,:,:) , tst3(j,i) , T2mt3(j,i) , q2mt3(j,i) , &
|
|
u10mt3(j,i) , v10mt3(j,i) , PSFCt3(j,i) , seaicet3(j,i) , j , i)
|
|
endif
|
|
else
|
|
Tt1(j,i,:) = -9999.9
|
|
qt1(j,i,:) = -9999.9
|
|
cfract1(j,i,:) = -9999.9
|
|
ozonet1(j,i,:) = -9999.9
|
|
cloudt1(j,i,:,:) = -9999.9
|
|
tst1(j,i) = -9999.9
|
|
T2mt1(j,i) = -9999.9
|
|
q2mt1(j,i) = -9999.9
|
|
u10mt1(j,i) = -9999.9
|
|
v10mt1(j,i) = -9999.9
|
|
PSFCt1(j,i) = -9999.9
|
|
seaicet1(j,i) = -9999.9
|
|
|
|
Tt2(j,i,:) = -9999.9
|
|
qt2(j,i,:) = -9999.9
|
|
cfract2(j,i,:) = -9999.9
|
|
ozonet2(j,i,:) = -9999.9
|
|
cloudt2(j,i,:,:) = -9999.9
|
|
tst2(j,i) = -9999.9
|
|
T2mt2(j,i) = -9999.9
|
|
q2mt2(j,i) = -9999.9
|
|
u10mt2(j,i) = -9999.9
|
|
v10mt2(j,i) = -9999.9
|
|
PSFCt2(j,i) = -9999.9
|
|
seaicet2(j,i) = -9999.9
|
|
if(sattype==polar)then
|
|
Tt3(j,i,:) = -9999.9
|
|
qt3(j,i,:) = -9999.9
|
|
cfract3(j,i,:) = -9999.9
|
|
ozonet3(j,i,:) = -9999.9
|
|
cloudt3(j,i,:,:) = -9999.9
|
|
tst3(j,i) = -9999.9
|
|
T2mt3(j,i) = -9999.9
|
|
q2mt3(j,i) = -9999.9
|
|
u10mt3(j,i) = -9999.9
|
|
v10mt3(j,i) = -9999.9
|
|
PSFCt3(j,i) = -9999.9
|
|
seaicet3(j,i) = -9999.9
|
|
endif
|
|
endif
|
|
|
|
!print*,'----- thread=',OMP_GET_THREAD_NUM(),' i,j=',i,j,'t2m=',T2mt1(j,i)
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
if(dtime2(1,i)<0)then
|
|
fractime2(:,i) = dtime1(:,i)/period
|
|
fractime1(:,i) = abs(dtime2(:,i)/period)
|
|
T(:,i,:) = Tt1(:,i,:) * fractime1(1,i) + Tt2(:,i,:) * fractime2(1,i)
|
|
q(:,i,:) = qt1(:,i,:) * fractime1(1,i) + qt2(:,i,:) * fractime2(1,i)
|
|
cfrac(:,i,:) = cfract1(:,i,:) * fractime1(1,i) + cfract2(:,i,:) * fractime2(1,i)
|
|
ozone(:,i,:) = ozonet1(:,i,:) * fractime1(1,i) + ozonet2(:,i,:) * fractime2(1,i)
|
|
cloud(:,i,:,:) = cloudt1(:,i,:,:) * fractime1(1,i) + cloudt2(:,i,:,:) * fractime2(1,i)
|
|
ts(:,i) = tst1(:,i) * fractime1(1,i) + tst2(:,i) * fractime2(1,i)
|
|
T2m(:,i) = T2mt1(:,i) * fractime1(1,i) + T2mt2(:,i) * fractime2(1,i)
|
|
q2m(:,i) = q2mt1(:,i) * fractime1(1,i) + q2mt2(:,i) * fractime2(1,i)
|
|
u10m(:,i) = u10mt1(:,i) * fractime1(1,i) + u10mt2(:,i) * fractime2(1,i)
|
|
v10m(:,i) = v10mt1(:,i) * fractime1(1,i) + v10mt2(:,i) * fractime2(1,i)
|
|
PSFC(:,i) = PSFCt1(:,i) * fractime1(1,i) + PSFCt2(:,i) * fractime2(1,i)
|
|
seaice(:,i) = seaicet1(:,i)* fractime1(1,i) + seaicet2(:,i) * fractime2(1,i)
|
|
else
|
|
if(sattype==polar)then
|
|
fractime2(:,i) = dtime2(:,i)/period
|
|
fractime1(:,i) = abs(dtime3(:,i)/period)
|
|
T(:,i,:) = Tt2(:,i,:) * fractime1(1,i) + Tt3(:,i,:) * fractime2(1,i)
|
|
q(:,i,:) = qt2(:,i,:) * fractime1(1,i) + qt3(:,i,:) * fractime2(1,i)
|
|
cfrac(:,i,:) = cfract2(:,i,:) * fractime1(1,i) + cfract3(:,i,:) * fractime2(1,i)
|
|
ozone(:,i,:) = ozonet2(:,i,:) * fractime1(1,i) + ozonet3(:,i,:) * fractime2(1,i)
|
|
cloud(:,i,:,:) = cloudt2(:,i,:,:) * fractime1(1,i) + cloudt3(:,i,:,:) * fractime2(1,i)
|
|
ts(:,i) = tst2(:,i) * fractime1(1,i) + tst3(:,i) * fractime2(1,i)
|
|
T2m(:,i) = T2mt2(:,i) * fractime1(1,i) + T2mt3(:,i) * fractime2(1,i)
|
|
q2m(:,i) = q2mt2(:,i) * fractime1(1,i) + q2mt3(:,i) * fractime2(1,i)
|
|
u10m(:,i) = u10mt2(:,i) * fractime1(1,i) + u10mt3(:,i) * fractime2(1,i)
|
|
v10m(:,i) = v10mt2(:,i) * fractime1(1,i) + v10mt3(:,i) * fractime2(1,i)
|
|
PSFC(:,i) = PSFCt2(:,i) * fractime1(1,i) + PSFCt3(:,i) * fractime2(1,i)
|
|
seaice(:,i) = seaicet2(:,i)* fractime1(1,i) + seaicet3(:,i) * fractime2(1,i)
|
|
endif
|
|
endif
|
|
enddo
|
|
endif
|
|
|
|
print*,'---------------------------After match to satelite!---------------------------------'
|
|
!================================================================================================================
|
|
|
|
|
|
deallocate(index_lat)
|
|
deallocate(index_lon)
|
|
deallocate(dtime1)
|
|
deallocate(dtime2)
|
|
deallocate(dtime3)
|
|
deallocate(fractime1)
|
|
deallocate(fractime2)
|
|
deallocate(fgeo1)
|
|
deallocate(fgeo2)
|
|
deallocate(fgeo3)
|
|
deallocate(fgeo4)
|
|
deallocate(s1)
|
|
deallocate(s2)
|
|
deallocate(s3)
|
|
deallocate(s4)
|
|
deallocate(sumgeo)
|
|
|
|
|
|
deallocate(Tt1 )
|
|
deallocate(ozonet1 )
|
|
deallocate(qt1 )
|
|
deallocate(psfct1 )
|
|
deallocate(seaicet1 )
|
|
deallocate(T2mt1 )
|
|
deallocate(tst1 )
|
|
deallocate(q2mt1 )
|
|
deallocate(u10mt1 )
|
|
deallocate(cloudt1 )
|
|
deallocate(cfract1 )
|
|
deallocate(v10mt1 )
|
|
|
|
deallocate(Tt2 )
|
|
deallocate(ozonet2 )
|
|
deallocate(qt2 )
|
|
deallocate(psfct2 )
|
|
deallocate(seaicet2 )
|
|
deallocate(T2mt2 )
|
|
deallocate(tst2 )
|
|
deallocate(q2mt2 )
|
|
deallocate(u10mt2 )
|
|
deallocate(cloudt2 )
|
|
deallocate(cfract2 )
|
|
deallocate(v10mt2 )
|
|
|
|
if(sattype==polar)then
|
|
deallocate(Tt3 )
|
|
deallocate(ozonet3 )
|
|
deallocate(qt3 )
|
|
deallocate(psfct3 )
|
|
deallocate(seaicet3 )
|
|
deallocate(T2mt3 )
|
|
deallocate(tst3 )
|
|
deallocate(q2mt3 )
|
|
deallocate(u10mt3 )
|
|
deallocate(cloudt3 )
|
|
deallocate(cfract3 )
|
|
deallocate(v10mt3 )
|
|
endif
|
|
|
|
return
|
|
end subroutine interpolation
|
|
!====================================================================================================================================
|
|
!====================================================================================================================================
|
|
subroutine select_obs_region(indexx0, indexy0, indexx1, indexy1, fgeo1, fgeo2, fgeo3, fgeo4, &
|
|
Tin , qin , cfracin , ozonein , cloudin , tsin , T2min , q2min , u10min , v10min , PSFCin , seaicein , &
|
|
Tout, qout, cfracout, ozoneout, cloudout, tsout, T2mout, q2mout, u10mout, v10mout, PSFCout, seaiceout, &
|
|
ix, iy)
|
|
|
|
!====================================================================================================================================
|
|
|
|
integer(kind=4) :: indexx0,indexy0,indexx1,indexy1
|
|
real(kind=4) :: fgeo1 ,fgeo2 ,fgeo3 ,fgeo4
|
|
integer :: i,j,k,ix,iy
|
|
!====================================================================================================================================
|
|
real(kind=4) :: Tin(n_lon, n_lat,nlevel)
|
|
real(kind=4) :: qin(n_lon, n_lat,nlevel)
|
|
real(kind=4) :: cloudin(n_lon, n_lat,6,nlevel-1)
|
|
real(kind=4) :: cfracin(n_lon, n_lat,nlevel-1)
|
|
real(kind=4) :: ozonein(n_lon, n_lat,nlevel)
|
|
real(kind=4) :: tsin(n_lon, n_lat)
|
|
real(kind=4) :: T2min(n_lon, n_lat)
|
|
real(kind=4) :: q2min(n_lon, n_lat)
|
|
real(kind=4) :: u10min(n_lon, n_lat)
|
|
real(kind=4) :: v10min(n_lon, n_lat)
|
|
real(kind=4) :: PSFCin(n_lon, n_lat)
|
|
real(kind=4) :: seaicein(n_lon, n_lat)
|
|
!====================================================================================================================================
|
|
real(kind=4) :: Tout(nlevel)
|
|
real(kind=4) :: qout(nlevel)
|
|
real(kind=4) :: cloudout(6,nlevel-1)
|
|
real(kind=4) :: cfracout(nlevel-1)
|
|
real(kind=4) :: ozoneout(nlevel)
|
|
real(kind=4) :: tsout
|
|
real(kind=4) :: T2mout
|
|
real(kind=4) :: q2mout
|
|
real(kind=4) :: u10mout
|
|
real(kind=4) :: v10mout
|
|
real(kind=4) :: PSFCout
|
|
real(kind=4) :: seaiceout
|
|
|
|
|
|
|
|
Tout(:) = Tin(indexx0,indexy0,:)*fgeo1 + Tin(indexx0,indexy1,:)*fgeo2 &
|
|
+ Tin(indexx1,indexy1,:)*fgeo3 + Tin(indexx1,indexy0,:)*fgeo4
|
|
qout(:) = qin(indexx0,indexy0,:)*fgeo1 + qin(indexx0,indexy1,:)*fgeo2 &
|
|
+ qin(indexx1,indexy1,:)*fgeo3 + qin(indexx1,indexy0,:)*fgeo4
|
|
if(minval(Tin(indexx0:indexx1,indexy0:indexy1,:))<-999.)then
|
|
Tout(:) = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
if(minval(qin(indexx0:indexx1,indexy0:indexy1,:))<-999.)then
|
|
qout(:) = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
! print*,'%%%%%%%Tout=',Tout(:)
|
|
! print*,'Tin=',Tin(indexx0,indexy0,:)
|
|
cfracout(:) = cfracin(indexx0,indexy0,1:nlevel-1)*fgeo1 + cfracin(indexx0,indexy1,1:nlevel-1)*fgeo2 &
|
|
+ cfracin(indexx1,indexy1,1:nlevel-1)*fgeo3 + cfracin(indexx1,indexy0,1:nlevel-1)*fgeo4
|
|
if(minval(cfracin(indexx0:indexx1,indexy0:indexy1,:))<-999.)then
|
|
cfracout(:) = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
do j=1,6
|
|
cloudout(j,:) = cloudin(indexx0,indexy0,j,1:nlevel-1)*fgeo1 + cloudin(indexx0,indexy1,j,1:nlevel-1)*fgeo2 &
|
|
+ cloudin(indexx1,indexy1,j,1:nlevel-1)*fgeo3 + cloudin(indexx1,indexy0,j,1:nlevel-1)*fgeo4
|
|
enddo
|
|
if(minval(cloudin(indexx0:indexx1,indexy0:indexy1,:,:))<-999.)then
|
|
cloudout(:,:) = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
ozoneout(:) = ozonein(indexx0,indexy0,:)*fgeo1 + ozonein(indexx0,indexy1,:)*fgeo2 &
|
|
+ ozonein(indexx1,indexy1,:)*fgeo3 + ozonein(indexx1,indexy0,:)*fgeo4
|
|
if(minval(ozonein(indexx0:indexx1,indexy0:indexy1,:))<-999.)then
|
|
ozoneout(:) = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
tsout = tsin(indexx0,indexy0)*fgeo1 + tsin(indexx0,indexy1)*fgeo2 &
|
|
+ tsin(indexx1,indexy1)*fgeo3 + tsin(indexx1,indexy0)*fgeo4
|
|
if(minval(tsin(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
tsout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
T2mout = T2min(indexx0,indexy0)*fgeo1 + T2min(indexx0,indexy1)*fgeo2 &
|
|
+ T2min(indexx1,indexy1)*fgeo3 + T2min(indexx1,indexy0)*fgeo4
|
|
if(minval(T2min(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
T2mout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
q2mout = q2min(indexx0,indexy0)*fgeo1 + q2min(indexx0,indexy1)*fgeo2 &
|
|
+ q2min(indexx1,indexy1)*fgeo3 + q2min(indexx1,indexy0)*fgeo4
|
|
if(minval(q2min(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
q2mout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
u10mout = u10min(indexx0,indexy0)*fgeo1 + u10min(indexx0,indexy1)*fgeo2 &
|
|
+ u10min(indexx1,indexy1)*fgeo3 + u10min(indexx1,indexy0)*fgeo4
|
|
if(minval(u10min(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
u10mout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
v10mout = v10min(indexx0,indexy0)*fgeo1 + v10min(indexx0,indexy1)*fgeo2 &
|
|
+ v10min(indexx1,indexy1)*fgeo3 + v10min(indexx1,indexy0)*fgeo4
|
|
if(minval(v10min(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
v10mout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
PSFCout = PSFCin(indexx0,indexy0)*fgeo1 + PSFCin(indexx0,indexy1)*fgeo2 &
|
|
+ PSFCin(indexx1,indexy1)*fgeo3 + PSFCin(indexx1,indexy0)*fgeo4
|
|
if(minval(PSFCin(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
PSFCout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
seaiceout = seaicein(indexx0,indexy0)*fgeo1 + seaicein(indexx0,indexy1)*fgeo2 &
|
|
+ seaicein(indexx1,indexy1)*fgeo3 + seaicein(indexx1,indexy0)*fgeo4
|
|
if(minval(seaicein(indexx0:indexx1,indexy0:indexy1))<-999.)then
|
|
seaiceout = -9999.9
|
|
qcmark(ix,iy) = 1
|
|
endif
|
|
|
|
return
|
|
end subroutine select_obs_region
|
|
|
|
end module module_NWP_match_to_sat
|