78 logical fexist, opened
79 integer fsize, ncid, error, id_dim, nx, ny
80 character(len=256) :: outgrid =
"none" 81 character(len=256) :: inputorog =
"none" 82 integer :: mtnres,im,jm,nm,nr,nf0,nf1,efac,blat,nw
84 READ(5,*) mtnres,im,jm,nm,nr,nf0,nf1,efac,blat
98 print*,
"INPUTOROG=", trim(inputorog)
99 print*,
"IM,JM=", im, jm
107 print*, mtnres,im,jm,nm,nr,nf0,nf1,efac,blat
108 nw=(nm+1)*((nr+1)*nm+2)
111 print *,
' Starting terr12 mtnlm7_slm30.f IMN,JMN:',imn,jmn
114 if( trim(outgrid) .NE.
"none" )
then 115 inquire(file=trim(outgrid), exist=fexist)
116 if(.not. fexist)
then 117 print*,
"file "//trim(outgrid)//
" does not exist" 121 inquire( ncid,opened=opened )
122 if( .NOT.opened )
exit 125 print*,
"outgrid=", trim(outgrid)
126 error=nf__open(trim(outgrid),nf_nowrite,fsize,ncid)
127 call netcdf_err(error,
'Open file '//trim(outgrid) )
128 error=nf_inq_dimid(ncid,
'nx', id_dim)
129 call netcdf_err(error,
'inquire dimension nx from file '//
131 error=nf_inq_dimlen(ncid,id_dim,nx)
132 call netcdf_err(error,
'inquire dimension nx length '//
133 &
'from file '//trim(outgrid) )
135 error=nf_inq_dimid(ncid,
'ny', id_dim)
136 call netcdf_err(error,
'inquire dimension ny from file '//
138 error=nf_inq_dimlen(ncid,id_dim,ny)
139 call netcdf_err(error,
'inquire dimension ny length '//
140 &
'from file '//trim(outgrid) )
142 if(im .ne. nx/2)
then 143 print*,
"IM=",im,
" /= grid file nx/2=",nx/2
144 print*,
"Set IM = ", nx/2
147 if(jm .ne. ny/2)
then 148 print*,
"JM=",jm,
" /= grid file ny/2=",ny/2
149 print*,
"Set JM = ", ny/2
153 call netcdf_err(error,
'close file '//trim(outgrid) )
157 CALL tersub(imn,jmn,im,jm,nm,nr,nf0,nf1,nw,efac,blat,
182 SUBROUTINE tersub(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
187 integer :: IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW
188 character(len=*),
intent(in) :: OUTGRID
189 character(len=*),
intent(in) :: INPUTOROG
191 real,
parameter :: MISSING_VALUE=-9999.
192 real,
PARAMETER :: PI=3.1415926535897931
193 integer,
PARAMETER :: NMT=14
195 integer :: efac,blat,zsave1,zsave2,itopo,kount
196 integer :: kount2,islmx,jslmx,oldslm,msksrc,mskocn,notocn
197 integer :: i,j,nx,ny,ncid,js,jn,iw,ie,k,it,jt,i1,error,id_dim
198 integer :: id_var,nx_in,ny_in,fsize,wgta,IN,INW,INE,IS,ISW,ISE
199 integer :: M,N,IMT,IRET,ios,iosg,latg2,istat,itest,jtest
200 integer :: i_south_pole,j_south_pole,i_north_pole,j_north_pole
201 integer :: maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
205 integer,
allocatable :: JST(:),JEN(:),KPDS(:),KGDS(:),numi(:)
206 integer,
allocatable :: lonsperlat(:)
208 integer,
allocatable :: IST(:,:),IEN(:,:),ZSLMX(:,:)
209 integer,
allocatable :: ZAVG(:,:),ZSLM(:,:)
210 integer(1),
allocatable :: UMD(:,:)
211 integer(2),
allocatable :: glob(:,:)
213 integer,
allocatable :: IWORK(:,:,:)
215 real :: DEGRAD,maxlat, minlat,timef,tbeg,tend,tbeg1
216 real :: PHI,DELXN,RS,RN,slma,oroa,vara,var4a,xn,XS,FFF,WWW
217 real :: sumdif,avedif
219 real,
allocatable :: COSCLT(:),WGTCLT(:),RCLT(:),XLAT(:),DIFFX(:)
220 real,
allocatable :: XLON(:),ORS(:),oaa(:),ola(:),GLAT(:)
222 real,
allocatable :: GEOLON(:,:),GEOLON_C(:,:),DX(:,:)
223 real,
allocatable :: GEOLAT(:,:),GEOLAT_C(:,:),DY(:,:)
224 real,
allocatable :: SLM(:,:),ORO(:,:),VAR(:,:),ORF(:,:)
225 real,
allocatable :: land_frac(:,:)
226 real,
allocatable :: THETA(:,:),GAMMA(:,:),SIGMA(:,:),ELVMAX(:,:)
227 real,
allocatable :: VAR4(:,:),SLMI(:,:)
228 real,
allocatable :: WORK1(:,:),WORK2(:,:),WORK3(:,:),WORK4(:,:)
229 real,
allocatable :: WORK5(:,:),WORK6(:,:)
230 real,
allocatable :: tmpvar(:,:)
231 real,
allocatable :: slm_in(:,:), lon_in(:,:), lat_in(:,:)
232 real(4),
allocatable:: GICE(:,:),OCLSM(:,:)
234 real,
allocatable :: OA(:,:,:),OL(:,:,:),HPRIME(:,:,:)
235 real,
allocatable :: oa_in(:,:,:), ol_in(:,:,:)
238 complex :: ffj(im/2+1)
240 logical :: grid_from_file,output_binary,fexist,opened
241 logical :: SPECTR, REVLAT, FILTER
243 logical :: is_south_pole(IM,JM), is_north_pole(IM,JM)
246 output_binary = .false.
251 allocate (jst(jm),jen(jm),kpds(200),kgds(200),numi(jm))
252 allocate (lonsperlat(jm/2))
253 allocate (ist(im,jm),ien(im,jm),zslmx(2700,1350))
254 allocate (glob(imn,jmn))
257 allocate (cosclt(jm),wgtclt(jm),rclt(jm),xlat(jm),diffx(jm/2))
258 allocate (xlon(im),ors(nw),oaa(4),ola(4),glat(jmn))
260 allocate (zavg(imn,jmn))
261 allocate (zslm(imn,jmn))
262 allocate (umd(imn,jmn))
280 print *,
' In TERSUB, ITOPO=',itopo
281 if (mskocn .eq. 1)
then 282 print *,
' Ocean Model LSM Present and ' 283 print *,
' Overrides OCEAN POINTS in LSM: mskocn=',mskocn
284 if (notocn .eq. 1)
then 285 print *,
' Ocean LSM Reversed: NOTOCN=',notocn
303 IF (msksrc .eq. 0 )
then 304 READ(14,12,iostat=ios) zslmx
308 print *,
' navy10 lake mask rd fail -- ios,MSKSRC:',ios,msksrc
311 print *,
' Attempt to open/read UMD 30" slmsk MSKSRC=',msksrc
317 open(10,file=
"landcover30.fixed",
318 & recl=43200*21600, access=
'direct',iostat=istat)
322 print *,
' UMD lake mask open failed -- ios,MSKSRC:',istat,msksrc
325 read(10, rec=1,iostat=istat) umd
326 print *,
' UMD lake mask opened OK -- ios,MSKSRC:',istat,msksrc
332 print *,
' UMD lake mask rd err -- trying navy 10',istat
334 print *,
' ***** MSKSRC set to 0 MSKSRC=',msksrc
335 if (msksrc .eq. 0 )
then 336 READ(14,12,iostat=ios) zslmx
339 print *,
' navy10 lake mask rd fail - ios,MSKSRC:',ios,msksrc
343 print *,
' UMD lake, UMD(50,50)=',umd(50,50),msksrc
351 print *,
' About to call read_g, ITOPO=',itopo
352 if ( itopo .ne. 0 )
call read_g(glob,itopo)
371 print *,
' After read_g, glob(500,500)=',glob(500,500)
375 print*,
' IM, JM, NM, NR, NF0, NF1, EFAC, BLAT' 376 print*, im,jm,nm,nr,nf0,nf1,efac,blat
377 print *,
' imn,jmn,glob(imn,jmn)=',imn,jmn,glob(imn,jmn)
378 print *,
' UBOUND ZAVG=',ubound(zavg)
379 print *,
' UBOUND glob=',ubound(glob)
380 print *,
' UBOUND ZSLM=',ubound(zslm)
381 print *,
' UBOUND GICE=',imn+1,3601
382 print *,
' UBOUND OCLSM=',im,jm
416 if ( umd(i,j) .eq. 0 ) zslm(i,j) = 0
424 print *,
' NAVY 10 (8) slmsk for lakes, MSKSRC=',msksrc
432 if ( glob(i,j) .eq. -9999 )
then 438 if ( zslmx(islmx,jslmx) .eq. 0 )
then 439 if ( j .gt. 8 .and. j .lt. jmn-8 )
then 440 if (i1 .gt. imn ) i1 = i1 - imn
442 if(zslm(i,j).eq.1 .and. oldslm .eq. 1 .and. zslm(i1,j).eq.1)
then 443 if (i .ne. 1) oldslm = zslm(i,j)
454 print *,
' ***** set slm from 30" glob, MSKSRC=',msksrc
461 if ( glob(i,j) .eq. -9999 )
then 469 deallocate (zslmx,umd,glob)
484 read(20,*,iostat=ios) latg2,lonsperlat
485 if(ios.ne.0.or.2*latg2.ne.jm)
then 489 print *,ios,latg2,
'COMPUTE TERRAIN ON A FULL GAUSSIAN GRID' 492 numi(j)=lonsperlat(j)
495 numi(j)=lonsperlat(jm+1-j)
497 print *,ios,latg2,
'COMPUTE TERRAIN ON A REDUCED GAUSSIAN GRID',
507 print *,
' SPECTR=',spectr,
' REVLAT=',revlat,
' ** with GICE-07 **' 509 CALL splat(4,jm,cosclt,wgtclt)
511 rclt(j) = acos(cosclt(j))
514 phi = rclt(j) * degrad
516 xlat(jm-j+1) = phi - 90.
519 CALL splat(0,jm,cosclt,wgtclt)
521 rclt(j) = acos(cosclt(j))
522 xlat(j) = 90.0 - rclt(j) * degrad
526 allocate (gice(imn+1,3601))
530 diffx(j) = xlat(j) - xlat(j-1)
531 sumdif = sumdif + diffx(j)
533 avedif=sumdif/(float(jm/2))
537 write (6,106) (xlat(j),j=jm,1,-1)
538 106
format( 10(f7.3,1x))
539 107
format( 10(f9.5,1x))
544 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
547 &
' Before GICE ZAVG(1,2)=',zavg(1,2),zslm(1,2)
549 &
' Before GICE ZAVG(1,12)=',zavg(1,12),zslm(1,12)
551 &
' Before GICE ZAVG(1,52)=',zavg(1,52),zslm(1,52)
553 &
' Before GICE ZAVG(1,112)=',zavg(1,jmn-112),zslm(1,112)
557 READ(15,iostat=iosg) gice
558 if(iosg .ne. 0 )
then 559 print *,
' *** Err on reading GICE record, iosg=',iosg
560 print *,
' exec continues but NO GICE correction done ' 563 print *,
' GICE 30" Antarctica RAMP orog 43200x3616 read OK' 564 print *,
' Processing! ' 565 print *,
' Processing! ' 566 print *,
' Processing! ' 571 if( gice(i,j) .ne. -99. .and. gice(i,j) .ne. -1.0 )
then 572 if ( gice(i,j) .gt. 0.)
then 573 zavg(i,j) = int( gice(i,j) + 0.5 )
579 152
format(1x,
' ZAVG(i=',i4,
' j=',i4,
')=',i5,i3,
580 &
' orig:',i5,i4,
' Lat=',f7.3,f8.2,
'E',
' GICE=',f8.1)
587 allocate (oclsm(im,jm),slmi(im,jm))
594 if (mskocn .eq. 1)
then 598 open(25,form=
'formatted',iostat=istat)
601 print *,
' Ocean lsm file Open failure: mskocn,istat=',mskocn,istat
604 print *,
' Ocean lsm file Opened OK: mskocn,istat=',mskocn,istat
610 read(25,*,iostat=ios)oclsm
615 print *,
' Rd fail: Ocean lsm - continue, mskocn,ios=',mskocn,ios
618 print *,
' Rd OK: ocean lsm: mskocn,ios=',mskocn,ios
622 if ( mskocn .eq. 1 )
then 625 if ( notocn .eq. 0 )
then 626 slmi(i,j) = float(nint(oclsm(i,j)))
628 if ( nint(oclsm(i,j)) .eq. 0)
then 636 print *,
' OCLSM',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
637 print *,
' SLMI:',slmi(1,1),slmi(50,50),slmi(75,75),slmi(im,jm)
645 print *,
' Not using Ocean model land sea mask' 648 if (mskocn .eq. 1)
then 649 print *,
' LSM:',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
652 allocate (geolon(im,jm),geolon_c(im+1,jm+1),dx(im,jm))
653 allocate (geolat(im,jm),geolat_c(im+1,jm+1),dy(im,jm))
654 allocate (slm(im,jm),oro(im,jm),var(im,jm),var4(im,jm))
655 allocate (land_frac(im,jm))
658 grid_from_file = .false.
659 is_south_pole = .false.
660 is_north_pole = .false.
665 if( trim(outgrid) .NE.
"none" )
then 666 grid_from_file = .true.
667 inquire(file=trim(outgrid), exist=fexist)
668 if(.not. fexist)
then 669 print*,
"file "//trim(outgrid)//
" does not exist" 673 inquire( ncid,opened=opened )
674 if( .NOT.opened )
exit 677 print*,
"outgrid=", trim(outgrid)
678 error=nf__open(trim(outgrid),nf_nowrite,fsize,ncid)
679 call netcdf_err(error,
'Open file '//trim(outgrid) )
680 error=nf_inq_dimid(ncid,
'nx', id_dim)
681 call netcdf_err(error,
'inquire dimension nx from file '//
704 print*,
"Read the grid from file "//trim(outgrid)
706 allocate(tmpvar(nx+1,ny+1))
708 error=nf_inq_varid(ncid,
'x', id_var)
709 call netcdf_err(error,
'inquire varid of x from file ' 711 error=nf_get_var_double(ncid, id_var, tmpvar)
712 call netcdf_err(error,
'inquire data of x from file ' 715 do j = 1,ny+1;
do i = 1,nx+1
716 if(tmpvar(i,j) .NE. missing_value)
then 717 if(tmpvar(i,j) .GT. 360) tmpvar(i,j) = tmpvar(i,j) - 360
718 if(tmpvar(i,j) .LT. 0) tmpvar(i,j) = tmpvar(i,j) + 360
722 geolon(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
723 geolon_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
725 error=nf_inq_varid(ncid,
'y', id_var)
726 call netcdf_err(error,
'inquire varid of y from file ' 728 error=nf_get_var_double(ncid, id_var, tmpvar)
729 call netcdf_err(error,
'inquire data of y from file ' 731 geolat(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
732 geolat_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
741 do j = 1, ny+1;
do i = 1, nx+1
742 if( tmpvar(i,j) > maxlat )
then 747 if( tmpvar(i,j) < minlat )
then 754 if(maxlat < 89.9 )
then 758 if(minlat > -89.9 )
then 762 print*,
"minlat=", minlat,
"maxlat=", maxlat
763 print*,
"north pole supergrid index is ",
764 & i_north_pole, j_north_pole
765 print*,
"south pole supergrid index is ",
766 & i_south_pole, j_south_pole
769 if(i_south_pole >0 .and. j_south_pole > 0)
then 770 if(mod(i_south_pole,2)==0)
then 771 do j = 1, jm;
do i = 1, im
772 if(i==i_south_pole/2 .and. (j==j_south_pole/2
773 & .or. j==j_south_pole/2+1) )
then 774 is_south_pole(i,j) = .true.
775 print*,
"south pole at i,j=", i, j
779 do j = 1, jm;
do i = 1, im
780 if((i==i_south_pole/2 .or. i==i_south_pole/2+1)
781 & .and. (j==j_south_pole/2 .or.
782 & j==j_south_pole/2+1) )
then 783 is_south_pole(i,j) = .true.
784 print*,
"south pole at i,j=", i, j
789 if(i_north_pole >0 .and. j_north_pole > 0)
then 790 if(mod(i_north_pole,2)==0)
then 791 do j = 1, jm;
do i = 1, im
792 if(i==i_north_pole/2 .and. (j==j_north_pole/2 .or.
793 & j==j_north_pole/2+1) )
then 794 is_north_pole(i,j) = .true.
795 print*,
"north pole at i,j=", i, j
799 do j = 1, jm;
do i = 1, im
800 if((i==i_north_pole/2 .or. i==i_north_pole/2+1)
801 & .and. (j==j_north_pole/2 .or.
802 & j==j_north_pole/2+1) )
then 803 is_north_pole(i,j) = .true.
804 print*,
"north pole at i,j=", i, j
811 allocate(tmpvar(nx,ny))
812 error=nf_inq_varid(ncid,
'area', id_var)
813 call netcdf_err(error,
'inquire varid of area from file ' 815 error=nf_get_var_double(ncid, id_var, tmpvar)
816 call netcdf_err(error,
'inquire data of area from file ' 821 dx(i,j) = sqrt(tmpvar(2*i-1,2*j-1)+tmpvar(2*i,2*j-1)
822 & +tmpvar(2*i-1,2*j )+tmpvar(2*i,2*j ))
848 write(6,*)
' Timer 1 time= ',tend-tbeg
850 if(grid_from_file)
then 852 CALL makemt2(zavg,zslm,oro,slm,land_frac,var,var4,glat,
853 & im,jm,imn,jmn,geolon_c,geolat_c)
855 write(6,*)
' MAKEMT2 time= ',tend-tbeg
857 CALL makemt(zavg,zslm,oro,slm,var,var4,glat,
858 & ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
861 call minmxj(im,jm,oro,
' ORO')
862 call minmxj(im,jm,slm,
' SLM')
863 call minmxj(im,jm,var,
' VAR')
864 call minmxj(im,jm,var4,
' VAR4')
885 allocate (theta(im,jm),gamma(im,jm),sigma(im,jm),elvmax(im,jm))
886 if(grid_from_file)
then 888 CALL makepc2(zavg,zslm,theta,gamma,sigma,glat,
889 1 im,jm,imn,jmn,geolon_c,geolat_c)
891 write(6,*)
' MAKEPC2 time= ',tend-tbeg
893 CALL makepc(zavg,zslm,theta,gamma,sigma,glat,
894 1 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
897 call minmxj(im,jm,theta,
' THETA')
898 call minmxj(im,jm,gamma,
' GAMMA')
899 call minmxj(im,jm,sigma,
' SIGMA')
908 allocate (iwork(im,jm,4))
909 allocate (oa(im,jm,4),ol(im,jm,4),hprime(im,jm,14))
910 allocate (work1(im,jm),work2(im,jm),work3(im,jm),work4(im,jm))
911 allocate (work5(im,jm),work6(im,jm))
913 call minmxj(im,jm,oro,
' ORO')
914 print*,
"inputorog=", trim(inputorog)
915 if(grid_from_file)
then 916 if(trim(inputorog) ==
"none")
then 917 print*,
"calling MAKEOA2 to compute OA, OL" 919 CALL makeoa2(zavg,zslm,var,glat,oa,ol,iwork,elvmax,oro,
920 1 work1,work2,work3,work4,work5,work6,
921 2 im,jm,imn,jmn,geolon_c,geolat_c,
922 3 geolon,geolat,dx,dy,is_south_pole,is_north_pole)
924 write(6,*)
' MAKEOA2 time= ',tend-tbeg
927 error=nf__open(trim(inputorog),nf_nowrite,fsize,ncid)
928 call netcdf_err(error,
'Open file '//trim(inputorog) )
929 error=nf_inq_dimid(ncid,
'lon', id_dim)
930 call netcdf_err(error,
'inquire dimension lon from file '//
932 error=nf_inq_dimlen(ncid,id_dim,nx_in)
933 call netcdf_err(error,
'inquire dimension lon length '//
934 &
'from file '//trim(inputorog) )
935 error=nf_inq_dimid(ncid,
'lat', id_dim)
936 call netcdf_err(error,
'inquire dimension lat from file '//
938 error=nf_inq_dimlen(ncid,id_dim,ny_in)
939 call netcdf_err(error,
'inquire dimension lat length '//
940 &
'from file '//trim(inputorog) )
942 print*,
"extrapolate OA, OL from Gaussian grid with nx=",
943 & nx_in,
", ny=", ny_in
944 allocate(oa_in(nx_in,ny_in,4), ol_in(nx_in,ny_in,4))
945 allocate(slm_in(nx_in,ny_in) )
946 allocate(lon_in(nx_in,ny_in), lat_in(nx_in,ny_in) )
948 error=nf_inq_varid(ncid,
'oa1', id_var)
949 call netcdf_err(error,
'inquire varid of oa1 from file ' 950 & //trim(inputorog) )
951 error=nf_get_var_double(ncid, id_var, oa_in(:,:,1))
952 call netcdf_err(error,
'inquire data of oa1 from file ' 953 & //trim(inputorog) )
954 error=nf_inq_varid(ncid,
'oa2', id_var)
955 call netcdf_err(error,
'inquire varid of oa2 from file ' 956 & //trim(inputorog) )
957 error=nf_get_var_double(ncid, id_var, oa_in(:,:,2))
958 call netcdf_err(error,
'inquire data of oa2 from file ' 959 & //trim(inputorog) )
960 error=nf_inq_varid(ncid,
'oa3', id_var)
961 call netcdf_err(error,
'inquire varid of oa3 from file ' 962 & //trim(inputorog) )
963 error=nf_get_var_double(ncid, id_var, oa_in(:,:,3))
964 call netcdf_err(error,
'inquire data of oa3 from file ' 965 & //trim(inputorog) )
966 error=nf_inq_varid(ncid,
'oa4', id_var)
967 call netcdf_err(error,
'inquire varid of oa4 from file ' 968 & //trim(inputorog) )
969 error=nf_get_var_double(ncid, id_var, oa_in(:,:,4))
970 call netcdf_err(error,
'inquire data of oa4 from file ' 971 & //trim(inputorog) )
973 error=nf_inq_varid(ncid,
'ol1', id_var)
974 call netcdf_err(error,
'inquire varid of ol1 from file ' 975 & //trim(inputorog) )
976 error=nf_get_var_double(ncid, id_var, ol_in(:,:,1))
977 call netcdf_err(error,
'inquire data of ol1 from file ' 978 & //trim(inputorog) )
979 error=nf_inq_varid(ncid,
'ol2', id_var)
980 call netcdf_err(error,
'inquire varid of ol2 from file ' 981 & //trim(inputorog) )
982 error=nf_get_var_double(ncid, id_var, ol_in(:,:,2))
983 call netcdf_err(error,
'inquire data of ol2 from file ' 984 & //trim(inputorog) )
985 error=nf_inq_varid(ncid,
'ol3', id_var)
986 call netcdf_err(error,
'inquire varid of ol3 from file ' 987 & //trim(inputorog) )
988 error=nf_get_var_double(ncid, id_var, ol_in(:,:,3))
989 call netcdf_err(error,
'inquire data of ol3 from file ' 990 & //trim(inputorog) )
991 error=nf_inq_varid(ncid,
'ol4', id_var)
992 call netcdf_err(error,
'inquire varid of ol4 from file ' 993 & //trim(inputorog) )
994 error=nf_get_var_double(ncid, id_var, ol_in(:,:,4))
995 call netcdf_err(error,
'inquire data of ol4 from file ' 996 & //trim(inputorog) )
998 error=nf_inq_varid(ncid,
'slmsk', id_var)
999 call netcdf_err(error,
'inquire varid of slmsk from file ' 1000 & //trim(inputorog) )
1001 error=nf_get_var_double(ncid, id_var, slm_in)
1002 call netcdf_err(error,
'inquire data of slmsk from file ' 1003 & //trim(inputorog) )
1005 error=nf_inq_varid(ncid,
'geolon', id_var)
1006 call netcdf_err(error,
'inquire varid of geolon from file ' 1007 & //trim(inputorog) )
1008 error=nf_get_var_double(ncid, id_var, lon_in)
1009 call netcdf_err(error,
'inquire data of geolon from file ' 1010 & //trim(inputorog) )
1012 error=nf_inq_varid(ncid,
'geolat', id_var)
1013 call netcdf_err(error,
'inquire varid of geolat from file ' 1014 & //trim(inputorog) )
1015 error=nf_get_var_double(ncid, id_var, lat_in)
1016 call netcdf_err(error,
'inquire data of geolat from file ' 1017 & //trim(inputorog) )
1020 do j=1,ny_in;
do i=1,nx_in
1021 if(slm_in(i,j) == 2) slm_in(i,j) = 0
1024 error=nf_close(ncid)
1026 & //trim(inputorog) )
1028 print*,
"calling MAKEOA3 to compute OA, OL" 1029 CALL makeoa3(zavg,zslm,var,glat,oa,ol,iwork,elvmax,oro,slm,
1030 1 work1,work2,work3,work4,work5,work6,
1031 2 im,jm,imn,jmn,geolon_c,geolat_c,
1032 3 geolon,geolat,is_south_pole,is_north_pole,nx_in,ny_in,
1033 4 oa_in,ol_in,slm_in,lon_in,lat_in)
1035 deallocate(oa_in,ol_in,slm_in,lon_in,lat_in)
1039 CALL makeoa(zavg,var,glat,oa,ol,iwork,elvmax,oro,
1040 1 work1,work2,work3,work4,
1042 3 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
1047 deallocate (zslm,zavg)
1049 deallocate (work2,work3,work4,work5,work6)
1055 call minmxj(im,jm,oa,
' OA')
1056 call minmxj(im,jm,ol,
' OL')
1057 call minmxj(im,jm,elvmax,
' ELVMAX')
1058 call minmxj(im,jm,oro,
' ORO')
1078 if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
1079 if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
1080 if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
1081 if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
1082 if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
1083 if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
1086 print *,
' MAXC3:',maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
1091 print *,
' ===> Replacing ELVMAX with ELVMAX-ORO <=== ' 1092 print *,
' ===> if ELVMAX<=ORO replace with proxy <=== ' 1093 print *,
' ===> the sum of mean orog (ORO) and std dev <=== ' 1096 if (elvmax(i,j) .lt. oro(i,j) )
then 1098 elvmax(i,j) = max( 3. * var(i,j),0.)
1100 elvmax(i,j) = max( elvmax(i,j) - oro(i,j),0.)
1112 if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
1113 if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
1114 if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
1115 if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
1116 if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
1117 if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
1120 print *,
' after MAXC 3-6 km:',maxc3,maxc4,maxc5,maxc6
1122 call mnmxja(im,jm,elvmax,itest,jtest,
' ELVMAX')
1127 print *,
' Testing at point (itest,jtest)=',itest,jtest
1128 print *,
' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
1129 print *,
' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
1132 IF(slm(i,j).EQ.0.)
THEN 1161 if ( mskocn .eq. 1 )
then 1165 if (abs(oro(i,j)) .lt. 1. )
then 1166 slm(i,j) = slmi(i,j)
1168 if ( slmi(i,j) .eq. 1. .and. slm(i,j) .eq. 1) slm(i,j) = 1
1169 if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
1170 if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 1) slm(i,j) = 0
1171 if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
1176 print *,
' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
1177 print *,
' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
1185 rn=
REAL(NUMI(JN))/
REAL(NUMI(J))
1186 rs=
REAL(NUMI(JS))/
REAL(NUMI(J))
1190 slma=slm(iw,j)+slm(ie,j)
1191 oroa=oro(iw,j)+oro(ie,j)
1192 vara=var(iw,j)+var(ie,j)
1193 var4a=var4(iw,j)+var4(ie,j)
1195 oaa(k)=oa(iw,j,k)+oa(ie,j,k)
1197 ola(k)=ol(iw,j,k)+ol(ie,j,k)
1201 IF(abs(xn-nint(xn)).LT.1.e-2)
THEN 1202 in=mod(nint(xn)-1,numi(jn))+1
1203 inw=mod(in+numi(jn)-2,numi(jn))+1
1204 ine=mod(in,numi(jn))+1
1205 slma=slma+slm(inw,jn)+slm(in,jn)+slm(ine,jn)
1206 oroa=oroa+oro(inw,jn)+oro(in,jn)+oro(ine,jn)
1207 vara=vara+var(inw,jn)+var(in,jn)+var(ine,jn)
1208 var4a=var4a+var4(inw,jn)+var4(in,jn)+var4(ine,jn)
1210 oaa(k)=oaa(k)+oa(inw,jn,k)+oa(in,jn,k)+oa(ine,jn,k)
1211 ola(k)=ola(k)+ol(inw,jn,k)+ol(in,jn,k)+ol(ine,jn,k)
1216 ine=mod(inw,numi(jn))+1
1217 slma=slma+slm(inw,jn)+slm(ine,jn)
1218 oroa=oroa+oro(inw,jn)+oro(ine,jn)
1219 vara=vara+var(inw,jn)+var(ine,jn)
1220 var4a=var4a+var4(inw,jn)+var4(ine,jn)
1222 oaa(k)=oaa(k)+oa(inw,jn,k)+oa(ine,jn,k)
1223 ola(k)=ola(k)+ol(inw,jn,k)+ol(ine,jn,k)
1228 IF(abs(xs-nint(xs)).LT.1.e-2)
THEN 1229 is=mod(nint(xs)-1,numi(js))+1
1230 isw=mod(is+numi(js)-2,numi(js))+1
1231 ise=mod(is,numi(js))+1
1232 slma=slma+slm(isw,js)+slm(is,js)+slm(ise,js)
1233 oroa=oroa+oro(isw,js)+oro(is,js)+oro(ise,js)
1234 vara=vara+var(isw,js)+var(is,js)+var(ise,js)
1235 var4a=var4a+var4(isw,js)+var4(is,js)+var4(ise,js)
1237 oaa(k)=oaa(k)+oa(isw,js,k)+oa(is,js,k)+oa(ise,js,k)
1238 ola(k)=ola(k)+ol(isw,js,k)+ol(is,js,k)+ol(ise,js,k)
1243 ise=mod(isw,numi(js))+1
1244 slma=slma+slm(isw,js)+slm(ise,js)
1245 oroa=oroa+oro(isw,js)+oro(ise,js)
1246 vara=vara+var(isw,js)+var(ise,js)
1247 var4a=var4a+var4(isw,js)+var4(ise,js)
1249 oaa(k)=oaa(k)+oa(isw,js,k)+oa(ise,js,k)
1250 ola(k)=ola(k)+ol(isw,js,k)+ol(ise,js,k)
1261 IF(slm(i,j).EQ.0..AND.slma.EQ.wgta)
THEN 1262 print
'("SEA ",2F8.0," MODIFIED TO LAND",2F8.0," AT ",2I8)',
1263 & oro(i,j),var(i,j),oroa,vara,i,j
1272 ELSEIF(slm(i,j).EQ.1..AND.slma.EQ.0.)
THEN 1273 print
'("LAND",2F8.0," MODIFIED TO SEA ",2F8.0," AT ",2I8)',
1274 & oro(i,j),var(i,j),oroa,vara,i,j
1287 print *,
' after isolated points removed' 1288 call minmxj(im,jm,oro,
' ORO')
1290 print *,
' ORO(itest,jtest)=',oro(itest,jtest)
1291 print *,
' VAR(itest,jtest)=',var(itest,jtest)
1292 print *,
' VAR4(itest,jtest)=',var4(itest,jtest)
1293 print *,
' OA(itest,jtest,1)=',oa(itest,jtest,1)
1294 print *,
' OA(itest,jtest,2)=',oa(itest,jtest,2)
1295 print *,
' OA(itest,jtest,3)=',oa(itest,jtest,3)
1296 print *,
' OA(itest,jtest,4)=',oa(itest,jtest,4)
1297 print *,
' OL(itest,jtest,1)=',ol(itest,jtest,1)
1298 print *,
' OL(itest,jtest,2)=',ol(itest,jtest,2)
1299 print *,
' OL(itest,jtest,3)=',ol(itest,jtest,3)
1300 print *,
' OL(itest,jtest,4)=',ol(itest,jtest,4)
1301 print *,
' Testing at point (itest,jtest)=',itest,jtest
1302 print *,
' THETA(itest,jtest)=',theta(itest,jtest)
1303 print *,
' GAMMA(itest,jtest)=',gamma(itest,jtest)
1304 print *,
' SIGMA(itest,jtest)=',sigma(itest,jtest)
1305 print *,
' ELVMAX(itest,jtest)=',elvmax(itest,jtest)
1306 print *,
' EFAC=',efac
1310 oro(i,j) = oro(i,j) + efac*var(i,j)
1311 hprime(i,j,1) = var(i,j)
1312 hprime(i,j,2) = var4(i,j)
1313 hprime(i,j,3) = oa(i,j,1)
1314 hprime(i,j,4) = oa(i,j,2)
1315 hprime(i,j,5) = oa(i,j,3)
1316 hprime(i,j,6) = oa(i,j,4)
1317 hprime(i,j,7) = ol(i,j,1)
1318 hprime(i,j,8) = ol(i,j,2)
1319 hprime(i,j,9) = ol(i,j,3)
1320 hprime(i,j,10)= ol(i,j,4)
1321 hprime(i,j,11)= theta(i,j)
1322 hprime(i,j,12)= gamma(i,j)
1323 hprime(i,j,13)= sigma(i,j)
1324 hprime(i,j,14)= elvmax(i,j)
1328 call mnmxja(im,jm,elvmax,itest,jtest,
' ELVMAX')
1337 allocate (orf(im,jm))
1338 IF ( nf1 - nf0 .eq. 0 ) filter=.false.
1339 print *,
' NF1, NF0, FILTER=',nf1,nf0,filter
1343 if(numi(j).lt.im)
then 1345 call spfft1(numi(j),im/2+1,numi(j),1,ffj,oro(1,j),-1)
1346 call spfft1(im,im/2+1,im,1,ffj,oro(1,j),+1)
1349 CALL sptez(nr,nm,4,im,jm,ors,oro,-1)
1351 print *,
' about to apply spectral filter ' 1357 www=max(1.-fff*(n-nf0)**2,0.)
1358 ors(i+1)=ors(i+1)*www
1359 ors(i+2)=ors(i+2)*www
1365 CALL sptez(nr,nm,4,im,jm,ors,orf,+1)
1367 if(numi(j).lt.im)
then 1368 call spfft1(im,im/2+1,im,1,ffj,orf(1,j),-1)
1369 call spfft1(numi(j),im/2+1,numi(j),1,ffj,orf(1,j),+1)
1375 CALL revers(im, jm, numi, slm, work1)
1376 CALL revers(im, jm, numi, oro, work1)
1378 CALL revers(im, jm, numi, hprime(1,1,imt), work1)
1387 call mnmxja(im,jm,elvmax,itest,jtest,
' ELVMAX')
1388 print *,
' ELVMAX(',itest,jtest,
')=',elvmax(itest,jtest)
1389 print *,
' after spectral filter is applied' 1390 call minmxj(im,jm,oro,
' ORO')
1391 call minmxj(im,jm,orf,
' ORF')
1394 call rg2gg(im,jm,numi,slm)
1395 call rg2gg(im,jm,numi,oro)
1396 call rg2gg(im,jm,numi,orf)
1399 call rg2gg(im,jm,numi,hprime(1,1,imt))
1402 print *,
' after nearest neighbor interpolation applied ' 1403 call minmxj(im,jm,oro,
' ORO')
1404 call minmxj(im,jm,orf,
' ORF')
1405 call mnmxja(im,jm,elvmax,itest,jtest,
' ELVMAX')
1406 print *,
' ORO,ORF(itest,jtest),itest,jtest:',
1407 & oro(itest,jtest),orf(itest,jtest),itest,jtest
1408 print *,
' ELVMAX(',itest,jtest,
')=',elvmax(itest,jtest)
1414 if ( i .le. 21 .and. i .ge. 1 )
then 1415 if (j .eq. jm )
write(6,153)i,j,oro(i,j),elvmax(i,j),slm(i,j)
1416 153
format(1x,
' ORO,ELVMAX(i=',i4,
' j=',i4,
')=',2e14.5,f5.1)
1421 write(6,*)
' Timer 5 time= ',tend-tbeg
1422 if (output_binary)
then 1425 print *,
' OUTPUT BINARY FIELDS' 1426 WRITE(51)
REAL(SLM,4)
1427 WRITE(52)
REAL(ORF,4)
1428 WRITE(53)
REAL(HPRIME,4)
1429 WRITE(54)
REAL(ORS,4)
1430 WRITE(55)
REAL(ORO,4)
1431 WRITE(66)
REAL(THETA,4)
1432 WRITE(67)
REAL(GAMMA,4)
1433 WRITE(68)
REAL(SIGMA,4)
1435 if ( mskocn .eq. 1 )
then 1437 WRITE(27,iostat=ios) oclsm
1438 print *,
' write OCLSM input:',ios
1442 call minmxj(im,jm,oro,
' ORO')
1443 print *,
' IM=',im,
' JM=',jm,
' SPECTR=',spectr
1445 WRITE(71)
REAL(SLM,4)
1447 WRITE(71)
REAL(HPRIME(:,:,IMT),4)
1448 print *,
' HPRIME(',itest,jtest,imt,
')=',hprime(itest,jtest,imt)
1450 WRITE(71)
REAL(ORO,4)
1452 WRITE(71)
REAL(ORF,4) 1477 kgds(4)=90000-180000/pi*rclt(1)
1479 kgds(7)=180000/pi*rclt(1)-90000
1480 kgds(8)=-nint(360000./im)
1481 kgds(9)=nint(360000./im)
1485 CALL baopen(56,
'fort.56',iret)
1486 if (iret .ne. 0) print *,
' BAOPEN ERROR UNIT 56: IRET=',iret
1487 CALL putgb(56,im*jm,kpds,kgds,lb,slm,iret)
1488 print *,
' SLM: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1489 if (iret .ne. 0) print *,
' SLM PUTGB ERROR: UNIT 56: IRET=',iret
1490 print *,
' SLM: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1502 CALL baopen(57,
'fort.57',iret)
1503 CALL putgb(57,im*jm,kpds,kgds,lb,orf,iret)
1504 print *,
' ORF (ORO): putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1511 CALL baopen(58,
'fort.58',iret)
1512 CALL putgb(58,im*jm,kpds,kgds,lb,theta,iret)
1513 print *,
' THETA: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1520 CALL baopen(60,
'fort.60',iret)
1521 CALL putgb(60,im*jm,kpds,kgds,lb,sigma,iret)
1522 print *,
' SIGMA: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1529 CALL baopen(59,
'fort.59',iret)
1530 CALL putgb(59,im*jm,kpds,kgds,lb,gamma,iret)
1531 print *,
' GAMMA: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1535 CALL baopen(61,
'fort.61',iret)
1536 CALL putgb(61,im*jm,kpds,kgds,lb,hprime,iret)
1537 print *,
' HPRIME: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1542 CALL baopen(62,
'fort.62',iret)
1543 CALL putgb(62,im*jm,kpds,kgds,lb,elvmax,iret)
1544 print *,
' ELVMAX: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1549 xlon(i) = delxn*(i-1)
1551 IF(trim(outgrid) ==
"none")
THEN 1554 geolon(i,j) = xlon(i)
1555 geolat(i,j) = xlat(j)
1560 xlat(j) = geolat(1,j)
1563 xlon(i) = geolon(i,1)
1567 write(6,*)
' Binary output time= ',tend-tbeg
1569 CALL write_netcdf(im,jm,slm,land_frac,oro,orf,hprime,1,1,
1570 1 geolon(1:im,1:jm),geolat(1:im,1:jm), xlon,xlat)
1572 write(6,*)
' WRITE_NETCDF time= ',tend-tbeg
1573 print *,
' wrote netcdf file out.oro.tile?.nc' 1575 print *,
' ===== Deallocate Arrays and ENDING MTN VAR OROG program' 1578 deallocate(jst,jen,kpds,kgds,numi,lonsperlat)
1579 deallocate(cosclt,wgtclt,rclt,xlat,diffx,xlon,ors,oaa,ola,glat)
1583 deallocate (geolon,geolon_c,geolat,geolat_c)
1584 deallocate (slm,oro,var,orf,land_frac)
1585 deallocate (theta,gamma,sigma,elvmax)
1589 write(6,*)
' Total runtime time= ',tend-tbeg1
1622 SUBROUTINE makemt(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
1623 1 GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
1624 dimension glat(jmn),xlat(jm)
1626 INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
1627 DIMENSION ORO(IM,JM),SLM(IM,JM),VAR(IM,JM),VAR4(IM,JM)
1628 dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
1629 INTEGER mskocn,isave
1637 print *,
' _____ SUBROUTINE MAKEMT ' 1644 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1654 faclon = delx / delxn
1655 ist(i,j) = faclon * float(i-1) - faclon * 0.5 + 1
1656 ien(i,j) = faclon * float(i) - faclon * 0.5 + 1
1660 IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
1661 IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
1670 print *,
' DELX=',delx,
' DELXN=',delxn
1674 xxlat = (xlat(j)+xlat(j+1))/2.
1675 IF(flag.AND.glat(j1).GT.xxlat)
THEN 1683 jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
1684 jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
1702 DO ii1 = 1, ien(i,j) - ist(i,j) + 1
1703 i1 = ist(i,j) + ii1 - 1
1704 IF(i1.LE.0.) i1 = i1 + imn
1705 IF(i1.GT.imn) i1 = i1 - imn
1712 xland = xland + float(zslm(i1,j1))
1713 xwatr = xwatr + float(1-zslm(i1,j1))
1715 height = float(zavg(i1,j1))
1717 IF(height.LT.-990.) height = 0.0
1718 xl1 = xl1 + height * float(zslm(i1,j1))
1719 xs1 = xs1 + height * float(1-zslm(i1,j1))
1721 xw2 = xw2 + height ** 2
1732 IF(xnsum.GT.1.)
THEN 1737 slm(i,j) = float(nint(xland/xnsum))
1738 IF(slm(i,j).NE.0.)
THEN 1739 oro(i,j)= xl1 / xland
1741 oro(i,j)= xs1 / xwatr
1743 var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1744 DO ii1 = 1, ien(i,j) - ist(i,j) + 1
1745 i1 = ist(i,j) + ii1 - 1
1746 IF(i1.LE.0.) i1 = i1 + imn
1747 IF(i1.GT.imn) i1 = i1 - imn
1749 height = float(zavg(i1,j1))
1750 IF(height.LT.-990.) height = 0.0
1751 xw4 = xw4 + (height-oro(i,j)) ** 4
1754 IF(var(i,j).GT.1.)
THEN 1758 var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
1763 WRITE(6,*)
"! MAKEMT ORO SLM VAR VAR4 DONE" 1788 SUBROUTINE get_index(IMN,JMN,npts,lonO,latO,DELXN,
1789 & jst,jen,ilist,numx)
1791 integer,
intent(in) :: IMN,JMN
1793 real,
intent(in) :: LONO(npts), LATO(npts)
1794 real,
intent(in) :: DELXN
1795 integer,
intent(out) :: jst,jen
1796 integer,
intent(out) :: ilist(IMN)
1797 integer,
intent(out) :: numx
1798 real minlat,maxlat,minlon,maxlon
1799 integer i2, ii, ist, ien
1801 minlat = minval(lato)
1802 maxlat = maxval(lato)
1803 minlon = minval(lono)
1804 maxlon = maxval(lono)
1805 ist = minlon/delxn+1
1806 ien = maxlon/delxn+1
1807 jst = (minlat+90)/delxn+1
1808 jen = (maxlat+90)/delxn
1813 if(jen>jmn) jen = jmn
1816 if((jst == 1 .OR. jen == jmn) .and.
1817 & (ien-ist+1 > imn/2) )
then 1822 else if( ien-ist+1 > imn/2 )
then 1828 ii = lono(i2)/delxn+1
1829 if(ii <0 .or. ii>imn) print*,
"ii=",ii,imn,lono(i2),delxn
1830 if( ii < imn/2 )
then 1832 else if( ii > imn/2 )
then 1836 if(ist<1 .OR. ist>imn)
then 1837 print*, .or.
"ist<1 ist>IMN" 1840 if(ien<1 .OR. ien>imn)
then 1841 print*, .or.
"iend<1 iend>IMN" 1845 numx = imn - ien + 1
1847 ilist(i2) = ien + (i2-1)
1856 ilist(i2) = ist + (i2-1)
1882 SUBROUTINE makemt2(ZAVG,ZSLM,ORO,SLM,land_frac,VAR,VAR4,
1883 1 GLAT,IM,JM,IMN,JMN,lon_c,lat_c)
1885 real,
parameter :: D2R = 3.14159265358979/180.
1886 integer,
parameter :: MAXSUM=20000000
1887 real,
dimension(:),
allocatable :: hgt_1d
1888 integer IM, JM, IMN, JMN
1889 real GLAT(JMN), GLON(IMN)
1890 INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
1891 real land_frac(IM,JM)
1892 real ORO(IM,JM),SLM(IM,JM),VAR(IM,JM),VAR4(IM,JM)
1893 integer IST,IEN,JST, JEN
1894 real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
1895 INTEGER mskocn,isave
1897 real LONO(4),LATO(4),LONI,LATI
1899 integer JM1,i,j,nsum,ii,jj,i1,numx,i2
1901 real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1,XW2,XW4
1903 real :: xnsum_j,xland_j,xwatr_j
1904 logical inside_a_polygon
1911 print *,
' _____ SUBROUTINE MAKEMT2 ' 1912 allocate(hgt_1d(maxsum))
1919 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1922 glon(i) = 0. + (i-1) * delxn + delxn * 0.5
1925 land_frac(:,:) = 0.0
1950 lono(1) = lon_c(i,j)
1951 lono(2) = lon_c(i+1,j)
1952 lono(3) = lon_c(i+1,j+1)
1953 lono(4) = lon_c(i,j+1)
1954 lato(1) = lat_c(i,j)
1955 lato(2) = lat_c(i+1,j)
1956 lato(3) = lat_c(i+1,j+1)
1957 lato(4) = lat_c(i,j+1)
1958 call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
1959 do jj = jst, jen;
do i2 = 1, numx
1962 lati = -90 + jj*delxn
1963 if(inside_a_polygon(loni*d2r,lati*d2r,4,
1964 & lono*d2r,lato*d2r))
then 1966 xland = xland + float(zslm(ii,jj))
1967 xwatr = xwatr + float(1-zslm(ii,jj))
1969 height = float(zavg(ii,jj))
1971 if(nsum > maxsum)
then 1972 print*,
"nsum is greater than MAXSUM, increase MAXSUM" 1975 hgt_1d(nsum) = height
1976 IF(height.LT.-990.) height = 0.0
1977 xl1 = xl1 + height * float(zslm(ii,jj))
1978 xs1 = xs1 + height * float(1-zslm(ii,jj))
1980 xw2 = xw2 + height ** 2
1985 IF(xnsum.GT.1.)
THEN 1989 land_frac(i,j) = xland/xnsum
1990 slm(i,j) = float(nint(xland/xnsum))
1991 IF(slm(i,j).NE.0.)
THEN 1992 oro(i,j)= xl1 / xland
1994 oro(i,j)= xs1 / xwatr
1996 var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1998 xw4 = xw4 + (hgt_1d(i1) - oro(i,j)) ** 4
2001 IF(var(i,j).GT.1.)
THEN 2002 var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
2008 WRITE(6,*)
"! MAKEMT2 ORO SLM VAR VAR4 DONE" 2042 SUBROUTINE makepc(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2043 1 GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
2047 parameter(rearth=6.3712e+6)
2048 dimension glat(jmn),xlat(jm),deltax(jmn)
2049 INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2050 DIMENSION ORO(IM,JM),SLM(IM,JM),HL(IM,JM),HK(IM,JM)
2051 DIMENSION HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
2052 dimension theta(im,jm),gamma(im,jm),sigma2(im,jm),sigma(im,jm)
2053 dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
2058 pi = 4.0 * atan(1.0)
2064 deltay = certh / float(jmn)
2065 print *,
'MAKEPC: DELTAY=',deltay
2068 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2069 deltax(j) = deltay * cos(glat(j)*pi/180.0)
2078 faclon = delx / delxn
2079 ist(i,j) = faclon * float(i-1) - faclon * 0.5
2080 ist(i,j) = ist(i,j) + 1
2081 ien(i,j) = faclon * float(i) - faclon * 0.5
2088 IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2089 IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2091 if ( i .lt. 10 .and. j .lt. 10 )
2092 1 print*,
' I j IST IEN ',i,j,ist(i,j),ien(i,j)
2094 IF (ien(i,j) .LT. ist(i,j))
2095 1 print *,
' MAKEPC: IEN < IST: I,J,IST(I,J),IEN(I,J)',
2096 2 i,j,ist(i,j),ien(i,j)
2102 xxlat = (xlat(j)+xlat(j+1))/2.
2103 IF(flag.AND.glat(j1).GT.xxlat)
THEN 2110 jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2111 jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2113 print*,
' IST,IEN(1,1-numi(1,JM))',ist(1,1),ien(1,1),
2114 1 ist(numi(jm),jm),ien(numi(jm),jm), numi(jm)
2115 print*,
' JST,JEN(1,JM) ',jst(1),jen(1),jst(jm),jen(jm)
2144 DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2145 i1 = ist(i,j) + ii1 - 1
2146 IF(i1.LE.0.) i1 = i1 + imn
2147 IF(i1.GT.imn) i1 = i1 - imn
2152 if (i1 - 1 .le. 0 ) i0 = i0 + imn
2153 if (i1 - 1 .gt. imn) i0 = i0 - imn
2156 if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2157 if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2161 if ( i1 .eq. ist(i,j) .and. j1 .eq. jst(j) )
2162 1 print*,
' J, J1,IST,JST,DELTAX,GLAT ',
2163 2 j,j1,ist(i,j),jst(j),deltax(j1),glat(j1)
2164 if ( i1 .eq. ien(i,j) .and. j1 .eq. jen(j) )
2165 1 print*,
' J, J1,IEN,JEN,DELTAX,GLAT ',
2166 2 j,j1,ien(i,j),jen(j),deltax(j1),glat(j1)
2168 xland = xland + float(zslm(i1,j1))
2169 xwatr = xwatr + float(1-zslm(i1,j1))
2172 height = float(zavg(i1,j1))
2173 hi0 = float(zavg(i0,j1))
2174 hip1 = float(zavg(ip1,j1))
2176 IF(height.LT.-990.) height = 0.0
2177 if(hi0 .lt. -990.) hi0 = 0.0
2178 if(hip1 .lt. -990.) hip1 = 0.0
2180 xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2181 xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2185 if ( j1 .ne. jst(jm) .and. j1 .ne. jen(1) )
then 2186 hj0 = float(zavg(i1,j1-1))
2187 hjp1 = float(zavg(i1,j1+1))
2188 if(hj0 .lt. -990.) hj0 = 0.0
2189 if(hjp1 .lt. -990.) hjp1 = 0.0
2191 yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2192 yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2198 elseif ( j1 .eq. jst(jm) )
then 2200 if (ijax .le. 0 ) ijax = ijax + imn
2201 if (ijax .gt. imn) ijax = ijax - imn
2203 hijax = float(zavg(ijax,j1))
2204 hi1j1 = float(zavg(i1,j1))
2205 if(hijax .lt. -990.) hijax = 0.0
2206 if(hi1j1 .lt. -990.) hi1j1 = 0.0
2208 yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2209 yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2215 elseif ( j1 .eq. jen(1) )
then 2217 if (ijax .le. 0 ) ijax = ijax + imn
2218 if (ijax .gt. imn) ijax = ijax - imn
2219 hijax = float(zavg(ijax,j1))
2220 hi1j1 = float(zavg(i1,j1))
2221 if(hijax .lt. -990.) hijax = 0.0
2222 if(hi1j1 .lt. -990.) hi1j1 = 0.0
2223 if ( i1 .lt. 5 )print *,
' S.Pole i1,j1 :',i1,j1,hijax,hi1j1
2225 yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2226 yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2233 xfpyfp = xfpyfp + xfp * yfp
2234 xl1 = xl1 + height * float(zslm(i1,j1))
2235 xs1 = xs1 + height * float(1-zslm(i1,j1))
2245 IF(xnsum.GT.1.)
THEN 2246 slm(i,j) = float(nint(xland/xnsum))
2247 IF(slm(i,j).NE.0.)
THEN 2248 oro(i,j)= xl1 / xland
2249 hx2(i,j) = xfp2 / xland
2250 hy2(i,j) = yfp2 / xland
2251 hxy(i,j) = xfpyfp / xland
2253 oro(i,j)= xs1 / xwatr
2257 print *,
" I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2259 print *,
" xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2260 print *,
" HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2266 hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2267 hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2268 hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2269 IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. )
THEN 2271 theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) * 180.0 / pi
2275 sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2276 if ( sigma2(i,j) .GE. 0. )
then 2277 sigma(i,j) = sqrt(sigma2(i,j) )
2278 if (sigma2(i,j) .ne. 0. .and.
2279 & hk(i,j) .GE. hlprim(i,j) )
2280 1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2286 print *,
" I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2287 1 sigma(i,j),gamma(i,j)
2288 print *,
" HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2292 WRITE(6,*)
"! MAKE Principal Coord DONE" 2316 SUBROUTINE makepc2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2317 1 GLAT,IM,JM,IMN,JMN,lon_c,lat_c)
2322 real,
parameter :: REARTH=6.3712e+6
2323 real,
parameter :: D2R = 3.14159265358979/180.
2324 integer :: im,jm,imn,jmn
2325 real :: GLAT(JMN),DELTAX(JMN)
2326 INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2327 real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
2328 real ORO(IM,JM),SLM(IM,JM),HL(IM,JM),HK(IM,JM)
2329 real HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
2330 real THETA(IM,JM),GAMMA(IM,JM),SIGMA2(IM,JM),SIGMA(IM,JM)
2331 real PI,CERTH,DELXN,DELTAY,XNSUM,XLAND,XWATR,XL1,XS1
2332 real xfp,yfp,xfpyfp,xfp2,yfp2,HEIGHT
2333 real hi0,hip1,hj0,hjp1,hijax,hi1j1
2334 real LONO(4),LATO(4),LONI,LATI
2335 integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
2337 logical inside_a_polygon
2342 pi = 4.0 * atan(1.0)
2347 deltay = certh / float(jmn)
2348 print *,
'MAKEPC2: DELTAY=',deltay
2351 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2352 deltax(j) = deltay * cos(glat(j)*d2r)
2390 lono(1) = lon_c(i,j)
2391 lono(2) = lon_c(i+1,j)
2392 lono(3) = lon_c(i+1,j+1)
2393 lono(4) = lon_c(i,j+1)
2394 lato(1) = lat_c(i,j)
2395 lato(2) = lat_c(i+1,j)
2396 lato(3) = lat_c(i+1,j+1)
2397 lato(4) = lat_c(i,j+1)
2398 call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
2400 do j1 = jst, jen;
do i2 = 1, numx
2403 lati = -90 + j1*delxn
2404 if(inside_a_polygon(loni*d2r,lati*d2r,4,
2405 & lono*d2r,lato*d2r))
then 2410 if (i1 - 1 .le. 0 ) i0 = i0 + imn
2411 if (i1 - 1 .gt. imn) i0 = i0 - imn
2414 if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2415 if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2417 xland = xland + float(zslm(i1,j1))
2418 xwatr = xwatr + float(1-zslm(i1,j1))
2421 height = float(zavg(i1,j1))
2422 hi0 = float(zavg(i0,j1))
2423 hip1 = float(zavg(ip1,j1))
2425 IF(height.LT.-990.) height = 0.0
2426 if(hi0 .lt. -990.) hi0 = 0.0
2427 if(hip1 .lt. -990.) hip1 = 0.0
2429 xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2430 xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2434 if ( j1 .ne. 1 .and. j1 .ne. jmn )
then 2435 hj0 = float(zavg(i1,j1-1))
2436 hjp1 = float(zavg(i1,j1+1))
2437 if(hj0 .lt. -990.) hj0 = 0.0
2438 if(hjp1 .lt. -990.) hjp1 = 0.0
2440 yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2441 yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2447 elseif ( j1 .eq. 1 )
then 2449 if (ijax .le. 0 ) ijax = ijax + imn
2450 if (ijax .gt. imn) ijax = ijax - imn
2452 hijax = float(zavg(ijax,j1))
2453 hi1j1 = float(zavg(i1,j1))
2454 if(hijax .lt. -990.) hijax = 0.0
2455 if(hi1j1 .lt. -990.) hi1j1 = 0.0
2457 yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2458 yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2464 elseif ( j1 .eq. jmn )
then 2466 if (ijax .le. 0 ) ijax = ijax + imn
2467 if (ijax .gt. imn) ijax = ijax - imn
2468 hijax = float(zavg(ijax,j1))
2469 hi1j1 = float(zavg(i1,j1))
2470 if(hijax .lt. -990.) hijax = 0.0
2471 if(hi1j1 .lt. -990.) hi1j1 = 0.0
2472 if ( i1 .lt. 5 )print *,
' S.Pole i1,j1 :',i1,j1,
2475 yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2476 yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2483 xfpyfp = xfpyfp + xfp * yfp
2484 xl1 = xl1 + height * float(zslm(i1,j1))
2485 xs1 = xs1 + height * float(1-zslm(i1,j1))
2496 IF(xnsum.GT.1.)
THEN 2497 slm(i,j) = float(nint(xland/xnsum))
2498 IF(slm(i,j).NE.0.)
THEN 2499 oro(i,j)= xl1 / xland
2500 hx2(i,j) = xfp2 / xland
2501 hy2(i,j) = yfp2 / xland
2502 hxy(i,j) = xfpyfp / xland
2504 oro(i,j)= xs1 / xwatr
2508 print *,
" I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2510 print *,
" xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2511 print *,
" HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2517 hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2518 hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2519 hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2520 IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. )
THEN 2522 theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) / d2r
2526 sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2527 if ( sigma2(i,j) .GE. 0. )
then 2528 sigma(i,j) = sqrt(sigma2(i,j) )
2529 if (sigma2(i,j) .ne. 0. .and.
2530 & hk(i,j) .GE. hlprim(i,j) )
2531 1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2537 print *,
" I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2538 1 sigma(i,j),gamma(i,j)
2539 print *,
" HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2544 WRITE(6,*)
"! MAKE Principal Coord DONE" 2547 END SUBROUTINE makepc2
2590 SUBROUTINE makeoa(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2591 1 ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
2592 2 IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
2593 dimension glat(jmn),xlat(jm)
2594 INTEGER ZAVG(IMN,JMN)
2595 DIMENSION ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
2596 DIMENSION OA4(IM,JM,4),IOA4(IM,JM,4)
2597 DIMENSION IST(IM,jm),IEN(IM,jm),JST(JM),JEN(JM)
2598 dimension xnsum(im,jm),xnsum1(im,jm),xnsum2(im,jm)
2599 dimension xnsum3(im,jm),xnsum4(im,jm)
2600 dimension var(im,jm),ol(im,jm,4),numi(jm)
2610 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2612 print *,
' IM=',im,
' JM=',jm,
' IMN=',imn,
' JMN=',jmn
2619 faclon = delx / delxn
2621 ist(i,j) = faclon * float(i-1) - faclon * 0.5
2622 ist(i,j) = ist(i,j) + 1
2623 ien(i,j) = faclon * float(i) - faclon * 0.5
2626 IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2627 IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2629 if ( i .lt. 3 .and. j .lt. 3 )
2630 1print*,
' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2631 if ( i .lt. 3 .and. j .ge. jm-1 )
2632 1print*,
' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2635 print *,
'MAKEOA: DELXN,DELX,FACLON',delxn,delx,faclon
2636 print *,
' ***** ready to start JST JEN section ' 2641 xxlat = (xlat(j)+xlat(j+1))/2.
2642 IF(flag.AND.glat(j1).GT.xxlat)
THEN 2647 1print*,
' MAKEOA: XX j JST JEN ',j,jst(j),jen(j)
2651 1print*,
' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2653 1print*,
' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2662 jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2663 jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2664 print *,
' ***** JST(1) JEN(1) ',jst(1),jen(1)
2665 print *,
' ***** JST(JM) JEN(JM) ',jst(jm),jen(jm)
2670 elvmax(i,j) = oro(i,j)
2679 DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2680 i1 = ist(i,j) + ii1 - 1
2682 IF(i1.LE.0.) i1 = i1 + imn
2683 IF (i1 .GT. imn) i1 = i1 - imn
2685 height = float(zavg(i1,j1))
2686 IF(height.LT.-990.) height = 0.0
2687 IF ( height .gt. oro(i,j) )
then 2688 if ( height .gt. zmax(i,j) )zmax(i,j) = height
2689 xnsum(i,j) = xnsum(i,j) + 1
2693 if ( i .lt. 5 .and. j .ge. jm-5 )
then 2694 print *,
' I,J,ORO(I,J),XNSUM(I,J),ZMAX(I,J):',
2695 1 i,j,oro(i,j),xnsum(i,j),zmax(i,j)
2706 oro1(i,j) = oro(i,j)
2707 elvmax(i,j) = zmax(i,j)
2740 hc = 1116.2 - 0.878 * var(i,j)
2742 DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2743 i1 = ist(i,j) + ii1 - 1
2748 IF(i1.GT.imn) i1 = i1 - imn
2750 IF(float(zavg(i1,j1)) .GT. hc)
2751 1 xnsum1(i,j) = xnsum1(i,j) + 1
2752 xnsum2(i,j) = xnsum2(i,j) + 1
2756 inci = nint((ien(i,j)-ist(i,j)) * 0.5)
2757 isttt = min(max(ist(i,j)-inci,1),imn)
2758 ieddd = min(max(ien(i,j)-inci,1),imn)
2760 incj = nint((jen(j)-jst(j)) * 0.5)
2761 jsttt = min(max(jst(j)-incj,1),jmn)
2762 jeddd = min(max(jen(j)-incj,1),jmn)
2772 IF(float(zavg(i1,j1)) .GT. hc)
2773 1 xnsum3(i,j) = xnsum3(i,j) + 1
2774 xnsum4(i,j) = xnsum4(i,j) + 1
2800 IF (ii .GT. numi(j)) ii = ii - numi(j)
2801 xnpu = xnsum(i,j) + xnsum(i,j+1)
2802 xnpd = xnsum(ii,j) + xnsum(ii,j+1)
2803 IF (xnpd .NE. xnpu) oa4(ii,j+1,1) = 1. - xnpd / max(xnpu , 1.)
2804 ol(ii,j+1,1) = (xnsum3(i,j+1)+xnsum3(ii,j+1))/
2805 1 (xnsum4(i,j+1)+xnsum4(ii,j+1))
2820 IF (ii .GT. numi(j)) ii = ii - numi(j)
2821 xnpu = xnsum(i,j+1) + xnsum(ii,j+1)
2822 xnpd = xnsum(i,j) + xnsum(ii,j)
2823 IF (xnpd .NE. xnpu) oa4(ii,j+1,2) = 1. - xnpd / max(xnpu , 1.)
2824 ol(ii,j+1,2) = (xnsum3(ii,j)+xnsum3(ii,j+1))/
2825 1 (xnsum4(ii,j)+xnsum4(ii,j+1))
2831 IF (ii .GT. numi(j)) ii = ii - numi(j)
2832 xnpu = xnsum(i,j+1) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2833 xnpd = xnsum(ii,j) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2834 IF (xnpd .NE. xnpu) oa4(ii,j+1,3) = 1. - xnpd / max(xnpu , 1.)
2835 ol(ii,j+1,3) = (xnsum1(ii,j)+xnsum1(i,j+1))/
2836 1 (xnsum2(ii,j)+xnsum2(i,j+1))
2842 IF (ii .GT. numi(j)) ii = ii - numi(j)
2843 xnpu = xnsum(i,j) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2844 xnpd = xnsum(ii,j+1) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2845 IF (xnpd .NE. xnpu) oa4(ii,j+1,4) = 1. - xnpd / max(xnpu , 1.)
2846 ol(ii,j+1,4) = (xnsum1(i,j)+xnsum1(ii,j+1))/
2847 1 (xnsum2(i,j)+xnsum2(ii,j+1))
2853 ol(i,1,kwd) = ol(i,2,kwd)
2854 ol(i,jm,kwd) = ol(i,jm-1,kwd)
2862 oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
2877 t = abs( oa4(i,j,kwd) )
2881 ELSE IF(t .GT. 0. .AND. t .LE. 1.)
THEN 2884 ELSE IF(t .GT. 1. .AND. t .LE. 10.)
THEN 2887 ELSE IF(t .GT. 10. .AND. t .LE. 100.)
THEN 2890 ELSE IF(t .GT. 100. .AND. t .LE. 1000.)
THEN 2893 ELSE IF(t .GT. 1000. .AND. t .LE. 10000.)
THEN 2896 ELSE IF(t .GT. 10000.)
THEN 2904 WRITE(6,*)
"! MAKEOA EXIT" 2907 END SUBROUTINE makeoa
2920 real dx, lat, degrad
2923 real,
parameter :: radius = 6371200
2925 get_lon_angle = 2*asin( sin(dx/radius*0.5)/cos(lat) )*degrad
2942 real,
parameter :: radius = 6371200
2984 SUBROUTINE makeoa2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2985 1 ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
2986 2 IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,dx,dy,
2987 3 is_south_pole,is_north_pole )
2989 real,
parameter :: MISSING_VALUE = -9999.
2990 real,
parameter :: d2r = 3.14159265358979/180.
2991 real,
PARAMETER :: R2D=180./3.14159265358979
2992 integer im,jm,imn,jmn
2994 INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2995 real ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
2997 integer IOA4(IM,JM,4)
2998 real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
2999 real lon_t(IM,JM), lat_t(IM,JM)
3000 real dx(IM,JM), dy(IM,JM)
3001 logical is_south_pole(IM,JM), is_north_pole(IM,JM)
3002 real XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
3003 real XNSUM3(IM,JM),XNSUM4(IM,JM)
3004 real VAR(IM,JM),OL(IM,JM,4)
3006 integer i,j,ilist(IMN),numx,i1,j1,ii1
3008 real LONO(4),LATO(4),LONI,LATI
3009 real DELXN,HC,HEIGHT,XNPU,XNPD,T
3010 integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
3011 logical inside_a_polygon
3012 real lon,lat,dlon,dlat,dlat_old
3013 real lon1,lat1,lon2,lat2
3014 real xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
3015 real HC_11, HC_12, HC_21, HC_22
3016 real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
3017 real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
3018 real get_lon_angle, get_lat_angle, get_xnsum
3019 integer ist, ien, jst, jen
3020 real xland,xwatr,xl1,xs1,oroavg,slm
3027 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3029 print *,
' IM=',im,
' JM=',jm,
' IMN=',imn,
' JMN=',jmn
3037 elvmax(i,j) = oro(i,j)
3045 oro1(i,j) = oro(i,j)
3046 elvmax(i,j) = zmax(i,j)
3058 hc = 1116.2 - 0.878 * var(i,j)
3059 lono(1) = lon_c(i,j)
3060 lono(2) = lon_c(i+1,j)
3061 lono(3) = lon_c(i+1,j+1)
3062 lono(4) = lon_c(i,j+1)
3063 lato(1) = lat_c(i,j)
3064 lato(2) = lat_c(i+1,j)
3065 lato(3) = lat_c(i+1,j+1)
3066 lato(4) = lat_c(i,j+1)
3067 call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3068 do j1 = jst, jen;
do ii1 = 1, numx
3071 lati = -90 + j1*delxn
3072 if(inside_a_polygon(loni*d2r,lati*d2r,4,
3073 & lono*d2r,lato*d2r))
then 3075 height = float(zavg(i1,j1))
3076 IF(height.LT.-990.) height = 0.0
3077 IF ( height .gt. oro(i,j) )
then 3078 if ( height .gt. zmax(i,j) )zmax(i,j) = height
3091 oro1(i,j) = oro(i,j)
3092 elvmax(i,j) = zmax(i,j)
3126 if(is_north_pole(i,j))
then 3127 print*,
"set oa1 = 0 and ol=0 at i,j=", i,j
3132 else if(is_south_pole(i,j))
then 3133 print*,
"set oa1 = 0 and ol=1 at i,j=", i,j
3141 dlon = get_lon_angle(dx(i,j), lat*d2r, r2d )
3142 dlat = get_lat_angle(dy(i,j), r2d)
3144 if( lat-dlat*0.5<-90.)
then 3145 print*,
"at i,j =", i,j, lat, dlat, lat-dlat*0.5
3146 print*,
"ERROR: lat-dlat*0.5<-90." 3149 if( lat+dlat*2 > 90.)
then 3152 print*,
"at i,j=",i,j,
" adjust dlat from ",
3153 & dlat_old,
" to ", dlat
3161 if(lat1<-90 .or. lat2>90)
then 3162 print*,
"at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3164 xnsum11 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3172 if(lat1<-90 .or. lat2>90)
then 3173 print*,
"at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3175 xnsum12 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3183 if(lat1<-90 .or. lat2>90)
then 3184 print*,
"at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3186 xnsum21 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3194 if(lat1<-90 .or. lat2>90)
then 3195 print*,
"at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3198 xnsum22 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3201 xnpu = xnsum11 + xnsum12
3202 xnpd = xnsum21 + xnsum22
3203 IF (xnpd .NE. xnpu) oa4(i,j,1) = 1. - xnpd / max(xnpu , 1.)
3205 xnpu = xnsum11 + xnsum21
3206 xnpd = xnsum12 + xnsum22
3207 IF (xnpd .NE. xnpu) oa4(i,j,2) = 1. - xnpd / max(xnpu , 1.)
3209 xnpu = xnsum11 + (xnsum12+xnsum21)*0.5
3210 xnpd = xnsum22 + (xnsum12+xnsum21)*0.5
3211 IF (xnpd .NE. xnpu) oa4(i,j,3) = 1. - xnpd / max(xnpu , 1.)
3213 xnpu = xnsum12 + (xnsum11+xnsum22)*0.5
3214 xnpd = xnsum21 + (xnsum11+xnsum22)*0.5
3215 IF (xnpd .NE. xnpu) oa4(i,j,4) = 1. - xnpd / max(xnpu , 1.)
3224 if(lat1<-90 .or. lat2>90)
then 3225 print*,
"at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3227 call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3228 & zavg,zslm,delxn, xnsum1_11, xnsum2_11, hc_11)
3235 if(lat1<-90 .or. lat2>90)
then 3236 print*,
"at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3238 call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3239 & zavg,zslm,delxn, xnsum1_12, xnsum2_12, hc_12)
3246 if(lat1<-90 .or. lat2>90)
then 3247 print*,
"at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3249 call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3250 & zavg,zslm,delxn, xnsum1_21, xnsum2_21, hc_21)
3257 if(lat1<-90 .or. lat2>90)
then 3258 print*,
"at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3260 call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3261 & zavg,zslm,delxn, xnsum1_22, xnsum2_22, hc_22)
3263 ol(i,j,3) = (xnsum1_22+xnsum1_11)/(xnsum2_22+xnsum2_11)
3264 ol(i,j,4) = (xnsum1_12+xnsum1_21)/(xnsum2_12+xnsum2_21)
3272 if(lat1<-90 .or. lat2>90)
then 3273 print*,
"at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3275 call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3276 & zavg,zslm,delxn, xnsum1_11, xnsum2_11, hc_11)
3283 if(lat1<-90 .or. lat2>90)
then 3284 print*,
"at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3287 call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3288 & zavg,zslm,delxn, xnsum1_12, xnsum2_12, hc_12)
3295 if(lat1<-90 .or. lat2>90)
then 3296 print*,
"at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3298 call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3299 & zavg,zslm,delxn, xnsum1_21, xnsum2_21, hc_21)
3306 if(lat1<-90 .or. lat2>90)
then 3307 print*,
"at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3310 call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3311 & zavg,zslm,delxn, xnsum1_22, xnsum2_22, hc_22)
3313 ol(i,j,1) = (xnsum1_11+xnsum1_21)/(xnsum2_11+xnsum2_21)
3314 ol(i,j,2) = (xnsum1_21+xnsum1_22)/(xnsum2_21+xnsum2_22)
3323 oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3338 t = abs( oa4(i,j,kwd) )
3342 ELSE IF(t .GT. 0. .AND. t .LE. 1.)
THEN 3345 ELSE IF(t .GT. 1. .AND. t .LE. 10.)
THEN 3348 ELSE IF(t .GT. 10. .AND. t .LE. 100.)
THEN 3351 ELSE IF(t .GT. 100. .AND. t .LE. 1000.)
THEN 3354 ELSE IF(t .GT. 1000. .AND. t .LE. 10000.)
THEN 3357 ELSE IF(t .GT. 10000.)
THEN 3365 WRITE(6,*)
"! MAKEOA2 EXIT" 3369 END SUBROUTINE makeoa2
3381 real,
intent(in) :: theta1, phi1, theta2, phi2
3384 if(theta1 == theta2 .and. phi1 == phi2)
then 3389 dot = cos(phi1)*cos(phi2)*cos(theta1-theta2) + sin(phi1)*sin(phi2)
3390 if(dot > 1. ) dot = 1.
3391 if(dot < -1.) dot = -1.
3415 & bitmap_in,num_out, lon_out,lat_out, iindx, jindx )
3416 integer,
intent(in) :: im_in, jm_in, num_out
3417 real,
intent(in) :: geolon_in(im_in,jm_in)
3418 real,
intent(in) :: geolat_in(im_in,jm_in)
3419 logical*1,
intent(in) :: bitmap_in(im_in,jm_in)
3420 real,
intent(in) :: lon_out(num_out), lat_out(num_out)
3421 integer,
intent(out):: iindx(num_out), jindx(num_out)
3422 real,
parameter :: MAX_DIST = 1.e+20
3423 integer,
parameter :: NUMNBR = 20
3424 integer :: i_c,j_c,jstart,jend
3427 print*,
"im_in,jm_in = ", im_in, jm_in
3428 print*,
"num_out = ", num_out
3429 print*,
"max and min of lon_in is", minval(geolon_in),
3431 print*,
"max and min of lat_in is", minval(geolat_in),
3433 print*,
"max and min of lon_out is", minval(lon_out),
3435 print*,
"max and min of lat_out is", minval(lat_out),
3437 print*,
"count(bitmap_in)= ", count(bitmap_in), max_dist
3446 if(lat .LE. geolat_in(1,j) .and.
3447 & lat .GE. geolat_in(1,j+1))
then 3451 if(lat > geolat_in(1,1)) j_c = 1
3452 if(lat < geolat_in(1,jm_in)) j_c = jm_in
3455 jstart = max(j_c-numnbr,1)
3456 jend = min(j_c+numnbr,jm_in)
3461 do j = jstart, jend;
do i = 1,im_in
3462 if(bitmap_in(i,j) )
then 3465 & geolon_in(i,j), geolat_in(i,j))
3467 iindx(n) = i; jindx(n) = j
3472 if(iindx(n) ==0)
then 3473 print*,
"lon,lat=", lon,lat
3474 print*,
"jstart, jend=", jstart, jend, dist
3475 print*,
"ERROR in get mismatch_index: not find nearest points" 3496 & num_out, data_out, iindx, jindx)
3497 integer,
intent(in) :: im_in, jm_in, num_out
3498 real,
intent(in) :: data_in(im_in,jm_in)
3499 real,
intent(out):: data_out(num_out)
3500 integer,
intent(in) :: iindx(num_out), jindx(num_out)
3503 data_out(n) = data_in(iindx(n),jindx(n))
3552 SUBROUTINE makeoa3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
3553 1 ORO,SLM,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
3554 2 IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,
3555 3 is_south_pole,is_north_pole,IMI,JMI,OA_IN,OL_IN,
3556 4 slm_in,lon_in,lat_in)
3564 real,
parameter :: MISSING_VALUE = -9999.
3565 real,
parameter :: d2r = 3.14159265358979/180.
3566 real,
PARAMETER :: R2D=180./3.14159265358979
3567 integer im,jm,imn,jmn,imi,jmi
3569 INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
3571 real ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
3573 integer IOA4(IM,JM,4)
3574 real OA_IN(IMI,JMI,4), OL_IN(IMI,JMI,4)
3575 real slm_in(IMI,JMI)
3576 real lon_in(IMI,JMI), lat_in(IMI,JMI)
3577 real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
3578 real lon_t(IM,JM), lat_t(IM,JM)
3579 logical is_south_pole(IM,JM), is_north_pole(IM,JM)
3580 real XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
3581 real XNSUM3(IM,JM),XNSUM4(IM,JM)
3582 real VAR(IM,JM),OL(IM,JM,4)
3584 integer i,j,ilist(IMN),numx,i1,j1,ii1
3586 real LONO(4),LATO(4),LONI,LATI
3587 real DELXN,HC,HEIGHT,XNPU,XNPD,T
3588 integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
3589 logical inside_a_polygon
3590 real lon,lat,dlon,dlat,dlat_old
3591 real lon1,lat1,lon2,lat2
3592 real xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
3593 real HC_11, HC_12, HC_21, HC_22
3594 real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
3595 real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
3596 real get_lon_angle, get_lat_angle, get_xnsum
3597 integer ist, ien, jst, jen
3598 real xland,xwatr,xl1,xs1,oroavg
3599 integer int_opt, ipopt(20), kgds_input(200), kgds_output(200)
3600 integer count_land_output
3601 integer ij, ijmdl_output, iret, num_mismatch_land, num
3602 integer ibo(1), ibi(1)
3603 logical*1,
allocatable :: bitmap_input(:,:)
3604 logical*1,
allocatable :: bitmap_output(:,:)
3605 integer,
allocatable :: ijsav_land_output(:)
3606 real,
allocatable :: lats_land_output(:)
3607 real,
allocatable :: lons_land_output(:)
3608 real,
allocatable :: output_data_land(:,:)
3609 real,
allocatable :: lons_mismatch_output(:)
3610 real,
allocatable :: lats_mismatch_output(:)
3611 real,
allocatable :: data_mismatch_output(:)
3612 integer,
allocatable :: iindx(:), jindx(:)
3618 ijmdl_output = im*jm
3621 glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3623 print *,
' IM=',im,
' JM=',jm,
' IMN=',imn,
' JMN=',jmn
3631 elvmax(i,j) = oro(i,j)
3639 oro1(i,j) = oro(i,j)
3640 elvmax(i,j) = zmax(i,j)
3649 hc = 1116.2 - 0.878 * var(i,j)
3650 lono(1) = lon_c(i,j)
3651 lono(2) = lon_c(i+1,j)
3652 lono(3) = lon_c(i+1,j+1)
3653 lono(4) = lon_c(i,j+1)
3654 lato(1) = lat_c(i,j)
3655 lato(2) = lat_c(i+1,j)
3656 lato(3) = lat_c(i+1,j+1)
3657 lato(4) = lat_c(i,j+1)
3658 call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3659 do j1 = jst, jen;
do ii1 = 1, numx
3662 lati = -90 + j1*delxn
3663 if(inside_a_polygon(loni*d2r,lati*d2r,4,
3664 & lono*d2r,lato*d2r))
then 3666 height = float(zavg(i1,j1))
3667 IF(height.LT.-990.) height = 0.0
3668 IF ( height .gt. oro(i,j) )
then 3669 if ( height .gt. zmax(i,j) )zmax(i,j) = height
3682 oro1(i,j) = oro(i,j)
3683 elvmax(i,j) = zmax(i,j)
3703 kgds_input(4) = 90000
3706 kgds_input(7) = -90000
3707 kgds_input(8) = nint(-360000./imi)
3708 kgds_input(9) = nint((360.0 / float(imi))*1000.0)
3710 kgds_input(10) = jmi /2
3711 kgds_input(12) = 255
3712 kgds_input(20) = 255
3719 kgds_output(4) = 90000
3721 kgds_output(6) = 128
3722 kgds_output(7) = -90000
3723 kgds_output(8) = nint(-360000./im)
3724 kgds_output(9) = nint((360.0 / float(im))*1000.0)
3726 kgds_output(10) = jm /2
3727 kgds_output(12) = 255
3728 kgds_output(20) = 255
3731 print*,
"sum(slm) = ", sum(slm)
3732 do ij = 1, ijmdl_output
3734 i = mod(ij-1,im) + 1
3735 if (slm(i,j) > 0.0)
then 3736 count_land_output=count_land_output+1
3739 allocate(bitmap_input(imi,jmi))
3740 bitmap_input=.false.
3741 print*,
"number of land input=", sum(slm_in)
3742 where(slm_in > 0.0) bitmap_input=.true.
3743 print*,
"count(bitmap_input)", count(bitmap_input)
3745 allocate(bitmap_output(count_land_output,1))
3746 allocate(output_data_land(count_land_output,1))
3747 allocate(ijsav_land_output(count_land_output))
3748 allocate(lats_land_output(count_land_output))
3749 allocate(lons_land_output(count_land_output))
3754 do ij = 1, ijmdl_output
3756 i = mod(ij-1,im) + 1
3757 if (slm(i,j) > 0.0)
then 3758 count_land_output=count_land_output+1
3759 ijsav_land_output(count_land_output)=ij
3760 lats_land_output(count_land_output)=lat_t(i,j)
3761 lons_land_output(count_land_output)=lon_t(i,j)
3770 bitmap_output = .false.
3771 output_data_land = 0.0
3772 call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3773 & (imi*jmi), count_land_output,
3774 & 1, ibi, bitmap_input, oa_in(:,:,kwd),
3775 & count_land_output, lats_land_output,
3776 & lons_land_output, ibo,
3777 & bitmap_output, output_data_land, iret)
3779 print*,
'- ERROR IN IPOLATES ',iret
3783 num_mismatch_land = 0
3784 do ij = 1, count_land_output
3785 if (bitmap_output(ij,1))
then 3786 j = (ijsav_land_output(ij)-1)/im + 1
3787 i = mod(ijsav_land_output(ij)-1,im) + 1
3788 oa4(i,j,kwd)=output_data_land(ij,1)
3790 num_mismatch_land = num_mismatch_land + 1
3793 print*,
"num_mismatch_land = ", num_mismatch_land
3795 if(.not.
allocated(data_mismatch_output))
then 3796 allocate(lons_mismatch_output(num_mismatch_land))
3797 allocate(lats_mismatch_output(num_mismatch_land))
3798 allocate(data_mismatch_output(num_mismatch_land))
3799 allocate(iindx(num_mismatch_land))
3800 allocate(jindx(num_mismatch_land))
3803 do ij = 1, count_land_output
3804 if (.not. bitmap_output(ij,1))
then 3806 lons_mismatch_output(num) = lons_land_output(ij)
3807 lats_mismatch_output(num) = lats_land_output(ij)
3813 print*,
"before get_mismatch_index", count(bitmap_input)
3815 & lat_in*d2r,bitmap_input,num_mismatch_land,
3816 & lons_mismatch_output*d2r,lats_mismatch_output*d2r,
3820 data_mismatch_output = 0
3822 & num_mismatch_land,data_mismatch_output,iindx,jindx)
3825 do ij = 1, count_land_output
3826 if (.not. bitmap_output(ij,1))
then 3828 j = (ijsav_land_output(ij)-1)/im + 1
3829 i = mod(ijsav_land_output(ij)-1,im) + 1
3830 oa4(i,j,kwd) = data_mismatch_output(num)
3831 if(i==372 .and. j== 611)
then 3832 print*,
"ij=",ij, num, oa4(i,j,kwd),iindx(num),jindx(num)
3842 bitmap_output = .false.
3843 output_data_land = 0.0
3844 call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3845 & (imi*jmi), count_land_output,
3846 & 1, ibi, bitmap_input, ol_in(:,:,kwd),
3847 & count_land_output, lats_land_output,
3848 & lons_land_output, ibo,
3849 & bitmap_output, output_data_land, iret)
3851 print*,
'- ERROR IN IPOLATES ',iret
3855 num_mismatch_land = 0
3856 do ij = 1, count_land_output
3857 if (bitmap_output(ij,1))
then 3858 j = (ijsav_land_output(ij)-1)/im + 1
3859 i = mod(ijsav_land_output(ij)-1,im) + 1
3860 ol(i,j,kwd)=output_data_land(ij,1)
3862 num_mismatch_land = num_mismatch_land + 1
3865 print*,
"num_mismatch_land = ", num_mismatch_land
3867 data_mismatch_output = 0
3869 & num_mismatch_land,data_mismatch_output,iindx,jindx)
3872 do ij = 1, count_land_output
3873 if (.not. bitmap_output(ij,1))
then 3875 j = (ijsav_land_output(ij)-1)/im + 1
3876 i = mod(ijsav_land_output(ij)-1,im) + 1
3877 ol(i,j,kwd) = data_mismatch_output(num)
3878 if(i==372 .and. j== 611)
then 3879 print*,
"ij=",ij, num, ol(i,j,kwd),iindx(num),jindx(num)
3887 deallocate(lons_mismatch_output,lats_mismatch_output)
3888 deallocate(data_mismatch_output,iindx,jindx)
3889 deallocate(bitmap_output,output_data_land)
3890 deallocate(ijsav_land_output,lats_land_output)
3891 deallocate(lons_land_output)
3897 oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3912 t = abs( oa4(i,j,kwd) )
3916 ELSE IF(t .GT. 0. .AND. t .LE. 1.)
THEN 3919 ELSE IF(t .GT. 1. .AND. t .LE. 10.)
THEN 3922 ELSE IF(t .GT. 10. .AND. t .LE. 100.)
THEN 3925 ELSE IF(t .GT. 100. .AND. t .LE. 1000.)
THEN 3928 ELSE IF(t .GT. 1000. .AND. t .LE. 10000.)
THEN 3931 ELSE IF(t .GT. 10000.)
THEN 3939 WRITE(6,*)
"! MAKEOA3 EXIT" 3942 END SUBROUTINE makeoa3
3954 SUBROUTINE revers(IM, JM, numi, F, WRK)
3956 REAL F(IM,JM), WRK(IM,JM)
3966 if (ir .gt. im) ir = ir - im
3991 subroutine rg2gg(im,jm,numi,a)
3993 integer,
intent(in):: im,jm,numi(jm)
3994 real,
intent(inout):: a(im,jm)
3998 r=
real(numi(j))/
real(im)
4000 ir=mod(nint((ig-1)*r),numi(j))+1
4017 subroutine gg2rg(im,jm,numi,a)
4019 integer,
intent(in):: im,jm,numi(jm)
4020 real,
intent(inout):: a(im,jm)
4024 r=
real(numi(j))/
real(im)
4043 SUBROUTINE minmxj(IM,JM,A,title)
4046 real A(IM,JM),rmin,rmax
4057 if(a(i,j).ge.rmax)rmax=a(i,j)
4058 if(a(i,j).le.rmin)rmin=a(i,j)
4061 write(6,150)rmin,rmax,title
4062 150
format(
'rmin=',e13.4,2x,
'rmax=',e13.4,2x,a8,
' ')
4078 SUBROUTINE mnmxja(IM,JM,A,imax,jmax,title)
4081 real A(IM,JM),rmin,rmax
4082 integer i,j,IM,JM,imax,jmax
4092 if(a(i,j).ge.rmax)
then 4097 if(a(i,j).le.rmin)rmin=a(i,j)
4100 write(6,150)rmin,rmax,title
4101 150
format(
'rmin=',e13.4,2x,
'rmax=',e13.4,2x,a8,
' ')
4140 SUBROUTINE spfft1(IMAX,INCW,INCG,KMAX,W,G,IDIR)
4142 INTEGER,
INTENT(IN):: IMAX,INCW,INCG,KMAX,IDIR
4143 COMPLEX,
INTENT(INOUT):: W(INCW,KMAX)
4144 REAL,
INTENT(INOUT):: G(INCG,KMAX)
4145 REAL:: AUX1(25000+INT(0.82*IMAX))
4146 REAL:: AUX2(20000+INT(0.57*IMAX))
4147 INTEGER:: NAUX1,NAUX2
4149 NAUX1=25000+int(0.82*imax)
4150 naux2=20000+int(0.57*imax)
4155 SELECT CASE(digits(1.))
4157 CALL scrft(1,w,incw,g,incg,imax,kmax,-1,1.,
4158 & aux1,naux1,aux2,naux2,0.,0)
4159 CALL scrft(0,w,incw,g,incg,imax,kmax,-1,1.,
4160 & aux1,naux1,aux2,naux2,0.,0)
4162 CALL dcrft(1,w,incw,g,incg,imax,kmax,-1,1.,
4163 & aux1,naux1,aux2,naux2,0.,0)
4164 CALL dcrft(0,w,incw,g,incg,imax,kmax,-1,1.,
4165 & aux1,naux1,aux2,naux2,0.,0)
4170 SELECT CASE(digits(1.))
4172 CALL srcft(1,g,incg,w,incw,imax,kmax,+1,1./imax,
4173 & aux1,naux1,aux2,naux2,0.,0)
4174 CALL srcft(0,g,incg,w,incw,imax,kmax,+1,1./imax,
4175 & aux1,naux1,aux2,naux2,0.,0)
4177 CALL drcft(1,g,incg,w,incw,imax,kmax,+1,1./imax,
4178 & aux1,naux1,aux2,naux2,0.,0)
4179 CALL drcft(0,g,incg,w,incw,imax,kmax,+1,1./imax,
4180 & aux1,naux1,aux2,naux2,0.,0)
4190 subroutine read_g(glob,ITOPO)
4193 integer*2 glob(360*120,180*120)
4198 parameter(ix=40*120,jx=50*120)
4199 parameter(ia=60*120,ja=30*120)
4201 integer*2 idat(ix,jx)
4206 real(kind=8) dloin,dlain,rlon,rlat
4208 open(235, file=
"./fort.235", access=
'direct', recl=43200*21600*2)
4213 call maxmin (glob,360*120*180*120,
'global0')
4219 rlon= -179.995833333333333333333333d0
4220 rlat= 89.995833333333333333333333d0
4241 subroutine maxmin(ia,len,tile)
4247 integer iaamax, iaamin, len, j, m, ja, kount
4248 integer(8) sum2,std,mean,isum
4249 integer i_count_notset,kount_9
4264 if ( ja .eq. -9999 )
then 4268 if ( ja .eq. -12345 ) i_count_notset=i_count_notset+1
4270 iaamax = max0( iaamax, ja )
4271 iaamin = min0( iaamin, ja )
4282 std = ifix(sqrt(float((sum2/(kount))-mean**2)))
4283 print*,tile,
' max=',iaamax,
' min=',iaamin,
' sum=',isum,
4284 &
' i_count_notset=',i_count_notset
4285 print*,tile,
' mean=',mean,
' std.dev=',std,
4286 &
' ko9s=',kount,kount_9,kount+kount_9
4299 SUBROUTINE minmaxj(IM,JM,A,title)
4302 real(kind=4) A(IM,JM),rmin,rmax,undef
4303 integer i,j,IM,JM,imax,jmax,imin,jmin,iundef
4304 character*8 title,chara
4320 if(a(i,j).ge.rmax)
then 4325 if(a(i,j).le.rmin)
then 4326 if ( a(i,j) .eq. undef )
then 4336 write(6,150)chara,rmin,imin,jmin,rmax,imax,jmax,iundef
4337 150
format(1x,a8,2x,
'rmin=',e13.4,2i6,2x,
'rmax=',e13.4,3i6)
4353 integer,
intent(in) :: siz
4354 real,
intent(in) :: lon(siz), lat(siz)
4355 real,
intent(out) :: x(siz), y(siz), z(siz)
4360 x(n) = cos(lat(n))*cos(lon(n))
4361 y(n) = cos(lat(n))*sin(lon(n))
4375 real,
parameter :: epsln30 = 1.e-30
4376 real,
parameter :: pi=3.1415926535897931
4377 real v1(3), v2(3), v3(3)
4380 real px, py, pz, qx, qy, qz, ddd;
4383 px = v1(2)*v2(3) - v1(3)*v2(2)
4384 py = v1(3)*v2(1) - v1(1)*v2(3)
4385 pz = v1(1)*v2(2) - v1(2)*v2(1)
4387 qx = v1(2)*v3(3) - v1(3)*v3(2);
4388 qy = v1(3)*v3(1) - v1(1)*v3(3);
4389 qz = v1(1)*v3(2) - v1(2)*v3(1);
4391 ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
4392 if ( ddd <= 0.0 )
then 4395 ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
4396 if( abs(ddd-1) < epsln30 ) ddd = 1;
4397 if( abs(ddd+1) < epsln30 ) ddd = -1;
4398 if ( ddd>1. .or. ddd<-1. )
then 4425 real,
parameter :: epsln10 = 1.e-10
4426 real,
parameter :: epsln8 = 1.e-8
4427 real,
parameter :: pi=3.1415926535897931
4428 real,
parameter :: range_check_criteria=0.05
4433 real lon2(npts), lat2(npts)
4434 real x2(npts), y2(npts), z2(npts)
4435 real lon1_1d(1), lat1_1d(1)
4436 real x1(1), y1(1), z1(1)
4437 real pnt0(3),pnt1(3),pnt2(3)
4439 real max_x2,min_x2,max_y2,min_y2,max_z2,min_z2
4441 call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
4444 call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
4447 if( x1(1) > max_x2+range_check_criteria )
return 4449 if( x1(1)+range_check_criteria < min_x2 )
return 4451 if( y1(1) > max_y2+range_check_criteria )
return 4453 if( y1(1)+range_check_criteria < min_y2 )
return 4455 if( z1(1) > max_z2+range_check_criteria )
return 4457 if( z1(1)+range_check_criteria < min_z2 )
return 4465 if(abs(x1(1)-x2(i)) < epsln10 .and.
4466 & abs(y1(1)-y2(i)) < epsln10 .and.
4467 & abs(z1(1)-z2(i)) < epsln10 )
then 4472 if(ip1>npts) ip1 = 1
4482 anglesum = anglesum + angle
4485 if(abs(anglesum-2*pi) < epsln8)
then 4518 function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
4519 & glat,zavg,zslm,delxn)
4524 real lon1,lat1,lon2,lat2,oro,delxn
4527 integer zavg(imn,jmn),zslm(imn,jmn)
4528 integer i, j, ist, ien, jst, jen, i1
4530 real xland,xwatr,xl1,xs1,slm,xnsum
4533 if( glat(j) .GT. lat1 )
then 4539 if( glat(j) .GT. lat2 )
then 4546 ist = lon1/delxn + 1
4548 if(ist .le.0) ist = ist + imn
4549 if(ien < ist) ien = ien + imn
4560 do i1 = 1, ien - ist + 1
4562 if( i .LE. 0) i = i + imn
4563 if( i .GT. imn) i = i - imn
4564 xland = xland + float(zslm(i,j))
4565 xwatr = xwatr + float(1-zslm(i,j))
4567 height = float(zavg(i,j))
4568 IF(height.LT.-990.) height = 0.0
4569 xl1 = xl1 + height * float(zslm(i,j))
4570 xs1 = xs1 + height * float(1-zslm(i,j))
4573 if( xnsum > 1.)
THEN 4574 slm = float(nint(xland/xnsum))
4586 if( i .LE. 0) i = i + imn
4587 if( i .GT. imn) i = i - imn
4588 height = float(zavg(i,j))
4589 IF(height.LT.-990.) height = 0.0
4625 subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
4626 & glat,zavg,zslm,delxn,xnsum1,xnsum2,HC)
4629 real,
intent(out) :: xnsum1,xnsum2,HC
4631 real lon1,lat1,lon2,lat2,oro,delxn
4634 integer zavg(IMN,JMN),zslm(IMN,JMN)
4635 integer i, j, ist, ien, jst, jen, i1
4637 real XW1,XW2,slm,xnsum
4640 if( glat(j) .GT. lat1 )
then 4646 if( glat(j) .GT. lat2 )
then 4653 ist = lon1/delxn + 1
4655 if(ist .le.0) ist = ist + imn
4656 if(ien < ist) ien = ien + imn
4664 do i1 = 1, ien - ist + 1
4666 if( i .LE. 0) i = i + imn
4667 if( i .GT. imn) i = i - imn
4669 height = float(zavg(i,j))
4670 IF(height.LT.-990.) height = 0.0
4672 xw2 = xw2 + height ** 2
4675 var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
4676 hc = 1116.2 - 0.878 * var
4682 if( i .LE. 0) i = i + imn
4683 if( i .GT. imn) i = i - imn
4684 height = float(zavg(i,j))
4685 IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4721 subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
4722 & glat,zavg,zslm,delxn,xnsum1,xnsum2,HC)
4725 real,
intent(out) :: xnsum1,xnsum2
4726 real lon1,lat1,lon2,lat2,oro,delxn
4729 integer zavg(IMN,JMN),zslm(IMN,JMN)
4730 integer i, j, ist, ien, jst, jen, i1
4732 real XW1,XW2,slm,xnsum
4738 if( glat(j) .GT. lat1 )
then 4744 if( glat(j) .GT. lat2 )
then 4751 ist = lon1/delxn + 1
4753 if(ist .le.0) ist = ist + imn
4754 if(ien < ist) ien = ien + imn
4762 if( i .LE. 0) i = i + imn
4763 if( i .GT. imn) i = i - imn
4764 height = float(zavg(i,j))
4765 IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4785 subroutine nanc(a,l,c)
4786 integer inan1,inan2,inan3,inan4,inaq1,inaq2,inaq3,inaq4
4789 equivalence (itest,word)
4792 data inan1/x
'7F800001'/
4793 data inan2/x
'7FBFFFFF'/
4794 data inan3/x
'FF800001'/
4795 data inan4/x
'FFBFFFFF'/
4799 data inaq1/x
'7FC00000'/
4800 data inaq2/x
'7FFFFFFF'/
4801 data inaq3/x
'FFC00000'/
4802 data inaq4/x
'FFFFFFFF'/
4804 real(kind=8)a(l),rtc,t1,t2
4811 if( (itest .GE. inan1 .AND. itest .LE. inan2) .OR.
4812 * (itest .GE. inan3 .AND. itest .LE. inan4) )
then 4813 print *,
' NaNs detected at word',k,
' ',c
4816 if( (itest .GE. inaq1 .AND. itest .LE. inaq2) .OR.
4817 * (itest .GE. inaq3 .AND. itest .LE. inaq4) )
then 4818 print *,
' NaNq detected at word',k,
' ',c
4826 102
format(
' time to check ',i9,
' words is ',f10.4,
' ',a24)
4834 real function timef()
4835 character(8) :: date
4836 character(10) :: time
4837 character(5) :: zone
4838 integer,
dimension(8) :: values
4841 call date_and_time(date,time,zone,values)
4842 total=(3600*values(5))+(60*values(6))
4844 elapsed=float(total) + (1.0e-3*float(values(8)))
subroutine get_index(IMN, JMN, npts, lonO, latO, DELXN, jst, jen, ilist, numx)
Determine the location of a cubed-sphere point within the high-resolution orography data...
real function timef()
Get the date/time for the system clock.
subroutine netcdf_err(err, string)
Check NetCDF error code and output the error message.
real function get_lat_angle(dy, DEGRAD)
Convert the 'y' direction distance of a cubed-sphere grid point to the corresponding distance in lati...
subroutine mnmxja(IM, JM, A, imax, jmax, title)
Print out the maximum and minimum values of an array.
subroutine read_g(glob, ITOPO)
Read input global 30-arc second orography data.
real function spherical_angle(v1, v2, v3)
Compute spherical angle.
subroutine minmaxj(IM, JM, A, title)
Print out the maximum and minimum values of an array and their i/j location.
subroutine get_xnsum3(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, zslm, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
subroutine get_xnsum2(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, zslm, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
subroutine rg2gg(im, jm, numi, a)
Convert from a reduced grid to a full grid.
subroutine write_netcdf(im, jm, slm, land_frac, oro, orf, hprime, ntiles, tile, geolon, geolat, lon, lat)
Write out orography file in netcdf format.
subroutine gg2rg(im, jm, numi, a)
Convert from a full grid to a reduced grid.
logical function inside_a_polygon(lon1, lat1, npts, lon2, lat2)
Check if a point is inside a polygon.
real function spherical_distance(theta1, phi1, theta2, phi2)
Compute a great circle distance between two points.
subroutine nanc(a, l, c)
Report NaNS and NaNQ within an address range for 8-byte real words.
subroutine get_mismatch_index(im_in, jm_in, geolon_in, geolat_in, bitmap_in, num_out, lon_out, lat_out, iindx, jindx)
For unmapped land points, find the nearest land point on the input data and pass back its i/j index...
subroutine maxmin(ia, len, tile)
Print the maximum, mininum, mean and standard deviation of an array.
subroutine interpolate_mismatch(im_in, jm_in, data_in, num_out, data_out, iindx, jindx)
Replace unmapped model land points with the nearest land point on the input grid. ...
real function get_lon_angle(dx, lat, DEGRAD)
Convert the 'x' direction distance of a cubed-sphere grid point to the corresponding distance in long...
subroutine latlon2xyz(siz, lon, lat, x, y, z)
Convert from latitude and longitude to x,y,z coordinates.
subroutine minmxj(IM, JM, A, title)
Print out the maximum and minimum values of an array.
real function get_xnsum(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, zslm, delxn)
Count the number of high-resolution orography points that are higher than the model grid box average ...