orog_mask_tools  1.13.0
mtnlm7_oclsm.F
Go to the documentation of this file.
1 C> @file
2 C> Terrain maker for global spectral model.
3 C> @author Mark Iredell @date 92-04-16
4 
5 C> This program creates 7 terrain-related files computed from the
6 C> GMTED2010 terrain dataset. The model physics grid parameters and
7 C> spectral truncation and filter parameters are read by this program as
8 C> input.
9 C>
10 C> The 7 files produced are:
11 C> 1. sea-land mask on model physics grid
12 C> 2. gridded orography on model physics grid
13 C> 3. mountain std dev on model physics grid
14 C> 4. spectral orography in spectral domain
15 C> 5. unfiltered gridded orography on model physics grid
16 C> 6. grib sea-land mask on model physics grid
17 C> 7. grib gridded orography on model physics grid
18 C>
19 C> The orography is only filtered for wavenumbers greater than nf0. For
20 C> wavenumbers n between nf0 and nf1, the orography is filtered by the
21 C> factor 1-((n-nf0)/(nf1-nf0))**2. The filtered orography will not have
22 C> information beyond wavenumber nf1.
23 C>
24 C> PROGRAM HISTORY LOG:
25 C> - 92-04-16 IREDELL
26 C> - 98-02-02 IREDELL FILTER
27 C> - 98-05-31 HONG Modified for subgrid orography used in Kim's scheme
28 C> - 98-12-31 HONG Modified for high-resolution GTOPO orography
29 C> - 99-05-31 HONG Modified for getting OL4 (mountain fraction)
30 C> - 00-02-10 Moorthi's modifications
31 C> - 00-04-11 HONG Modified for reduced grids
32 C> - 00-04-12 Iredell Modified for reduced grids
33 C> - 02-01-07 (*j*) modified for principal axes of orography
34 C> There are now 14 files, 4 additional for lm mb
35 C> - 04-04-04 (*j*) re-Test on IST/ilen calc for sea-land mask(*j*)
36 C> - 04-09-04 minus sign here in MAKEOA IST and IEN as in MAKEMT!
37 C> - 05-09-05 if test on HK and HLPRIM for GAMMA SQRT
38 C> - 07-08-07 replace 8' with 30" incl GICE, conintue w/ S-Y. lake slm
39 C> - 08-08-07 All input 30", UMD option, and filter as described below
40 C> Quadratic filter applied by default.
41 C> NF0 is normally set to an even value beyond the previous truncation,
42 C> for example, for jcap=382, NF0=254+2
43 C> NF1 is set as jcap+2 (and/or nearest even), eg., for t382, NF1=382+2=384
44 C> if no filter is desired then NF1=NF0=0 and ORF=ORO
45 C> but if no filter but spectral to grid (with gibbs) then NF1=jcap+2, and NF1=jcap+1
46 C>
47 C> INPUT FILES:
48 C> - UNIT5 - PHYSICS LONGITUDES (IM), PHYSICS LATITUDES (JM),
49 C> SPECTRAL TRUNCATION (NM), RHOMBOIDAL FLAG (NR),
50 C> AND FIRST AND SECOND FILTER PARAMETERS (NF0,NF1).
51 C> RESPECTIVELY READ IN FREE FORMAT.
52 C> - NCID - GMTED2010 USGS orography (NetCDF)
53 C> - NCID - 30" UMD land cover mask. (NetCDF)
54 C> - NCID - GICE Grumbine 30" RAMP Antarctica orog IMNx3601. (NetCDF)
55 C> - UNIT25 - Ocean land-sea mask on gaussian grid
56 C>
57 C> OUTPUT FILES:
58 C> - UNIT51 - SEA-LAND MASK (IM,JM)
59 C> - UNIT52 - GRIDDED OROGRAPHY (IM,JM)
60 C> - UNIT54 - SPECTRAL OROGRAPHY ((NM+1)*((NR+1)*NM+2))
61 C> - UNIT55 - UNFILTERED GRIDDED OROGRAPHY (IM,JM)
62 C> - UNIT57 - GRIB GRIDDED OROGRAPHY (IM,JM)
63 C>
64 C> SUBPROGRAMS CALLED:
65 C> - UNIQUE:
66 C> - TERSUB - MAIN SUBPROGRAM
67 C> - SPLAT - COMPUTE GAUSSIAN LATITUDES OR EQUALLY-SPACED LATITUDES
68 C> - LIBRARY:
69 C> - SPTEZ - SPHERICAL TRANSFORM
70 C> - GBYTES - UNPACK BITS
71 C>
72 C> @return 0 for success, error code otherwise.
73  include 'netcdf.inc'
74  logical fexist, opened
75  integer fsize, ncid, error, id_dim, nx, ny
76  character(len=256) :: outgrid = "none"
77  character(len=256) :: inputorog = "none"
78  character(len=256) :: merge_file = "none"
79  logical :: mask_only = .false.
80  integer :: mtnres,im,jm,nm,nr,nf0,nf1,efac,nw
81  fsize=65536
82  READ(5,*) outgrid
83  READ(5,*) mask_only
84  READ(5,*) merge_file
85  nm=0
86  nf0=0
87  nf1=0
88  efac=0
89  nr=0
90  print*, "INPUTOROG= ", trim(inputorog)
91  print*, "MASK_ONLY", mask_only
92  print*, "MERGE_FILE ", trim(merge_file)
93 ! --- MTNRES defines the input (highest) elev resolution
94 ! --- =1 is topo30 30" in units of 1/2 minute.
95 ! so MTNRES for old values must be *2.
96 ! =16 is now Song Yu's 8' orog the old ops standard
97 ! --- other possibilities are =8 for 4' and =4 for 2' see
98 ! HJ for T1000 test. Must set to 1 for now.
99  mtnres=1
100  print*, mtnres,nm,nr,nf0,nf1,efac
101  nw=(nm+1)*((nr+1)*nm+2)
102  imn = 360*120/mtnres
103  jmn = 180*120/mtnres
104  print *, ' Starting terr12 mtnlm7_slm30.f IMN,JMN:',imn,jmn
105 
106 ! --- read the grid resolution from OUTGRID.
107  inquire(file=trim(outgrid), exist=fexist)
108  if(.not. fexist) then
109  print*, "FATAL ERROR: file "//trim(outgrid)
110  print*, " does not exist."
111  CALL errexit(4)
112  endif
113  do ncid = 103, 512
114  inquire( ncid,opened=opened )
115  if( .NOT.opened )exit
116  end do
117 
118  print*, "READ outgrid=", trim(outgrid)
119  error=nf__open(trim(outgrid),nf_nowrite,fsize,ncid)
120  call netcdf_err(error, 'Open file '//trim(outgrid) )
121  error=nf_inq_dimid(ncid, 'nx', id_dim)
122  call netcdf_err(error, 'inquire dimension nx from file '//
123  & trim(outgrid) )
124  error=nf_inq_dimlen(ncid,id_dim,nx)
125  call netcdf_err(error, 'inquire dimension nx length '//
126  & 'from file '//trim(outgrid) )
127 
128  error=nf_inq_dimid(ncid, 'ny', id_dim)
129  call netcdf_err(error, 'inquire dimension ny from file '//
130  & trim(outgrid) )
131  error=nf_inq_dimlen(ncid,id_dim,ny)
132  call netcdf_err(error, 'inquire dimension ny length '//
133  & 'from file '//trim(outgrid) )
134  im = nx/2
135  jm = ny/2
136  print*, "nx, ny, im, jm = ", nx, ny, im, jm
137  error=nf_close(ncid)
138  call netcdf_err(error, 'close file '//trim(outgrid) )
139 
140  CALL tersub(imn,jmn,im,jm,nm,nr,nf0,nf1,nw,efac,
141  & outgrid,inputorog,mask_only,merge_file)
142  stop
143  END
144 
165  SUBROUTINE tersub(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,
166  & OUTGRID,INPUTOROG,MASK_ONLY,MERGE_FILE)
167  implicit none
168  include 'netcdf.inc'
169 C
170  integer :: IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW
171  character(len=*), intent(in) :: OUTGRID
172  character(len=*), intent(in) :: INPUTOROG
173  character(len=*), intent(in) :: MERGE_FILE
174 
175  logical, intent(in) :: mask_only
176 
177  real, parameter :: MISSING_VALUE=-9999.
178  real, PARAMETER :: PI=3.1415926535897931
179  integer, PARAMETER :: NMT=14
180 
181  integer :: efac,zsave1,zsave2
182  integer :: mskocn,notocn
183  integer :: i,j,nx,ny,ncid,js,jn,iw,ie,k,it,jt,error,id_dim
184  integer :: id_var,nx_in,ny_in,fsize,wgta,IN,INW,INE,IS,ISW,ISE
185  integer :: M,N,ios,istat,itest,jtest
186  integer :: i_south_pole,j_south_pole,i_north_pole,j_north_pole
187  integer :: maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
188  integer(1) :: i3save
189  integer(2) :: i2save
190 
191  integer, allocatable :: JST(:),JEN(:),numi(:)
192 
193  integer, allocatable :: IST(:,:),IEN(:,:),ZSLMX(:,:)
194  integer, allocatable :: ZAVG(:,:),ZSLM(:,:)
195  integer(1), allocatable :: UMD(:,:)
196  integer(2), allocatable :: glob(:,:)
197 
198  integer, allocatable :: IWORK(:,:,:)
199 
200  real :: DEGRAD,maxlat, minlat,timef,tbeg,tend,tbeg1
201  real :: PHI,DELXN,slma,oroa,vara,var4a,xn,XS,FFF,WWW
202  real :: sumdif,avedif
203 
204  real, allocatable :: COSCLT(:),WGTCLT(:),RCLT(:),XLAT(:),DIFFX(:)
205  real, allocatable :: XLON(:),ORS(:),oaa(:),ola(:),GLAT(:)
206 
207  real, allocatable :: GEOLON(:,:),GEOLON_C(:,:),DX(:,:)
208  real, allocatable :: GEOLAT(:,:),GEOLAT_C(:,:),DY(:,:)
209  real, allocatable :: SLM(:,:),ORO(:,:),VAR(:,:),ORF(:,:)
210  real, allocatable :: land_frac(:,:),lake_frac(:,:)
211  real, allocatable :: THETA(:,:),GAMMA(:,:),SIGMA(:,:),ELVMAX(:,:)
212  real, allocatable :: VAR4(:,:),SLMI(:,:)
213  real, allocatable :: WORK1(:,:),WORK2(:,:),WORK3(:,:),WORK4(:,:)
214  real, allocatable :: WORK5(:,:),WORK6(:,:)
215  real, allocatable :: tmpvar(:,:)
216  real, allocatable :: slm_in(:,:), lon_in(:,:), lat_in(:,:)
217  real(4), allocatable:: GICE(:,:),OCLSM(:,:)
218 
219  real, allocatable :: OA(:,:,:),OL(:,:,:),HPRIME(:,:,:)
220  real, allocatable :: oa_in(:,:,:), ol_in(:,:,:)
221 
222  logical :: grid_from_file,fexist,opened
223  logical :: SPECTR, FILTER
224  logical :: is_south_pole(IM,JM), is_north_pole(IM,JM)
225 
226  tbeg1=timef()
227  tbeg=timef()
228  fsize = 65536
229 ! integers
230  allocate (jst(jm),jen(jm),numi(jm))
231  allocate (ist(im,jm),ien(im,jm),zslmx(2700,1350))
232  allocate (glob(imn,jmn))
233 
234 ! reals
235  allocate (cosclt(jm),wgtclt(jm),rclt(jm),xlat(jm),diffx(jm/2))
236  allocate (xlon(im),ors(nw),oaa(4),ola(4),glat(jmn))
237 
238  allocate (zavg(imn,jmn))
239  allocate (zslm(imn,jmn))
240  allocate (umd(imn,jmn))
241 
242 !
243 ! SET CONSTANTS AND ZERO FIELDS
244 !
245  degrad = 180./pi
246  spectr = nm .GT. 0 ! if NM <=0 grid is assumed lat/lon
247  filter = .true. ! Spectr Filter defaults true and set by NF1 & NF0
248  mskocn = 1 ! Ocean land sea mask =1, =0 if not present
249  notocn = 1 ! =1 Ocean lsm input reverse: Ocean=1, land=0
250 ! --- The LSM Gaussian file from the ocean model sometimes arrives with
251 ! --- 0=Ocean and 1=Land or it arrives with 1=Ocean and 0=land without
252 ! --- metadata to distinguish its disposition. The AI below mitigates this.
253 
254  print *,' In TERSUB'
255  if (mskocn .eq. 1)then
256  print *,' Ocean Model LSM Present and '
257  print *, ' Overrides OCEAN POINTS in LSM: mskocn=',mskocn
258  if (notocn .eq. 1)then
259  print *,' Ocean LSM Reversed: NOTOCN=',notocn
260  endif
261  endif
262 
263  print *,' Attempt to open/read UMD 30sec slmsk.'
264 
265  error=nf__open("./landcover.umd.30s.nc",nf_nowrite,fsize,ncid)
266  call netcdf_err(error, 'Open file landcover.umd.30s.nc' )
267  error=nf_inq_varid(ncid, 'land_mask', id_var)
268  call netcdf_err(error, 'Inquire varid of land_mask')
269  error=nf_get_var_int1(ncid, id_var, umd)
270  call netcdf_err(error, 'Inquire data of land_mask')
271  error = nf_close(ncid)
272 
273  print *,' UMD lake, UMD(50,50)=',umd(50,50)
274 C
275 C- READ_G for global 30" terrain
276 C
277  print *,' Call read_g to read global topography'
278  call read_g(glob)
279 ! --- transpose even though glob 30" is from S to N and NCEP std is N to S
280  do j=1,jmn/2
281  do i=1,imn
282  jt=jmn - j + 1
283  i2save = glob(i,j)
284  glob(i,j)=glob(i,jt)
285  glob(i,jt) = i2save
286  enddo
287  enddo
288 ! --- transpose glob as USGS 30" is from dateline and NCEP std is 0
289  do j=1,jmn
290  do i=1,imn/2
291  it=imn/2 + i
292  i2save = glob(i,j)
293  glob(i,j)=glob(it,j)
294  glob(it,j) = i2save
295  enddo
296  enddo
297  print *,' After read_g, glob(500,500)=',glob(500,500)
298 !
299 
300 ! --- IMN,JMN
301  print*, ' IM, JM, NM, NR, NF0, NF1, EFAC'
302  print*, im,jm,nm,nr,nf0,nf1,efac
303  print *,' imn,jmn,glob(imn,jmn)=',imn,jmn,glob(imn,jmn)
304  print *,' UBOUND ZAVG=',ubound(zavg)
305  print *,' UBOUND glob=',ubound(glob)
306  print *,' UBOUND ZSLM=',ubound(zslm)
307  print *,' UBOUND GICE=',imn+1,3601
308  print *,' UBOUND OCLSM=',im,jm
309 !
310 ! --- 0 is ocean and 1 is land for slm
311 !
312 C
313 ! --- ZSLM initialize with all land 1, ocean 0
314  zslm=1
315 ! --- ZAVG initialize from glob
316  zavg=glob
317 
318 ! --- transpose mask even though glob 30" is from N to S and NCEP std is S to N
319  do j=1,jmn/2
320  do i=1,imn
321  jt=jmn - j + 1
322  i3save = umd(i,j)
323  umd(i,j)=umd(i,jt)
324  umd(i,jt) = i3save
325  enddo
326  enddo
327 ! --- transpose UMD as USGS 30" is from dateline and NCEP std is 0
328  do j=1,jmn
329  do i=1,imn/2
330  it=imn/2 + i
331  i3save = umd(i,j)
332  umd(i,j)=umd(it,j)
333  umd(it,j) = i3save
334  enddo
335  enddo
336 ! --- Non-land is 0.
337  do j=1,jmn
338  do i=1,imn
339  if ( umd(i,j) .eq. 0 ) zslm(i,j) = 0
340  enddo
341  enddo
342 
343  deallocate (zslmx,umd,glob)
344 ! ---
345 ! --- Fixing an error in the topo 30" data set at pole (-9999).
346  do i=1,imn
347  zslm(i,1)=0
348  zslm(i,jmn)=1
349  enddo
350 !
351 ! print *,' kount1,2,ZAVG(1,1),ZAVG(imn,jmn),ZAVG(500,500)',
352 ! & kount,kount2,ZAVG(1,1),ZAVG(imn,jmn),ZAVG(500,500)
353 ! --- The center of pixel (1,1) is 89.9958333N/179.9958333W with dx/dy
354 ! --- spacing of 1/120 degrees.
355 !
356 ! When the gaussian grid routines makemt, makepc and makeoa are
357 ! removed, numi can be removed.
358  do j=1,jm
359  numi(j)=im
360  enddo
361 !
362 ! This code assumes that lat runs from north to south for gg!
363 !
364 
365  print *,' SPECTR=',spectr,' ** with GICE-07 **'
366  IF (spectr) THEN
367  CALL splat(4,jm,cosclt,wgtclt)
368  DO j=1,jm/2
369  rclt(j) = acos(cosclt(j))
370  ENDDO
371  DO j = 1,jm/2
372  phi = rclt(j) * degrad
373  xlat(j) = 90. - phi
374  xlat(jm-j+1) = phi - 90.
375  ENDDO
376  ELSE
377  CALL splat(0,jm,cosclt,wgtclt)
378  DO j=1,jm
379  rclt(j) = acos(cosclt(j))
380  xlat(j) = 90.0 - rclt(j) * degrad
381  ENDDO
382  ENDIF
383 
384  allocate (gice(imn+1,3601))
385 !
386  sumdif = 0.
387  DO j = jm/2,2,-1
388  diffx(j) = xlat(j) - xlat(j-1)
389  sumdif = sumdif + diffx(j)
390  ENDDO
391  avedif=sumdif/(float(jm/2))
392 ! print *,' XLAT= avedif: ',avedif
393 ! write (6,107) (xlat(J)-xlat(j-1),J=JM,2,-1)
394  print *,' XLAT='
395  write (6,106) (xlat(j),j=jm,1,-1)
396  106 format( 10(f7.3,1x))
397  107 format( 10(f9.5,1x))
398 C
399  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
400 C
401  DO j=1,jmn
402  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
403  ENDDO
404  print *,
405  & ' Before GICE ZAVG(1,2)=',zavg(1,2),zslm(1,2)
406  print *,
407  & ' Before GICE ZAVG(1,12)=',zavg(1,12),zslm(1,12)
408  print *,
409  & ' Before GICE ZAVG(1,52)=',zavg(1,52),zslm(1,52)
410  print *,
411  & ' Before GICE ZAVG(1,112)=',zavg(1,jmn-112),zslm(1,112)
412 
413 ! Read 30-sec Antarctica RAMP data. Points scan from South
414 ! to North, and from Greenwich to Greenwich.
415 
416  error=nf__open("./topography.antarctica.ramp.30s.nc",
417  & nf_nowrite,fsize,ncid)
418  call netcdf_err(error, 'Opening RAMP topo file' )
419  error=nf_inq_varid(ncid, 'topo', id_var)
420  call netcdf_err(error, 'Inquire varid of RAMP topo')
421  error=nf_get_var_real(ncid, id_var, gice)
422  call netcdf_err(error, 'Inquire data of RAMP topo')
423  error = nf_close(ncid)
424 
425  print *,' GICE 30" Antarctica RAMP orog 43201x3601 read OK'
426  print *,' Processing! '
427  print *,' Processing! '
428  print *,' Processing! '
429  do j = 1, 3601
430  do i = 1, imn
431  zsave1 = zavg(i,j)
432  zsave2 = zslm(i,j)
433  if( gice(i,j) .ne. -99. .and. gice(i,j) .ne. -1.0 ) then
434  if ( gice(i,j) .gt. 0.) then
435  zavg(i,j) = int( gice(i,j) + 0.5 )
436 !! --- for GICE values less than or equal to 0 (0, -1, or -99) then
437 !! --- radar-sat (RAMP) values are not valid and revert back to old orog
438  zslm(i,j) = 1
439  endif
440  endif
441  152 format(1x,' ZAVG(i=',i4,' j=',i4,')=',i5,i3,
442  &' orig:',i5,i4,' Lat=',f7.3,f8.2,'E',' GICE=',f8.1)
443  enddo
444  enddo
445 
446  deallocate (gice)
447 
448  allocate (oclsm(im,jm),slmi(im,jm))
449 !C
450 C COMPUTE MOUNTAIN DATA : ORO SLM VAR (Std Dev) OC
451 C
452 ! --- The coupled ocean model is already on a Guasian grid if (IM,JM)
453 ! --- Attempt to Open the file if mskocn=1
454  istat=0
455  if (mskocn .eq. 1) then
456 ! open(25,form='unformatted',iostat=istat)
457 ! open(25,form='binary',iostat=istat)
458 ! --- open to fort.25 with link to file in script
459  open(25,form='formatted',iostat=istat)
460  if (istat.ne.0) then
461  mskocn = 0
462  print *,' Ocean lsm file Open failure: mskocn,istat=',mskocn,istat
463  else
464  mskocn = 1
465  print *,' Ocean lsm file Opened OK: mskocn,istat=',mskocn,istat
466  endif
467 ! --- Read it in
468  ios=0
469  oclsm=0.
470 ! read(25,iostat=ios)OCLSM
471  read(25,*,iostat=ios)oclsm
472  if (ios.ne.0) then
473  mskocn = 0
474 ! --- did not properly read Gaussian grid ocean land-sea mask, but
475 ! continue using ZSLMX
476  print *,' Rd fail: Ocean lsm - continue, mskocn,ios=',mskocn,ios
477  else
478  mskocn = 1
479  print *,' Rd OK: ocean lsm: mskocn,ios=',mskocn,ios
480 ! --- LSM initialized to ocean mask especially for case where Ocean
481 ! --- changed by ocean model to land to cope with its problems
482 ! --- remember, that lake mask is in zslm to be assigned in MAKEMT.
483  if ( mskocn .eq. 1 ) then
484  DO j = 1,jm
485  DO i = 1,im
486  if ( notocn .eq. 0 ) then
487  slmi(i,j) = float(nint(oclsm(i,j)))
488  else
489  if ( nint(oclsm(i,j)) .eq. 0) then
490  slmi(i,j) = 1
491  else
492  slmi(i,j) = 0
493  endif
494  endif
495  enddo
496  enddo
497  print *,' OCLSM',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
498  print *,' SLMI:',slmi(1,1),slmi(50,50),slmi(75,75),slmi(im,jm)
499 ! --- Diag
500 ! WRITE(27,iostat=ios) REAL(SLMI,4)
501 ! print *,' write SLMI/OCLSM diag input:',ios
502  endif
503  endif
504 
505  else
506  print *,' Not using Ocean model land sea mask'
507  endif
508 
509  if (mskocn .eq. 1)then
510  print *,' LSM:',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
511  endif
512 
513  allocate (geolon(im,jm),geolon_c(im+1,jm+1),dx(im,jm))
514  allocate (geolat(im,jm),geolat_c(im+1,jm+1),dy(im,jm))
515  allocate (slm(im,jm),oro(im,jm),var(im,jm),var4(im,jm))
516  allocate (land_frac(im,jm),lake_frac(im,jm))
517 
518 !--- reading grid file.
519  grid_from_file = .false.
520  is_south_pole = .false.
521  is_north_pole = .false.
522  i_south_pole = 0
523  j_south_pole = 0
524  i_north_pole = 0
525  j_north_pole = 0
526  if( trim(outgrid) .NE. "none" ) then
527  grid_from_file = .true.
528  inquire(file=trim(outgrid), exist=fexist)
529  if(.not. fexist) then
530  print*, "FATAL ERROR: file "//trim(outgrid)
531  print*, "does not exist."
532  CALL errexit(4)
533  endif
534  do ncid = 103, 512
535  inquire( ncid,opened=opened )
536  if( .NOT.opened )exit
537  end do
538 
539  print*, "outgrid=", trim(outgrid)
540  error=nf__open(trim(outgrid),nf_nowrite,fsize,ncid)
541  call netcdf_err(error, 'Open file '//trim(outgrid) )
542  error=nf_inq_dimid(ncid, 'nx', id_dim)
543  call netcdf_err(error, 'inquire dimension nx from file '//
544  & trim(outgrid) )
545  nx = 2*im
546  ny = 2*jm
547  print*, "Read the grid from file "//trim(outgrid)
548 
549  allocate(tmpvar(nx+1,ny+1))
550 
551  error=nf_inq_varid(ncid, 'x', id_var)
552  call netcdf_err(error, 'inquire varid of x from file '
553  & //trim(outgrid) )
554  error=nf_get_var_double(ncid, id_var, tmpvar)
555  call netcdf_err(error, 'inquire data of x from file '
556  & //trim(outgrid) )
557  !--- adjust lontitude to be between 0 and 360.
558  do j = 1,ny+1; do i = 1,nx+1
559  if(tmpvar(i,j) .NE. missing_value) then
560  if(tmpvar(i,j) .GT. 360) tmpvar(i,j) = tmpvar(i,j) - 360
561  if(tmpvar(i,j) .LT. 0) tmpvar(i,j) = tmpvar(i,j) + 360
562  endif
563  enddo; enddo
564 
565  geolon(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
566  geolon_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
567 
568  error=nf_inq_varid(ncid, 'y', id_var)
569  call netcdf_err(error, 'inquire varid of y from file '
570  & //trim(outgrid) )
571  error=nf_get_var_double(ncid, id_var, tmpvar)
572  call netcdf_err(error, 'inquire data of y from file '
573  & //trim(outgrid) )
574  geolat(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
575  geolat_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
576 
577  !--- figure out pole location.
578  maxlat = -90
579  minlat = 90
580  i_north_pole = 0
581  j_north_pole = 0
582  i_south_pole = 0
583  j_south_pole = 0
584  do j = 1, ny+1; do i = 1, nx+1
585  if( tmpvar(i,j) > maxlat ) then
586  i_north_pole=i
587  j_north_pole=j
588  maxlat = tmpvar(i,j)
589  endif
590  if( tmpvar(i,j) < minlat ) then
591  i_south_pole=i
592  j_south_pole=j
593  minlat = tmpvar(i,j)
594  endif
595  enddo ; enddo
596  !--- only when maxlat is close to 90. the point is north pole
597  if(maxlat < 89.9 ) then
598  i_north_pole = 0
599  j_north_pole = 0
600  endif
601  if(minlat > -89.9 ) then
602  i_south_pole = 0
603  j_south_pole = 0
604  endif
605  print*, "minlat=", minlat, "maxlat=", maxlat
606  print*, "north pole supergrid index is ",
607  & i_north_pole, j_north_pole
608  print*, "south pole supergrid index is ",
609  & i_south_pole, j_south_pole
610  deallocate(tmpvar)
611 
612  if(i_south_pole >0 .and. j_south_pole > 0) then
613  if(mod(i_south_pole,2)==0) then ! stretched grid
614  do j = 1, jm; do i = 1, im
615  if(i==i_south_pole/2 .and. (j==j_south_pole/2
616  & .or. j==j_south_pole/2+1) ) then
617  is_south_pole(i,j) = .true.
618  print*, "south pole at i,j=", i, j
619  endif
620  enddo; enddo
621  else
622  do j = 1, jm; do i = 1, im
623  if((i==i_south_pole/2 .or. i==i_south_pole/2+1)
624  & .and. (j==j_south_pole/2 .or.
625  & j==j_south_pole/2+1) ) then
626  is_south_pole(i,j) = .true.
627  print*, "south pole at i,j=", i, j
628  endif
629  enddo; enddo
630  endif
631  endif
632  if(i_north_pole >0 .and. j_north_pole > 0) then
633  if(mod(i_north_pole,2)==0) then ! stretched grid
634  do j = 1, jm; do i = 1, im
635  if(i==i_north_pole/2 .and. (j==j_north_pole/2 .or.
636  & j==j_north_pole/2+1) ) then
637  is_north_pole(i,j) = .true.
638  print*, "north pole at i,j=", i, j
639  endif
640  enddo; enddo
641  else
642  do j = 1, jm; do i = 1, im
643  if((i==i_north_pole/2 .or. i==i_north_pole/2+1)
644  & .and. (j==j_north_pole/2 .or.
645  & j==j_north_pole/2+1) ) then
646  is_north_pole(i,j) = .true.
647  print*, "north pole at i,j=", i, j
648  endif
649  enddo; enddo
650  endif
651  endif
652 
653 
654  allocate(tmpvar(nx,ny))
655  error=nf_inq_varid(ncid, 'area', id_var)
656  call netcdf_err(error, 'inquire varid of area from file '
657  & //trim(outgrid) )
658  error=nf_get_var_double(ncid, id_var, tmpvar)
659  call netcdf_err(error, 'inquire data of area from file '
660  & //trim(outgrid) )
661 
662  do j = 1, jm
663  do i = 1, im
664  dx(i,j) = sqrt(tmpvar(2*i-1,2*j-1)+tmpvar(2*i,2*j-1)
665  & +tmpvar(2*i-1,2*j )+tmpvar(2*i,2*j ))
666  dy(i,j) = dx(i,j)
667  enddo
668  enddo
669 ! allocate(tmpvar(nx,ny+1))
670 
671 ! error=nf_inq_varid(ncid, 'dx', id_var)
672 ! call netcdf_err(error, 'inquire varid of dx from file '
673 ! & //trim(OUTGRID) )
674 ! error=nf_get_var_double(ncid, id_var, tmpvar)
675 ! call netcdf_err(error, 'inquire data of dx from file '
676 ! & //trim(OUTGRID) )
677 ! dx(1:IM,1:JM+1) = tmpvar(2:nx:2,1:ny+1:2)
678 ! deallocate(tmpvar)
679 
680 ! allocate(tmpvar(nx+1,ny))
681 ! error=nf_inq_varid(ncid, 'dy', id_var)
682 ! call netcdf_err(error, 'inquire varid of dy from file '
683 ! & //trim(OUTGRID) )
684 ! error=nf_get_var_double(ncid, id_var, tmpvar)
685 ! call netcdf_err(error, 'inquire data of dy from file '
686 ! & //trim(OUTGRID) )
687 ! dy(1:IM+1,1:JM) = tmpvar(1:nx+1:2,2:ny:2)
688  deallocate(tmpvar)
689  endif
690  tend=timef()
691  write(6,*)' Timer 1 time= ',tend-tbeg
692  !
693  if(grid_from_file) then
694  tbeg=timef()
695 
696  IF (merge_file == 'none') then
697  CALL make_mask(zslm,slm,land_frac,glat,
698  & im,jm,imn,jmn,geolon_c,geolat_c)
699  lake_frac=9999.9
700  ELSE
701  print*,'Read in external mask ',merge_file
702  CALL read_mask(merge_file,slm,land_frac,lake_frac,im,jm)
703  ENDIF
704 
705  IF (mask_only) THEN
706  print*,'Computing mask only.'
707  CALL write_mask_netcdf(im,jm,slm,land_frac,
708  1 1,1,geolon,geolat)
709 
710  print*,' DONE.'
711  stop
712  END IF
713 
714  CALL makemt2(zavg,zslm,oro,slm,var,var4,glat,
715  & im,jm,imn,jmn,geolon_c,geolat_c,lake_frac,land_frac)
716 
717  tend=timef()
718  write(6,*)' MAKEMT2 time= ',tend-tbeg
719  else
720  CALL makemt(zavg,zslm,oro,slm,var,var4,glat,
721  & ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
722  endif
723 
724  call minmxj(im,jm,oro,' ORO')
725  call minmxj(im,jm,slm,' SLM')
726  call minmxj(im,jm,var,' VAR')
727  call minmxj(im,jm,var4,' VAR4')
728 !
729 C check antarctic pole
730 ! DO J = 1,JM
731 ! DO I = 1,IM
732 ! if ( i .le. 100 .and. i .ge. 1 )then
733 ! if (j .ge. JM-1 )then
734 ! if (height .eq. 0.) print *,'I,J,SLM:',I,J,SLM(I,J)
735 ! write(6,153)i,j,ORO(i,j),HEIGHT,SLM(i,j)
736 ! endif
737 ! endif
738 ! ENDDO
739 ! ENDDO
740 C
741 C === Compute mtn principal coord HTENSR: THETA,GAMMA,SIGMA
742 C
743  allocate (theta(im,jm),gamma(im,jm),sigma(im,jm),elvmax(im,jm))
744  if(grid_from_file) then
745  tbeg=timef()
746  CALL makepc2(zavg,zslm,theta,gamma,sigma,glat,
747  1 im,jm,imn,jmn,geolon_c,geolat_c,slm)
748  tend=timef()
749  write(6,*)' MAKEPC2 time= ',tend-tbeg
750  else
751  CALL makepc(zavg,zslm,theta,gamma,sigma,glat,
752  1 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
753  endif
754 
755  call minmxj(im,jm,theta,' THETA')
756  call minmxj(im,jm,gamma,' GAMMA')
757  call minmxj(im,jm,sigma,' SIGMA')
758 
759 C
760 C COMPUTE MOUNTAIN DATA : OA OL
761 C
762  allocate (iwork(im,jm,4))
763  allocate (oa(im,jm,4),ol(im,jm,4),hprime(im,jm,14))
764  allocate (work1(im,jm),work2(im,jm),work3(im,jm),work4(im,jm))
765  allocate (work5(im,jm),work6(im,jm))
766 
767  call minmxj(im,jm,oro,' ORO')
768  print*, "inputorog=", trim(inputorog)
769  if(grid_from_file) then
770  if(trim(inputorog) == "none") then
771  print*, "calling MAKEOA2 to compute OA, OL"
772  tbeg=timef()
773  CALL makeoa2(zavg,zslm,var,glat,oa,ol,iwork,elvmax,oro,
774  1 work1,work2,work3,work4,work5,work6,
775  2 im,jm,imn,jmn,geolon_c,geolat_c,
776  3 geolon,geolat,dx,dy,is_south_pole,is_north_pole)
777  tend=timef()
778  write(6,*)' MAKEOA2 time= ',tend-tbeg
779  else
780  !-- read the data from INPUTOROG file.
781  error=nf__open(trim(inputorog),nf_nowrite,fsize,ncid)
782  call netcdf_err(error, 'Open file '//trim(inputorog) )
783  error=nf_inq_dimid(ncid, 'lon', id_dim)
784  call netcdf_err(error, 'inquire dimension lon from file '//
785  & trim(inputorog) )
786  error=nf_inq_dimlen(ncid,id_dim,nx_in)
787  call netcdf_err(error, 'inquire dimension lon length '//
788  & 'from file '//trim(inputorog) )
789  error=nf_inq_dimid(ncid, 'lat', id_dim)
790  call netcdf_err(error, 'inquire dimension lat from file '//
791  & trim(inputorog) )
792  error=nf_inq_dimlen(ncid,id_dim,ny_in)
793  call netcdf_err(error, 'inquire dimension lat length '//
794  & 'from file '//trim(inputorog) )
795 
796  print*, "extrapolate OA, OL from Gaussian grid with nx=",
797  & nx_in, ", ny=", ny_in
798  allocate(oa_in(nx_in,ny_in,4), ol_in(nx_in,ny_in,4))
799  allocate(slm_in(nx_in,ny_in) )
800  allocate(lon_in(nx_in,ny_in), lat_in(nx_in,ny_in) )
801 
802  error=nf_inq_varid(ncid, 'oa1', id_var)
803  call netcdf_err(error, 'inquire varid of oa1 from file '
804  & //trim(inputorog) )
805  error=nf_get_var_double(ncid, id_var, oa_in(:,:,1))
806  call netcdf_err(error, 'inquire data of oa1 from file '
807  & //trim(inputorog) )
808  error=nf_inq_varid(ncid, 'oa2', id_var)
809  call netcdf_err(error, 'inquire varid of oa2 from file '
810  & //trim(inputorog) )
811  error=nf_get_var_double(ncid, id_var, oa_in(:,:,2))
812  call netcdf_err(error, 'inquire data of oa2 from file '
813  & //trim(inputorog) )
814  error=nf_inq_varid(ncid, 'oa3', id_var)
815  call netcdf_err(error, 'inquire varid of oa3 from file '
816  & //trim(inputorog) )
817  error=nf_get_var_double(ncid, id_var, oa_in(:,:,3))
818  call netcdf_err(error, 'inquire data of oa3 from file '
819  & //trim(inputorog) )
820  error=nf_inq_varid(ncid, 'oa4', id_var)
821  call netcdf_err(error, 'inquire varid of oa4 from file '
822  & //trim(inputorog) )
823  error=nf_get_var_double(ncid, id_var, oa_in(:,:,4))
824  call netcdf_err(error, 'inquire data of oa4 from file '
825  & //trim(inputorog) )
826 
827  error=nf_inq_varid(ncid, 'ol1', id_var)
828  call netcdf_err(error, 'inquire varid of ol1 from file '
829  & //trim(inputorog) )
830  error=nf_get_var_double(ncid, id_var, ol_in(:,:,1))
831  call netcdf_err(error, 'inquire data of ol1 from file '
832  & //trim(inputorog) )
833  error=nf_inq_varid(ncid, 'ol2', id_var)
834  call netcdf_err(error, 'inquire varid of ol2 from file '
835  & //trim(inputorog) )
836  error=nf_get_var_double(ncid, id_var, ol_in(:,:,2))
837  call netcdf_err(error, 'inquire data of ol2 from file '
838  & //trim(inputorog) )
839  error=nf_inq_varid(ncid, 'ol3', id_var)
840  call netcdf_err(error, 'inquire varid of ol3 from file '
841  & //trim(inputorog) )
842  error=nf_get_var_double(ncid, id_var, ol_in(:,:,3))
843  call netcdf_err(error, 'inquire data of ol3 from file '
844  & //trim(inputorog) )
845  error=nf_inq_varid(ncid, 'ol4', id_var)
846  call netcdf_err(error, 'inquire varid of ol4 from file '
847  & //trim(inputorog) )
848  error=nf_get_var_double(ncid, id_var, ol_in(:,:,4))
849  call netcdf_err(error, 'inquire data of ol4 from file '
850  & //trim(inputorog) )
851 
852  error=nf_inq_varid(ncid, 'slmsk', id_var)
853  call netcdf_err(error, 'inquire varid of slmsk from file '
854  & //trim(inputorog) )
855  error=nf_get_var_double(ncid, id_var, slm_in)
856  call netcdf_err(error, 'inquire data of slmsk from file '
857  & //trim(inputorog) )
858 
859  error=nf_inq_varid(ncid, 'geolon', id_var)
860  call netcdf_err(error, 'inquire varid of geolon from file '
861  & //trim(inputorog) )
862  error=nf_get_var_double(ncid, id_var, lon_in)
863  call netcdf_err(error, 'inquire data of geolon from file '
864  & //trim(inputorog) )
865 
866  error=nf_inq_varid(ncid, 'geolat', id_var)
867  call netcdf_err(error, 'inquire varid of geolat from file '
868  & //trim(inputorog) )
869  error=nf_get_var_double(ncid, id_var, lat_in)
870  call netcdf_err(error, 'inquire data of geolat from file '
871  & //trim(inputorog) )
872 
873  ! set slmsk=2 to be ocean (0)
874  do j=1,ny_in; do i=1,nx_in
875  if(slm_in(i,j) == 2) slm_in(i,j) = 0
876  enddo; enddo
877 
878  error=nf_close(ncid)
879  call netcdf_err(error, 'close file '
880  & //trim(inputorog) )
881 
882  print*, "calling MAKEOA3 to compute OA, OL"
883  CALL makeoa3(zavg,var,glat,oa,ol,iwork,elvmax,oro,slm,
884  1 work1,work2,work3,work4,work5,work6,
885  2 im,jm,imn,jmn,geolon_c,geolat_c,
886  3 geolon,geolat,nx_in,ny_in,
887  4 oa_in,ol_in,slm_in,lon_in,lat_in)
888 
889  deallocate(oa_in,ol_in,slm_in,lon_in,lat_in)
890 
891  endif
892  else
893  CALL makeoa(zavg,var,glat,oa,ol,iwork,elvmax,oro,
894  1 work1,work2,work3,work4,
895  2 work5,work6,
896  3 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
897  endif
898 
899 ! Deallocate 2d vars
900  deallocate(ist,ien)
901  deallocate (zslm,zavg)
902  deallocate (dx,dy)
903  deallocate (work2,work3,work4,work5,work6)
904 
905 ! Deallocate 3d vars
906  deallocate(iwork)
907 
908  tbeg=timef()
909  call minmxj(im,jm,oa,' OA')
910  call minmxj(im,jm,ol,' OL')
911  call minmxj(im,jm,elvmax,' ELVMAX')
912  call minmxj(im,jm,oro,' ORO')
913 
914  maxc3 = 0
915  maxc4 = 0
916  maxc5 = 0
917  maxc6 = 0
918  maxc7 = 0
919  maxc8 = 0
920  DO j = 1,jm
921  DO i = 1,im
922  if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
923  if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
924  if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
925  if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
926  if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
927  if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
928  ENDDO
929  ENDDO
930  print *,' MAXC3:',maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
931 !
932 c itest=151
933 c jtest=56
934 C
935  print *,' ===> Replacing ELVMAX with ELVMAX-ORO <=== '
936  print *,' ===> if ELVMAX<=ORO replace with proxy <=== '
937  print *,' ===> the sum of mean orog (ORO) and std dev <=== '
938  DO j = 1,jm
939  DO i = 1,im
940  if (elvmax(i,j) .lt. oro(i,j) ) then
941 C--- subtracting off ORO leaves std dev (this should never happen)
942  elvmax(i,j) = max( 3. * var(i,j),0.)
943  else
944  elvmax(i,j) = max( elvmax(i,j) - oro(i,j),0.)
945  endif
946  ENDDO
947  ENDDO
948  maxc3 = 0
949  maxc4 = 0
950  maxc5 = 0
951  maxc6 = 0
952  maxc7 = 0
953  maxc8 = 0
954  DO j = 1,jm
955  DO i = 1,im
956  if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
957  if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
958  if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
959  if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
960  if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
961  if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
962  ENDDO
963  ENDDO
964  print *,' after MAXC 3-6 km:',maxc3,maxc4,maxc5,maxc6
965 c
966  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
967 ! if (JM .gt. 0) stop
968 C
969 C ZERO OVER OCEAN
970 C
971  print *,' Testing at point (itest,jtest)=',itest,jtest
972  print *,' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
973  print *,' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
974  DO j = 1,jm
975  DO i = 1,im
976  IF(slm(i,j).EQ.0.) THEN
977 C VAR(I,J) = 0.
978  var4(i,j) = 0.
979  oa(i,j,1) = 0.
980  oa(i,j,2) = 0.
981  oa(i,j,3) = 0.
982  oa(i,j,4) = 0.
983  ol(i,j,1) = 0.
984  ol(i,j,2) = 0.
985  ol(i,j,3) = 0.
986  ol(i,j,4) = 0.
987 C THETA(I,J) =0.
988 C GAMMA(I,J) =0.
989 C SIGMA(I,J) =0.
990 C ELVMAX(I,J)=0.
991 ! --- the sub-grid scale parameters for mtn blocking and gwd retain
992 ! --- properties even if over ocean but there is elevation within the
993 ! --- gaussian grid box.
994  ENDIF
995  ENDDO
996  ENDDO
997 C
998 ! --- if mskocn=1 ocean land sea mask given, =0 if not present
999 ! --- OCLSM is real(*4) array with fractional values possible
1000 ! --- 0 is ocean and 1 is land for slm
1001 ! --- Step 1: Only change SLM after GFS SLM is applied
1002 ! --- SLM is only field that will be altered by OCLSM
1003 ! --- Ocean land sea mask ocean points made ocean in atm model
1004 ! --- Land and Lakes and all other atm elv moments remain unchanged.
1005 
1006  IF (merge_file == 'none') then
1007 
1008  msk_ocn : if ( mskocn .eq. 1 ) then
1009 
1010  DO j = 1,jm
1011  DO i = 1,im
1012  if (abs(oro(i,j)) .lt. 1. ) then
1013  slm(i,j) = slmi(i,j)
1014  else
1015  if ( slmi(i,j) .eq. 1. .and. slm(i,j) .eq. 1) slm(i,j) = 1
1016  if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
1017  if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 1) slm(i,j) = 0
1018  if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
1019  endif
1020  enddo
1021  enddo
1022  endif msk_ocn
1023  endif
1024  print *,' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
1025  print *,' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
1026 
1027  deallocate(slmi)
1028 
1029  IF (merge_file == 'none') then
1030 
1031 C REMOVE ISOLATED POINTS
1032  iso_loop : DO j=2,jm-1
1033  jn=j-1
1034  js=j+1
1035  DO i=1,im
1036  iw=mod(i+im-2,im)+1
1037  ie=mod(i,im)+1
1038  slma=slm(iw,j)+slm(ie,j)
1039  oroa=oro(iw,j)+oro(ie,j)
1040  vara=var(iw,j)+var(ie,j)
1041  var4a=var4(iw,j)+var4(ie,j)
1042  DO k=1,4
1043  oaa(k)=oa(iw,j,k)+oa(ie,j,k)
1044 ! --- (*j*) fix typo:
1045  ola(k)=ol(iw,j,k)+ol(ie,j,k)
1046  ENDDO
1047  wgta=2
1048  xn=(i-1)+1
1049  IF(abs(xn-nint(xn)).LT.1.e-2) THEN
1050  in=mod(nint(xn)-1,im)+1
1051  inw=mod(in+im-2,im)+1
1052  ine=mod(in,im)+1
1053  slma=slma+slm(inw,jn)+slm(in,jn)+slm(ine,jn)
1054  oroa=oroa+oro(inw,jn)+oro(in,jn)+oro(ine,jn)
1055  vara=vara+var(inw,jn)+var(in,jn)+var(ine,jn)
1056  var4a=var4a+var4(inw,jn)+var4(in,jn)+var4(ine,jn)
1057  DO k=1,4
1058  oaa(k)=oaa(k)+oa(inw,jn,k)+oa(in,jn,k)+oa(ine,jn,k)
1059  ola(k)=ola(k)+ol(inw,jn,k)+ol(in,jn,k)+ol(ine,jn,k)
1060  ENDDO
1061  wgta=wgta+3
1062  ELSE
1063  inw=int(xn)
1064  ine=mod(inw,im)+1
1065  slma=slma+slm(inw,jn)+slm(ine,jn)
1066  oroa=oroa+oro(inw,jn)+oro(ine,jn)
1067  vara=vara+var(inw,jn)+var(ine,jn)
1068  var4a=var4a+var4(inw,jn)+var4(ine,jn)
1069  DO k=1,4
1070  oaa(k)=oaa(k)+oa(inw,jn,k)+oa(ine,jn,k)
1071  ola(k)=ola(k)+ol(inw,jn,k)+ol(ine,jn,k)
1072  ENDDO
1073  wgta=wgta+2
1074  ENDIF
1075  xs=(i-1)+1
1076  IF(abs(xs-nint(xs)).LT.1.e-2) THEN
1077  is=mod(nint(xs)-1,im)+1
1078  isw=mod(is+im-2,im)+1
1079  ise=mod(is,im)+1
1080  slma=slma+slm(isw,js)+slm(is,js)+slm(ise,js)
1081  oroa=oroa+oro(isw,js)+oro(is,js)+oro(ise,js)
1082  vara=vara+var(isw,js)+var(is,js)+var(ise,js)
1083  var4a=var4a+var4(isw,js)+var4(is,js)+var4(ise,js)
1084  DO k=1,4
1085  oaa(k)=oaa(k)+oa(isw,js,k)+oa(is,js,k)+oa(ise,js,k)
1086  ola(k)=ola(k)+ol(isw,js,k)+ol(is,js,k)+ol(ise,js,k)
1087  ENDDO
1088  wgta=wgta+3
1089  ELSE
1090  isw=int(xs)
1091  ise=mod(isw,im)+1
1092  slma=slma+slm(isw,js)+slm(ise,js)
1093  oroa=oroa+oro(isw,js)+oro(ise,js)
1094  vara=vara+var(isw,js)+var(ise,js)
1095  var4a=var4a+var4(isw,js)+var4(ise,js)
1096  DO k=1,4
1097  oaa(k)=oaa(k)+oa(isw,js,k)+oa(ise,js,k)
1098  ola(k)=ola(k)+ol(isw,js,k)+ol(ise,js,k)
1099  ENDDO
1100  wgta=wgta+2
1101  ENDIF
1102  oroa=oroa/wgta
1103  vara=vara/wgta
1104  var4a=var4a/wgta
1105  DO k=1,4
1106  oaa(k)=oaa(k)/wgta
1107  ola(k)=ola(k)/wgta
1108  ENDDO
1109  IF(slm(i,j).EQ.0..AND.slma.EQ.wgta) THEN
1110  print '("SEA ",2F8.0," MODIFIED TO LAND",2F8.0," AT ",2I8)',
1111  & oro(i,j),var(i,j),oroa,vara,i,j
1112  slm(i,j)=1.
1113  oro(i,j)=oroa
1114  var(i,j)=vara
1115  var4(i,j)=var4a
1116  DO k=1,4
1117  oa(i,j,k)=oaa(k)
1118  ol(i,j,k)=ola(k)
1119  ENDDO
1120  ELSEIF(slm(i,j).EQ.1..AND.slma.EQ.0.) THEN
1121  print '("LAND",2F8.0," MODIFIED TO SEA ",2F8.0," AT ",2I8)',
1122  & oro(i,j),var(i,j),oroa,vara,i,j
1123  slm(i,j)=0.
1124  oro(i,j)=oroa
1125  var(i,j)=vara
1126  var4(i,j)=var4a
1127  DO k=1,4
1128  oa(i,j,k)=oaa(k)
1129  ol(i,j,k)=ola(k)
1130  ENDDO
1131  ENDIF
1132  ENDDO
1133  ENDDO iso_loop
1134 C--- print for testing after isolated points removed
1135  print *,' after isolated points removed'
1136  call minmxj(im,jm,oro,' ORO')
1137  print *,' ORO(itest,jtest)=',oro(itest,jtest)
1138  print *,' VAR(itest,jtest)=',var(itest,jtest)
1139  print *,' VAR4(itest,jtest)=',var4(itest,jtest)
1140  print *,' OA(itest,jtest,1)=',oa(itest,jtest,1)
1141  print *,' OA(itest,jtest,2)=',oa(itest,jtest,2)
1142  print *,' OA(itest,jtest,3)=',oa(itest,jtest,3)
1143  print *,' OA(itest,jtest,4)=',oa(itest,jtest,4)
1144  print *,' OL(itest,jtest,1)=',ol(itest,jtest,1)
1145  print *,' OL(itest,jtest,2)=',ol(itest,jtest,2)
1146  print *,' OL(itest,jtest,3)=',ol(itest,jtest,3)
1147  print *,' OL(itest,jtest,4)=',ol(itest,jtest,4)
1148  print *,' Testing at point (itest,jtest)=',itest,jtest
1149  print *,' THETA(itest,jtest)=',theta(itest,jtest)
1150  print *,' GAMMA(itest,jtest)=',gamma(itest,jtest)
1151  print *,' SIGMA(itest,jtest)=',sigma(itest,jtest)
1152  print *,' ELVMAX(itest,jtest)=',elvmax(itest,jtest)
1153  print *,' EFAC=',efac
1154 
1155  endif
1156 
1157 C
1158  DO j=1,jm
1159  DO i=1,im
1160  oro(i,j) = oro(i,j) + efac*var(i,j)
1161  hprime(i,j,1) = var(i,j)
1162  hprime(i,j,2) = var4(i,j)
1163  hprime(i,j,3) = oa(i,j,1)
1164  hprime(i,j,4) = oa(i,j,2)
1165  hprime(i,j,5) = oa(i,j,3)
1166  hprime(i,j,6) = oa(i,j,4)
1167  hprime(i,j,7) = ol(i,j,1)
1168  hprime(i,j,8) = ol(i,j,2)
1169  hprime(i,j,9) = ol(i,j,3)
1170  hprime(i,j,10)= ol(i,j,4)
1171  hprime(i,j,11)= theta(i,j)
1172  hprime(i,j,12)= gamma(i,j)
1173  hprime(i,j,13)= sigma(i,j)
1174  hprime(i,j,14)= elvmax(i,j)
1175  ENDDO
1176  ENDDO
1177 !
1178  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1179 ! --- Quadratic filter applied by default.
1180 ! --- NF0 is normally set to an even value beyond the previous truncation,
1181 ! --- for example, for jcap=382, NF0=254+2
1182 ! --- NF1 is set as jcap+2 (and/or nearest even), eg., for t382, NF1=382+2=384
1183 ! --- if no filter is desired then NF1=NF0=0 and ORF=ORO
1184 ! --- if no filter but spectral to grid (with gibbs) then NF1=jcap+2, and NF1=jcap+1
1185 !
1186  deallocate(var4)
1187  allocate (orf(im,jm))
1188  IF ( nf1 - nf0 .eq. 0 ) filter=.false.
1189  print *,' NF1, NF0, FILTER=',nf1,nf0,filter
1190  IF (filter) THEN
1191 C SPECTRALLY TRUNCATE AND FILTER OROGRAPHY
1192 
1193  CALL sptez(nr,nm,4,im,jm,ors,oro,-1)
1194 !
1195  print *,' about to apply spectral filter '
1196  fff=1./(nf1-nf0)**2
1197  i=0
1198  DO m=0,nm
1199  DO n=m,nm+nr*m
1200  IF(n.GT.nf0) THEN
1201  www=max(1.-fff*(n-nf0)**2,0.)
1202  ors(i+1)=ors(i+1)*www
1203  ors(i+2)=ors(i+2)*www
1204  ENDIF
1205  i=i+2
1206  ENDDO
1207  ENDDO
1208 !
1209  CALL sptez(nr,nm,4,im,jm,ors,orf,+1)
1210 
1211  ELSE
1212  ors=0.
1213  orf=oro
1214  ENDIF
1215 
1216  deallocate (work1)
1217 
1218  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1219  print *,' ELVMAX(',itest,jtest,')=',elvmax(itest,jtest)
1220  print *,' after spectral filter is applied'
1221  call minmxj(im,jm,oro,' ORO')
1222  call minmxj(im,jm,orf,' ORF')
1223 C
1224  print *,' after nearest neighbor interpolation applied '
1225  call minmxj(im,jm,oro,' ORO')
1226  call minmxj(im,jm,orf,' ORF')
1227  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1228  print *,' ORO,ORF(itest,jtest),itest,jtest:',
1229  & oro(itest,jtest),orf(itest,jtest),itest,jtest
1230  print *,' ELVMAX(',itest,jtest,')=',elvmax(itest,jtest)
1231 
1232 
1233 C check antarctic pole
1234  DO j = 1,jm
1235  DO i = 1,im
1236  if ( i .le. 21 .and. i .ge. 1 )then
1237  if (j .eq. jm )write(6,153)i,j,oro(i,j),elvmax(i,j),slm(i,j)
1238  153 format(1x,' ORO,ELVMAX(i=',i4,' j=',i4,')=',2e14.5,f5.1)
1239  endif
1240  ENDDO
1241  ENDDO
1242  tend=timef()
1243  write(6,*)' Timer 5 time= ',tend-tbeg
1244 C
1245  delxn = 360./im
1246  do i=1,im
1247  xlon(i) = delxn*(i-1)
1248  enddo
1249  IF(trim(outgrid) == "none") THEN
1250  do j=1,jm
1251  do i=1,im
1252  geolon(i,j) = xlon(i)
1253  geolat(i,j) = xlat(j)
1254  enddo
1255  enddo
1256  else
1257  do j = 1, jm
1258  xlat(j) = geolat(1,j)
1259  enddo
1260  do i = 1, im
1261  xlon(i) = geolon(i,1)
1262  enddo
1263  endif
1264 
1265  tbeg=timef()
1266  CALL write_netcdf(im,jm,slm,land_frac,oro,orf,hprime,1,1,
1267  1 geolon(1:im,1:jm),geolat(1:im,1:jm), xlon,xlat)
1268  tend=timef()
1269  write(6,*)' WRITE_NETCDF time= ',tend-tbeg
1270  print *,' wrote netcdf file out.oro.tile?.nc'
1271 
1272  print *,' ===== Deallocate Arrays and ENDING MTN VAR OROG program'
1273 
1274 ! Deallocate 1d vars
1275  deallocate(jst,jen,numi)
1276  deallocate(cosclt,wgtclt,rclt,xlat,diffx,xlon,ors,oaa,ola,glat)
1277 
1278 ! Deallocate 2d vars
1279  deallocate (oclsm)
1280  deallocate (geolon,geolon_c,geolat,geolat_c)
1281  deallocate (slm,oro,var,orf,land_frac)
1282  deallocate (theta,gamma,sigma,elvmax)
1283 
1284 
1285  tend=timef()
1286  write(6,*)' Total runtime time= ',tend-tbeg1
1287  RETURN
1288  END SUBROUTINE tersub
1289 
1319  SUBROUTINE makemt(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
1320  1 GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
1321  dimension glat(jmn),xlat(jm)
1322 ! REAL*4 OCLSM
1323  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
1324  DIMENSION ORO(IM,JM),SLM(IM,JM),VAR(IM,JM),VAR4(IM,JM)
1325  dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
1326  LOGICAL FLAG
1327 C
1328 ! ---- OCLSM holds the ocean (im,jm) grid
1329  print *,' _____ SUBROUTINE MAKEMT '
1330 C---- GLOBAL XLAT AND XLON ( DEGREE )
1331 C
1332  jm1 = jm - 1
1333  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1334 C
1335  DO j=1,jmn
1336  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1337  ENDDO
1338 C
1339 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1340 C
1341 C (*j*) for hard wired zero offset (lambda s =0) for terr05
1342  DO j=1,jm
1343  DO i=1,numi(j)
1344  im1 = numi(j) - 1
1345  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
1346  faclon = delx / delxn
1347  ist(i,j) = faclon * float(i-1) - faclon * 0.5 + 1
1348  ien(i,j) = faclon * float(i) - faclon * 0.5 + 1
1349 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
1350 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
1351 C
1352  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
1353  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
1354 !
1355 ! if ( I .lt. 10 .and. J .ge. JM-1 )
1356 ! 1 PRINT*,' MAKEMT: I j IST IEN ',I,j,IST(I,j),IEN(I,j)
1357  ENDDO
1358 ! if ( J .ge. JM-1 ) then
1359 ! print *,' *** FACLON=',FACLON, 'numi(j=',j,')=',numi(j)
1360 ! endif
1361  ENDDO
1362  print *,' DELX=',delx,' DELXN=',delxn
1363  DO j=1,jm-1
1364  flag=.true.
1365  DO j1=1,jmn
1366  xxlat = (xlat(j)+xlat(j+1))/2.
1367  IF(flag.AND.glat(j1).GT.xxlat) THEN
1368  jst(j) = j1
1369  jen(j+1) = j1 - 1
1370  flag = .false.
1371  ENDIF
1372  ENDDO
1373 CX PRINT*, ' J JST JEN ',J,JST(J),JEN(J),XLAT(J),GLAT(J1)
1374  ENDDO
1375  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
1376  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
1377 ! PRINT*, ' JM JST JEN=',JST(JM),JEN(JM),XLAT(JM),GLAT(JMN)
1378 C
1379 C...FIRST, AVERAGED HEIGHT
1380 C
1381  DO j=1,jm
1382  DO i=1,numi(j)
1383  oro(i,j) = 0.0
1384  var(i,j) = 0.0
1385  var4(i,j) = 0.0
1386  xnsum = 0.0
1387  xland = 0.0
1388  xwatr = 0.0
1389  xl1 = 0.0
1390  xs1 = 0.0
1391  xw1 = 0.0
1392  xw2 = 0.0
1393  xw4 = 0.0
1394  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
1395  i1 = ist(i,j) + ii1 - 1
1396  IF(i1.LE.0.) i1 = i1 + imn
1397  IF(i1.GT.imn) i1 = i1 - imn
1398 ! if ( i .le. 10 .and. i .ge. 1 ) then
1399 ! if (j .eq. JM )
1400 ! &print *,' J,JST,JEN,IST,IEN,I1=',
1401 ! &J,JST(j),JEN(J),IST(I,j),IEN(I,j),I1
1402 ! endif
1403  DO j1=jst(j),jen(j)
1404  xland = xland + float(zslm(i1,j1))
1405  xwatr = xwatr + float(1-zslm(i1,j1))
1406  xnsum = xnsum + 1.
1407  height = float(zavg(i1,j1))
1408 C.........
1409  IF(height.LT.-990.) height = 0.0
1410  xl1 = xl1 + height * float(zslm(i1,j1))
1411  xs1 = xs1 + height * float(1-zslm(i1,j1))
1412  xw1 = xw1 + height
1413  xw2 = xw2 + height ** 2
1414 C check antarctic pole
1415 ! if ( i .le. 10 .and. i .ge. 1 )then
1416 ! if (j .ge. JM-1 )then
1417 C=== degub testing
1418 ! print *," I,J,I1,J1,XL1,XS1,XW1,XW2:",I,J,I1,J1,XL1,XS1,XW1,XW2
1419 ! 153 format(1x,' ORO,ELVMAX(i=',i4,' j=',i4,')=',2E14.5,3f5.1)
1420 ! endif
1421 ! endif
1422  ENDDO
1423  ENDDO
1424  IF(xnsum.GT.1.) THEN
1425 ! --- SLM initialized with OCLSM calc from all land points except ....
1426 ! --- 0 is ocean and 1 is land for slm
1427 ! --- Step 1 is to only change SLM after GFS SLM is applied
1428 
1429  slm(i,j) = float(nint(xland/xnsum))
1430  IF(slm(i,j).NE.0.) THEN
1431  oro(i,j)= xl1 / xland
1432  ELSE
1433  oro(i,j)= xs1 / xwatr
1434  ENDIF
1435  var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1436  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
1437  i1 = ist(i,j) + ii1 - 1
1438  IF(i1.LE.0.) i1 = i1 + imn
1439  IF(i1.GT.imn) i1 = i1 - imn
1440  DO j1=jst(j),jen(j)
1441  height = float(zavg(i1,j1))
1442  IF(height.LT.-990.) height = 0.0
1443  xw4 = xw4 + (height-oro(i,j)) ** 4
1444  ENDDO
1445  ENDDO
1446  IF(var(i,j).GT.1.) THEN
1447 ! if ( I .lt. 20 .and. J .ge. JM-19 ) then
1448 ! print *,'I,J,XW4,XNSUM,VAR(I,J)',I,J,XW4,XNSUM,VAR(I,J)
1449 ! endif
1450  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
1451  ENDIF
1452  ENDIF
1453  ENDDO
1454  ENDDO
1455  WRITE(6,*) "! MAKEMT ORO SLM VAR VAR4 DONE"
1456 C
1457 
1458  RETURN
1459  END
1460 
1480  SUBROUTINE get_index(IMN,JMN,npts,lonO,latO,DELXN,
1481  & jst,jen,ilist,numx)
1482  implicit none
1483  integer, intent(in) :: IMN,JMN
1484  integer :: npts
1485  real, intent(in) :: LONO(npts), LATO(npts)
1486  real, intent(in) :: DELXN
1487  integer, intent(out) :: jst,jen
1488  integer, intent(out) :: ilist(IMN)
1489  integer, intent(out) :: numx
1490  real minlat,maxlat,minlon,maxlon
1491  integer i2, ii, ist, ien
1492 
1493  minlat = minval(lato)
1494  maxlat = maxval(lato)
1495  minlon = minval(lono)
1496  maxlon = maxval(lono)
1497  ist = minlon/delxn+1
1498  ien = maxlon/delxn+1
1499  jst = (minlat+90)/delxn+1
1500  jen = (maxlat+90)/delxn
1501  !--- add a few points to both ends of j-direction
1502  jst = jst - 5
1503  if(jst<1) jst = 1
1504  jen = jen + 5
1505  if(jen>jmn) jen = jmn
1506 
1507  !--- when around the pole, just search through all the points.
1508  if((jst == 1 .OR. jen == jmn) .and.
1509  & (ien-ist+1 > imn/2) )then
1510  numx = imn
1511  do i2 = 1, imn
1512  ilist(i2) = i2
1513  enddo
1514  else if( ien-ist+1 > imn/2 ) then ! cross longitude = 0
1515  !--- find the minimum that greater than IMN/2
1516  !--- and maximum that less than IMN/2
1517  ist = 0
1518  ien = imn+1
1519  do i2 = 1, npts
1520  ii = lono(i2)/delxn+1
1521  if(ii <0 .or. ii>imn) print*,"ii=",ii,imn,lono(i2),delxn
1522  if( ii < imn/2 ) then
1523  ist = max(ist,ii)
1524  else if( ii > imn/2 ) then
1525  ien = min(ien,ii)
1526  endif
1527  enddo
1528  if(ist<1 .OR. ist>imn) then
1529  print*, .or."FATAL ERROR: ist<1 ist>IMN"
1530  call abort()
1531  endif
1532  if(ien<1 .OR. ien>imn) then
1533  print*, .or."FATAL ERROR: iend<1 iend>IMN"
1534  call abort()
1535  endif
1536 
1537  numx = imn - ien + 1
1538  do i2 = 1, numx
1539  ilist(i2) = ien + (i2-1)
1540  enddo
1541  do i2 = 1, ist
1542  ilist(numx+i2) = i2
1543  enddo
1544  numx = numx+ist
1545  else
1546  numx = ien-ist+1
1547  do i2 = 1, numx
1548  ilist(i2) = ist + (i2-1)
1549  enddo
1550  endif
1551 
1552  END
1553 
1569  SUBROUTINE make_mask(zslm,SLM,land_frac,
1570  1 GLAT,IM,JM,IMN,JMN,lon_c,lat_c)
1571  implicit none
1572  real, parameter :: D2R = 3.14159265358979/180.
1573  integer, parameter :: MAXSUM=20000000
1574  integer im, jm, imn, jmn, jst, jen
1575  real GLAT(JMN), GLON(IMN)
1576  INTEGER ZSLM(IMN,JMN)
1577  real land_frac(IM,JM)
1578  real SLM(IM,JM)
1579  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
1580  real LONO(4),LATO(4),LONI,LATI
1581  integer JM1,i,j,nsum,nsum_all,ii,jj,numx,i2
1582  integer ilist(IMN)
1583  real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1
1584  real XNSUM_ALL,XLAND_ALL,XWATR_ALL
1585  logical inside_a_polygon
1586 C
1587 ! ---- OCLSM holds the ocean (im,jm) grid
1588  print *,' _____ SUBROUTINE MAKE_MASK '
1589 C---- GLOBAL XLAT AND XLON ( DEGREE )
1590 C
1591  jm1 = jm - 1
1592  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1593 C
1594  DO j=1,jmn
1595  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1596  ENDDO
1597  DO i=1,imn
1598  glon(i) = 0. + (i-1) * delxn + delxn * 0.5
1599  ENDDO
1600 
1601  land_frac(:,:) = 0.0
1602 C
1603 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1604 C
1605 C (*j*) for hard wired zero offset (lambda s =0) for terr05
1606 !$omp parallel do
1607 !$omp* private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,lono,
1608 !$omp* lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,
1609 !$omp* xnsum_all,xland_all,xwatr_all,nsum_all)
1610 !$omp*
1611  DO j=1,jm
1612 ! print*, "J=", J
1613  DO i=1,im
1614  xnsum = 0.0
1615  xland = 0.0
1616  xwatr = 0.0
1617  nsum = 0
1618  xnsum_all = 0.0
1619  xland_all = 0.0
1620  xwatr_all = 0.0
1621  nsum_all = 0
1622 
1623  lono(1) = lon_c(i,j)
1624  lono(2) = lon_c(i+1,j)
1625  lono(3) = lon_c(i+1,j+1)
1626  lono(4) = lon_c(i,j+1)
1627  lato(1) = lat_c(i,j)
1628  lato(2) = lat_c(i+1,j)
1629  lato(3) = lat_c(i+1,j+1)
1630  lato(4) = lat_c(i,j+1)
1631  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
1632  do jj = jst, jen; do i2 = 1, numx
1633  ii = ilist(i2)
1634  loni = ii*delxn
1635  lati = -90 + jj*delxn
1636 
1637  xland_all = xland_all + float(zslm(ii,jj))
1638  xwatr_all = xwatr_all + float(1-zslm(ii,jj))
1639  xnsum_all = xnsum_all + 1.
1640  nsum_all = nsum_all+1
1641  if(nsum_all > maxsum) then
1642  print*, "FATAL ERROR: nsum_all is greater than MAXSUM,"
1643  print*, "increase MAXSUM."
1644  call abort()
1645  endif
1646 
1647  if(inside_a_polygon(loni*d2r,lati*d2r,4,
1648  & lono*d2r,lato*d2r))then
1649 
1650  xland = xland + float(zslm(ii,jj))
1651  xwatr = xwatr + float(1-zslm(ii,jj))
1652  xnsum = xnsum + 1.
1653  nsum = nsum+1
1654  if(nsum > maxsum) then
1655  print*, "FATAL ERROR: nsum is greater than MAXSUM,"
1656  print*, "increase MAXSUM."
1657  call abort()
1658  endif
1659  endif
1660  enddo ; enddo
1661 
1662 
1663  IF(xnsum.GT.1.) THEN
1664 ! --- SLM initialized with OCLSM calc from all land points except ....
1665 ! --- 0 is ocean and 1 is land for slm
1666 ! --- Step 1 is to only change SLM after GFS SLM is applied
1667  land_frac(i,j) = xland/xnsum
1668  slm(i,j) = float(nint(xland/xnsum))
1669  ELSEIF(xnsum_all.GT.1.) THEN
1670  land_frac(i,j) = xland_all/xnsum _all
1671  slm(i,j) = float(nint(xland_all/xnsum_all))
1672  ELSE
1673  print*, "FATAL ERROR: no source points in MAKE_MASK."
1674  call abort()
1675  ENDIF
1676  ENDDO
1677  ENDDO
1678 !$omp end parallel do
1679  WRITE(6,*) "! MAKE_MASK DONE"
1680 C
1681  RETURN
1682  END SUBROUTINE make_mask
1704  SUBROUTINE makemt2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
1705  1 GLAT,IM,JM,IMN,JMN,lon_c,lat_c,lake_frac,land_frac)
1706  implicit none
1707  real, parameter :: D2R = 3.14159265358979/180.
1708  integer, parameter :: maxsum=20000000
1709  real, dimension(:), allocatable :: hgt_1d, hgt_1d_all
1710  integer IM, JM, IMN, JMN
1711  real GLAT(JMN), GLON(IMN)
1712  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
1713  real ORO(IM,JM),VAR(IM,JM),VAR4(IM,JM)
1714  real, intent(in) :: SLM(IM,JM), lake_frac(im,jm),land_frac(im,jm)
1715  integer JST, JEN
1716  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
1717  real LONO(4),LATO(4),LONI,LATI
1718  real HEIGHT
1719  integer JM1,i,j,nsum,nsum_all,ii,jj,i1,numx,i2
1720  integer ilist(IMN)
1721  real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1,XW2,XW4
1722  real XNSUM_ALL,XLAND_ALL,XWATR_ALL,HEIGHT_ALL
1723  real XL1_ALL,XS1_ALL,XW1_ALL,XW2_ALL,XW4_ALL
1724  logical inside_a_polygon
1725 C
1726 ! ---- OCLSM holds the ocean (im,jm) grid
1727 ! --- mskocn=1 Use ocean model sea land mask, OK and present,
1728 ! --- mskocn=0 dont use Ocean model sea land mask, not OK, not present
1729  print *,' _____ SUBROUTINE MAKEMT2 '
1730  allocate(hgt_1d(maxsum))
1731  allocate(hgt_1d_all(maxsum))
1732 C---- GLOBAL XLAT AND XLON ( DEGREE )
1733 C
1734  jm1 = jm - 1
1735  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1736 C
1737  DO j=1,jmn
1738  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1739  ENDDO
1740  DO i=1,imn
1741  glon(i) = 0. + (i-1) * delxn + delxn * 0.5
1742  ENDDO
1743 
1744 ! land_frac(:,:) = 0.0
1745 C
1746 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1747 C
1748 C (*j*) for hard wired zero offset (lambda s =0) for terr05
1749 !$omp parallel do
1750 !$omp* private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,xw2,xw4,lono,
1751 !$omp* lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,height,
1752 !$omp* hgt_1d,
1753 !$omp* xnsum_all,xland_all,xwatr_all,nsum_all,
1754 !$omp* xl1_all,xs1_all,xw1_all,xw2_all,xw4_all,
1755 !$omp* height_all,hgt_1d_all)
1756  DO j=1,jm
1757 ! print*, "J=", J
1758  DO i=1,im
1759  oro(i,j) = 0.0
1760  var(i,j) = 0.0
1761  var4(i,j) = 0.0
1762  xnsum = 0.0
1763  xland = 0.0
1764  xwatr = 0.0
1765  nsum = 0
1766  xl1 = 0.0
1767  xs1 = 0.0
1768  xw1 = 0.0
1769  xw2 = 0.0
1770  xw4 = 0.0
1771  xnsum_all = 0.0
1772  xland_all = 0.0
1773  xwatr_all = 0.0
1774  nsum_all = 0
1775  xl1_all = 0.0
1776  xs1_all = 0.0
1777  xw1_all = 0.0
1778  xw2_all = 0.0
1779  xw4_all = 0.0
1780 
1781  lono(1) = lon_c(i,j)
1782  lono(2) = lon_c(i+1,j)
1783  lono(3) = lon_c(i+1,j+1)
1784  lono(4) = lon_c(i,j+1)
1785  lato(1) = lat_c(i,j)
1786  lato(2) = lat_c(i+1,j)
1787  lato(3) = lat_c(i+1,j+1)
1788  lato(4) = lat_c(i,j+1)
1789  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
1790  do jj = jst, jen; do i2 = 1, numx
1791  ii = ilist(i2)
1792  loni = ii*delxn
1793  lati = -90 + jj*delxn
1794 
1795  xland_all = xland_all + float(zslm(ii,jj))
1796  xwatr_all = xwatr_all + float(1-zslm(ii,jj))
1797  xnsum_all = xnsum_all + 1.
1798  height_all = float(zavg(ii,jj))
1799  nsum_all = nsum_all+1
1800  if(nsum_all > maxsum) then
1801  print*, "FATAL ERROR: nsum_all is greater than MAXSUM,"
1802  print*, "increase MAXSUM."
1803  call abort()
1804  endif
1805  hgt_1d_all(nsum_all) = height_all
1806  IF(height_all.LT.-990.) height_all = 0.0
1807  xl1_all = xl1_all + height_all * float(zslm(ii,jj))
1808  xs1_all = xs1_all + height_all * float(1-zslm(ii,jj))
1809  xw1_all = xw1_all + height_all
1810  xw2_all = xw2_all + height_all ** 2
1811 
1812  if(inside_a_polygon(loni*d2r,lati*d2r,4,
1813  & lono*d2r,lato*d2r))then
1814 
1815  xland = xland + float(zslm(ii,jj))
1816  xwatr = xwatr + float(1-zslm(ii,jj))
1817  xnsum = xnsum + 1.
1818  height = float(zavg(ii,jj))
1819  nsum = nsum+1
1820  if(nsum > maxsum) then
1821  print*, "FATAL ERROR: nsum is greater than MAXSUM,"
1822  print*, "increase MAXSUM."
1823  call abort()
1824  endif
1825  hgt_1d(nsum) = height
1826  IF(height.LT.-990.) height = 0.0
1827  xl1 = xl1 + height * float(zslm(ii,jj))
1828  xs1 = xs1 + height * float(1-zslm(ii,jj))
1829  xw1 = xw1 + height
1830  xw2 = xw2 + height ** 2
1831  endif
1832  enddo ; enddo
1833 
1834  IF(xnsum.GT.1.) THEN
1835 ! --- SLM initialized with OCLSM calc from all land points except ....
1836 ! --- 0 is ocean and 1 is land for slm
1837 ! --- Step 1 is to only change SLM after GFS SLM is applied
1838 
1839  !IF(SLM(I,J).NE.0.) THEN
1840  IF(slm(i,j) .NE. 0. .OR. land_frac(i,j) > 0.) THEN
1841  IF (xland > 0) THEN
1842  oro(i,j)= xl1 / xland
1843  ELSE
1844  oro(i,j)= xs1 / xwatr
1845  ENDIF
1846  ELSE
1847  IF (xwatr > 0) THEN
1848  oro(i,j)= xs1 / xwatr
1849  ELSE
1850  oro(i,j)= xl1 / xland
1851  ENDIF
1852  ENDIF
1853 
1854  var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1855  do i1 = 1, nsum
1856  xw4 = xw4 + (hgt_1d(i1) - oro(i,j)) ** 4
1857  enddo
1858 
1859  IF(var(i,j).GT.1.) THEN
1860  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
1861  ENDIF
1862 
1863  ELSEIF(xnsum_all.GT.1.) THEN
1864 
1865  !IF(SLM(I,J).NE.0.) THEN
1866  IF(slm(i,j) .NE. 0. .OR. land_frac(i,j) > 0.) THEN
1867  IF (xland_all > 0) THEN
1868  oro(i,j)= xl1_all / xland_all
1869  ELSE
1870  oro(i,j)= xs1_all / xwatr_all
1871  ENDIF
1872  ELSE
1873  IF (xwatr_all > 0) THEN
1874  oro(i,j)= xs1_all / xwatr_all
1875  ELSE
1876  oro(i,j)= xl1_all / xland_all
1877  ENDIF
1878  ENDIF
1879 
1880  var(i,j)=sqrt(max(xw2_all/xnsum_all-
1881  & (xw1_all/xnsum_all)**2,0.))
1882  do i1 = 1, nsum_all
1883  xw4_all = xw4_all +
1884  & (hgt_1d_all(i1) - oro(i,j)) ** 4
1885  enddo
1886 
1887  IF(var(i,j).GT.1.) THEN
1888  var4(i,j) = min(xw4_all/xnsum_all/var(i,j) **4,10.)
1889  ENDIF
1890  ELSE
1891  print*, "FATAL ERROR: no source points in MAKEMT2."
1892  call abort()
1893  ENDIF
1894 
1895 ! set orog to 0 meters at ocean.
1896 ! IF (LAKE_FRAC(I,J) .EQ. 0. .AND. SLM(I,J) .EQ. 0.)THEN
1897  IF (lake_frac(i,j) .EQ. 0. .AND. land_frac(i,j) .EQ. 0.)THEN
1898  oro(i,j) = 0.0
1899  ENDIF
1900 
1901  ENDDO
1902  ENDDO
1903 !$omp end parallel do
1904  WRITE(6,*) "! MAKEMT2 ORO SLM VAR VAR4 DONE"
1905 C
1906  deallocate(hgt_1d)
1907  deallocate(hgt_1d_all)
1908  RETURN
1909  END
1910 
1939  SUBROUTINE makepc(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
1940  1 GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
1941 C
1942 C=== PC: principal coordinates of each Z avg orog box for L&M
1943 C
1944  parameter(rearth=6.3712e+6)
1945  dimension glat(jmn),xlat(jm),deltax(jmn)
1946  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
1947  DIMENSION ORO(IM,JM),SLM(IM,JM),HL(IM,JM),HK(IM,JM)
1948  DIMENSION HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
1949  DIMENSION THETA(IM,JM),GAMMA(IM,JM),SIGMA2(IM,JM),SIGMA(IM,JM)
1950  dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
1951  LOGICAL FLAG, DEBUG
1952 C=== DATA DEBUG/.TRUE./
1953  DATA debug/.false./
1954 C
1955  pi = 4.0 * atan(1.0)
1956  certh = pi * rearth
1957 C---- GLOBAL XLAT AND XLON ( DEGREE )
1958 C
1959  jm1 = jm - 1
1960  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1961  deltay = certh / float(jmn)
1962  print *, 'MAKEPC: DELTAY=',deltay
1963 C
1964  DO j=1,jmn
1965  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1966  deltax(j) = deltay * cos(glat(j)*pi/180.0)
1967  ENDDO
1968 C
1969 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1970 C
1971  DO j=1,jm
1972  DO i=1,numi(j)
1973 C IM1 = numi(j) - 1
1974  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
1975  faclon = delx / delxn
1976  ist(i,j) = faclon * float(i-1) - faclon * 0.5
1977  ist(i,j) = ist(i,j) + 1
1978  ien(i,j) = faclon * float(i) - faclon * 0.5
1979 C if (debug) then
1980 C if ( I .lt. 10 .and. J .lt. 10 )
1981 C 1 PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
1982 C endif
1983 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
1984 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
1985  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
1986  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
1987  if (debug) then
1988  if ( i .lt. 10 .and. j .lt. 10 )
1989  1 print*, ' I j IST IEN ',i,j,ist(i,j),ien(i,j)
1990  endif
1991  IF (ien(i,j) .LT. ist(i,j))
1992  1 print *,' MAKEPC: IEN < IST: I,J,IST(I,J),IEN(I,J)',
1993  2 i,j,ist(i,j),ien(i,j)
1994  ENDDO
1995  ENDDO
1996  DO j=1,jm-1
1997  flag=.true.
1998  DO j1=1,jmn
1999  xxlat = (xlat(j)+xlat(j+1))/2.
2000  IF(flag.AND.glat(j1).GT.xxlat) THEN
2001  jst(j) = j1
2002  jen(j+1) = j1 - 1
2003  flag = .false.
2004  ENDIF
2005  ENDDO
2006  ENDDO
2007  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2008  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2009  if (debug) then
2010  print*, ' IST,IEN(1,1-numi(1,JM))',ist(1,1),ien(1,1),
2011  1 ist(numi(jm),jm),ien(numi(jm),jm), numi(jm)
2012  print*, ' JST,JEN(1,JM) ',jst(1),jen(1),jst(jm),jen(jm)
2013  endif
2014 C
2015 C... DERIVITIVE TENSOR OF HEIGHT
2016 C
2017  DO j=1,jm
2018  DO i=1,numi(j)
2019  oro(i,j) = 0.0
2020  hx2(i,j) = 0.0
2021  hy2(i,j) = 0.0
2022  hxy(i,j) = 0.0
2023  xnsum = 0.0
2024  xland = 0.0
2025  xwatr = 0.0
2026  xl1 = 0.0
2027  xs1 = 0.0
2028  xfp = 0.0
2029  yfp = 0.0
2030  xfpyfp = 0.0
2031  xfp2 = 0.0
2032  yfp2 = 0.0
2033  hl(i,j) = 0.0
2034  hk(i,j) = 0.0
2035  hlprim(i,j) = 0.0
2036  theta(i,j) = 0.0
2037  gamma(i,j) = 0.
2038  sigma2(i,j) = 0.
2039  sigma(i,j) = 0.
2040 C
2041  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2042  i1 = ist(i,j) + ii1 - 1
2043  IF(i1.LE.0.) i1 = i1 + imn
2044  IF(i1.GT.imn) i1 = i1 - imn
2045 C
2046 C=== set the rest of the indexs for ave: 2pt staggered derivitive
2047 C
2048  i0 = i1 - 1
2049  if (i1 - 1 .le. 0 ) i0 = i0 + imn
2050  if (i1 - 1 .gt. imn) i0 = i0 - imn
2051 C
2052  ip1 = i1 + 1
2053  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2054  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2055 C
2056  DO j1=jst(j),jen(j)
2057  if (debug) then
2058  if ( i1 .eq. ist(i,j) .and. j1 .eq. jst(j) )
2059  1 print*, ' J, J1,IST,JST,DELTAX,GLAT ',
2060  2 j,j1,ist(i,j),jst(j),deltax(j1),glat(j1)
2061  if ( i1 .eq. ien(i,j) .and. j1 .eq. jen(j) )
2062  1 print*, ' J, J1,IEN,JEN,DELTAX,GLAT ',
2063  2 j,j1,ien(i,j),jen(j),deltax(j1),glat(j1)
2064  endif
2065  xland = xland + float(zslm(i1,j1))
2066  xwatr = xwatr + float(1-zslm(i1,j1))
2067  xnsum = xnsum + 1.
2068 C
2069  height = float(zavg(i1,j1))
2070  hi0 = float(zavg(i0,j1))
2071  hip1 = float(zavg(ip1,j1))
2072 C
2073  IF(height.LT.-990.) height = 0.0
2074  if(hi0 .lt. -990.) hi0 = 0.0
2075  if(hip1 .lt. -990.) hip1 = 0.0
2076 C........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
2077  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2078  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2079 C
2080 ! --- not at boundaries
2081 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
2082  if ( j1 .ne. jst(jm) .and. j1 .ne. jen(1) ) then
2083  hj0 = float(zavg(i1,j1-1))
2084  hjp1 = float(zavg(i1,j1+1))
2085  if(hj0 .lt. -990.) hj0 = 0.0
2086  if(hjp1 .lt. -990.) hjp1 = 0.0
2087 C....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
2088  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2089  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2090 C
2091 C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
2092 C === the NH pole: NB J1 goes from High at NP to Low toward SP
2093 C
2094 !RAB elseif ( J1 .eq. JST(1) ) then
2095  elseif ( j1 .eq. jst(jm) ) then
2096  ijax = i1 + imn/2
2097  if (ijax .le. 0 ) ijax = ijax + imn
2098  if (ijax .gt. imn) ijax = ijax - imn
2099 C..... at N pole we stay at the same latitude j1 but cross to opp side
2100  hijax = float(zavg(ijax,j1))
2101  hi1j1 = float(zavg(i1,j1))
2102  if(hijax .lt. -990.) hijax = 0.0
2103  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2104 C....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
2105  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2106  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2107  1 / deltay )**2
2108 C
2109 C === the SH pole: NB J1 goes from High at NP to Low toward SP
2110 C
2111 !RAB elseif ( J1 .eq. JEN(JM) ) then
2112  elseif ( j1 .eq. jen(1) ) then
2113  ijax = i1 + imn/2
2114  if (ijax .le. 0 ) ijax = ijax + imn
2115  if (ijax .gt. imn) ijax = ijax - imn
2116  hijax = float(zavg(ijax,j1))
2117  hi1j1 = float(zavg(i1,j1))
2118  if(hijax .lt. -990.) hijax = 0.0
2119  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2120  if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,hijax,hi1j1
2121 C..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
2122  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2123  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2124  1 / deltay )**2
2125  endif
2126 C
2127 C === The above does an average across the pole for the bndry in j.
2128 C23456789012345678901234567890123456789012345678901234567890123456789012......
2129 C
2130  xfpyfp = xfpyfp + xfp * yfp
2131  xl1 = xl1 + height * float(zslm(i1,j1))
2132  xs1 = xs1 + height * float(1-zslm(i1,j1))
2133 C
2134 C === average the HX2, HY2 and HXY
2135 C === This will be done over all land
2136 C
2137  ENDDO
2138  ENDDO
2139 C
2140 C === HTENSR
2141 C
2142  IF(xnsum.GT.1.) THEN
2143  slm(i,j) = float(nint(xland/xnsum))
2144  IF(slm(i,j).NE.0.) THEN
2145  oro(i,j)= xl1 / xland
2146  hx2(i,j) = xfp2 / xland
2147  hy2(i,j) = yfp2 / xland
2148  hxy(i,j) = xfpyfp / xland
2149  ELSE
2150  oro(i,j)= xs1 / xwatr
2151  ENDIF
2152 C=== degub testing
2153  if (debug) then
2154  print *," I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2155  1 xland,slm(i,j)
2156  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2157  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2158  ENDIF
2159 C
2160 C === make the principal axes, theta, and the degree of anisotropy,
2161 C === and sigma2, the slope parameter
2162 C
2163  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2164  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2165  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2166  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
2167 C
2168  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) * 180.0 / pi
2169 C === for testing print out in degrees
2170 C THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
2171  ENDIF
2172  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2173  if ( sigma2(i,j) .GE. 0. ) then
2174  sigma(i,j) = sqrt(sigma2(i,j) )
2175  if (sigma2(i,j) .ne. 0. .and.
2176  & hk(i,j) .GE. hlprim(i,j) )
2177  1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2178  else
2179  sigma(i,j)=0.
2180  endif
2181  ENDIF
2182  if (debug) then
2183  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2184  1 sigma(i,j),gamma(i,j)
2185  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2186  endif
2187  ENDDO
2188  ENDDO
2189  WRITE(6,*) "! MAKE Principal Coord DONE"
2190 C
2191  RETURN
2192  END
2193 
2214  SUBROUTINE makepc2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2215  1 GLAT,IM,JM,IMN,JMN,lon_c,lat_c,SLM)
2216 C
2217 C=== PC: principal coordinates of each Z avg orog box for L&M
2218 C
2219  implicit none
2220  real, parameter :: REARTH=6.3712e+6
2221  real, parameter :: D2R = 3.14159265358979/180.
2222  integer :: im,jm,imn,jmn
2223  real :: GLAT(JMN),DELTAX(JMN)
2224  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2225  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
2226  real, intent(in) :: SLM(IM,JM)
2227  real HL(IM,JM),HK(IM,JM)
2228  real HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
2229  real THETA(IM,JM),GAMMA(IM,JM),SIGMA2(IM,JM),SIGMA(IM,JM)
2230  real PI,CERTH,DELXN,DELTAY,XNSUM,XLAND
2231  real xfp,yfp,xfpyfp,xfp2,yfp2
2232  real hi0,hip1,hj0,hjp1,hijax,hi1j1
2233  real LONO(4),LATO(4),LONI,LATI
2234  integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
2235  integer ilist(IMN)
2236  logical inside_a_polygon
2237  LOGICAL DEBUG
2238 C=== DATA DEBUG/.TRUE./
2239  DATA debug/.false./
2240 C
2241  pi = 4.0 * atan(1.0)
2242  certh = pi * rearth
2243 C---- GLOBAL XLAT AND XLON ( DEGREE )
2244 C
2245  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2246  deltay = certh / float(jmn)
2247  print *, 'MAKEPC2: DELTAY=',deltay
2248 C
2249  DO j=1,jmn
2250  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2251  deltax(j) = deltay * cos(glat(j)*d2r)
2252  ENDDO
2253 C
2254 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2255 C
2256 
2257 C... DERIVITIVE TENSOR OF HEIGHT
2258 C
2259 !$omp parallel do
2260 !$omp* private (j,i,xnsum,xland,xfp,yfp,xfpyfp,
2261 !$omp* xfp2,yfp2,lono,lato,jst,jen,ilist,numx,j1,i2,i1,
2262 !$omp* loni,lati,i0,ip1,hi0,hip1,hj0,hjp1,ijax,
2263 !$omp* hijax,hi1j1)
2264  jloop : DO j=1,jm
2265 ! print*, "J=", J
2266  iloop : DO i=1,im
2267  hx2(i,j) = 0.0
2268  hy2(i,j) = 0.0
2269  hxy(i,j) = 0.0
2270  xnsum = 0.0
2271  xland = 0.0
2272  xfp = 0.0
2273  yfp = 0.0
2274  xfpyfp = 0.0
2275  xfp2 = 0.0
2276  yfp2 = 0.0
2277  hl(i,j) = 0.0
2278  hk(i,j) = 0.0
2279  hlprim(i,j) = 0.0
2280  theta(i,j) = 0.0
2281  gamma(i,j) = 0.
2282  sigma2(i,j) = 0.
2283  sigma(i,j) = 0.
2284 
2285  lono(1) = lon_c(i,j)
2286  lono(2) = lon_c(i+1,j)
2287  lono(3) = lon_c(i+1,j+1)
2288  lono(4) = lon_c(i,j+1)
2289  lato(1) = lat_c(i,j)
2290  lato(2) = lat_c(i+1,j)
2291  lato(3) = lat_c(i+1,j+1)
2292  lato(4) = lat_c(i,j+1)
2293  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
2294 
2295  do j1 = jst, jen; do i2 = 1, numx
2296  i1 = ilist(i2)
2297  loni = i1*delxn
2298  lati = -90 + j1*delxn
2299  inside : if(inside_a_polygon(loni*d2r,lati*d2r,4,
2300  & lono*d2r,lato*d2r))then
2301 
2302 C=== set the rest of the indexs for ave: 2pt staggered derivitive
2303 C
2304  i0 = i1 - 1
2305  if (i1 - 1 .le. 0 ) i0 = i0 + imn
2306  if (i1 - 1 .gt. imn) i0 = i0 - imn
2307 C
2308  ip1 = i1 + 1
2309  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2310  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2311 
2312  xland = xland + float(zslm(i1,j1))
2313  xnsum = xnsum + 1.
2314 C
2315  hi0 = float(zavg(i0,j1))
2316  hip1 = float(zavg(ip1,j1))
2317 C
2318  if(hi0 .lt. -990.) hi0 = 0.0
2319  if(hip1 .lt. -990.) hip1 = 0.0
2320 C........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
2321  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2322  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2323 C
2324 ! --- not at boundaries
2325 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
2326  if ( j1 .ne. 1 .and. j1 .ne. jmn ) then
2327  hj0 = float(zavg(i1,j1-1))
2328  hjp1 = float(zavg(i1,j1+1))
2329  if(hj0 .lt. -990.) hj0 = 0.0
2330  if(hjp1 .lt. -990.) hjp1 = 0.0
2331 C....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
2332  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2333  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2334 C
2335 C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
2336 C === the NH pole: NB J1 goes from High at NP to Low toward SP
2337 C
2338 !RAB elseif ( J1 .eq. JST(1) ) then
2339  elseif ( j1 .eq. 1 ) then
2340  ijax = i1 + imn/2
2341  if (ijax .le. 0 ) ijax = ijax + imn
2342  if (ijax .gt. imn) ijax = ijax - imn
2343 C..... at N pole we stay at the same latitude j1 but cross to opp side
2344  hijax = float(zavg(ijax,j1))
2345  hi1j1 = float(zavg(i1,j1))
2346  if(hijax .lt. -990.) hijax = 0.0
2347  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2348 C....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
2349  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2350  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2351  1 / deltay )**2
2352 C
2353 C === the SH pole: NB J1 goes from High at NP to Low toward SP
2354 C
2355 !RAB elseif ( J1 .eq. JEN(JM) ) then
2356  elseif ( j1 .eq. jmn ) then
2357  ijax = i1 + imn/2
2358  if (ijax .le. 0 ) ijax = ijax + imn
2359  if (ijax .gt. imn) ijax = ijax - imn
2360  hijax = float(zavg(ijax,j1))
2361  hi1j1 = float(zavg(i1,j1))
2362  if(hijax .lt. -990.) hijax = 0.0
2363  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2364  if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,
2365  & hijax,hi1j1
2366 C..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
2367  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2368  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2369  1 / deltay )**2
2370  endif
2371 C
2372 C === The above does an average across the pole for the bndry in j.
2373 C23456789012345678901234567890123456789012345678901234567890123456789012......
2374 C
2375  xfpyfp = xfpyfp + xfp * yfp
2376  ENDIF inside
2377 C
2378 C === average the HX2, HY2 and HXY
2379 C === This will be done over all land
2380 C
2381  ENDDO
2382  ENDDO
2383 C
2384 C === HTENSR
2385 C
2386  xnsum_gt_1 : IF(xnsum.GT.1.) THEN
2387  IF(slm(i,j).NE.0.) THEN
2388  IF (xland > 0) THEN
2389  hx2(i,j) = xfp2 / xland
2390  hy2(i,j) = yfp2 / xland
2391  hxy(i,j) = xfpyfp / xland
2392  ELSE
2393  hx2(i,j) = xfp2 / xnsum
2394  hy2(i,j) = yfp2 / xnsum
2395  hxy(i,j) = xfpyfp / xnsum
2396  ENDIF
2397  ENDIF
2398 C=== degub testing
2399  if (debug) then
2400  print *," I,J,i1,j1:", i,j,i1,j1,
2401  1 xland,slm(i,j)
2402  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2403  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2404  ENDIF
2405 C
2406 C === make the principal axes, theta, and the degree of anisotropy,
2407 C === and sigma2, the slope parameter
2408 C
2409  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2410  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2411  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2412  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
2413 C
2414  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) / d2r
2415 C === for testing print out in degrees
2416 C THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
2417  ENDIF
2418  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2419  if ( sigma2(i,j) .GE. 0. ) then
2420  sigma(i,j) = sqrt(sigma2(i,j) )
2421  if (sigma2(i,j) .ne. 0. .and.
2422  & hk(i,j) .GE. hlprim(i,j) )
2423  1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2424  else
2425  sigma(i,j)=0.
2426  endif
2427  ENDIF xnsum_gt_1
2428  if (debug) then
2429  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2430  1 sigma(i,j),gamma(i,j)
2431  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2432  endif
2433  ENDDO iloop
2434  ENDDO jloop
2435 !$omp end parallel do
2436  WRITE(6,*) "! MAKE Principal Coord DONE"
2437 C
2438  RETURN
2439  END SUBROUTINE makepc2
2440 
2482  SUBROUTINE makeoa(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2483  1 ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
2484  2 IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
2485  dimension glat(jmn),xlat(jm)
2486  INTEGER ZAVG(IMN,JMN)
2487  DIMENSION ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
2488  DIMENSION OA4(IM,JM,4),IOA4(IM,JM,4)
2489  DIMENSION IST(IM,jm),IEN(IM,jm),JST(JM),JEN(JM)
2490  DIMENSION XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
2491  dimension xnsum3(im,jm),xnsum4(im,jm)
2492  dimension var(im,jm),ol(im,jm,4),numi(jm)
2493  LOGICAL FLAG
2494 C
2495 C---- GLOBAL XLAT AND XLON ( DEGREE )
2496 C
2497 ! --- IM1 = IM - 1 removed (not used in this sub)
2498  jm1 = jm - 1
2499  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2500 C
2501  DO j=1,jmn
2502  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2503  ENDDO
2504  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
2505 C
2506 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2507 C
2508  DO j=1,jm
2509  DO i=1,numi(j)
2510  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
2511  faclon = delx / delxn
2512 C --- minus sign here in IST and IEN as in MAKEMT!
2513  ist(i,j) = faclon * float(i-1) - faclon * 0.5
2514  ist(i,j) = ist(i,j) + 1
2515  ien(i,j) = faclon * float(i) - faclon * 0.5
2516 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
2517 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
2518  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2519  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2520 cx PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
2521  if ( i .lt. 3 .and. j .lt. 3 )
2522  1print*,' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2523  if ( i .lt. 3 .and. j .ge. jm-1 )
2524  1print*,' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2525  ENDDO
2526  ENDDO
2527  print *,'MAKEOA: DELXN,DELX,FACLON',delxn,delx,faclon
2528  print *, ' ***** ready to start JST JEN section '
2529  DO j=1,jm-1
2530  flag=.true.
2531  DO j1=1,jmn
2532 ! --- XXLAT added as in MAKEMT and in next line as well
2533  xxlat = (xlat(j)+xlat(j+1))/2.
2534  IF(flag.AND.glat(j1).GT.xxlat) THEN
2535  jst(j) = j1
2536 ! --- JEN(J+1) = J1 - 1
2537  flag = .false.
2538  if ( j .eq. 1 )
2539  1print*,' MAKEOA: XX j JST JEN ',j,jst(j),jen(j)
2540  ENDIF
2541  ENDDO
2542  if ( j .lt. 3 )
2543  1print*,' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2544  if ( j .ge. jm-2 )
2545  1print*,' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2546 C FLAG=.TRUE.
2547 C DO J1=JST(J),JMN
2548 C IF(FLAG.AND.GLAT(J1).GT.XLAT(J)) THEN
2549 C JEN(J) = J1 - 1
2550 C FLAG = .FALSE.
2551 C ENDIF
2552 C ENDDO
2553  ENDDO
2554  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2555  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2556  print *,' ***** JST(1) JEN(1) ',jst(1),jen(1)
2557  print *,' ***** JST(JM) JEN(JM) ',jst(jm),jen(jm)
2558 C
2559  DO j=1,jm
2560  DO i=1,numi(j)
2561  xnsum(i,j) = 0.0
2562  elvmax(i,j) = oro(i,j)
2563  zmax(i,j) = 0.0
2564  ENDDO
2565  ENDDO
2566 !
2567 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
2568 ! --- to JM or to JM1
2569  DO j=1,jm
2570  DO i=1,numi(j)
2571  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2572  i1 = ist(i,j) + ii1 - 1
2573 ! --- next line as in makemt (I1 not II1) (*j*) 20070701
2574  IF(i1.LE.0.) i1 = i1 + imn
2575  IF (i1 .GT. imn) i1 = i1 - imn
2576  DO j1=jst(j),jen(j)
2577  height = float(zavg(i1,j1))
2578  IF(height.LT.-990.) height = 0.0
2579  IF ( height .gt. oro(i,j) ) then
2580  if ( height .gt. zmax(i,j) )zmax(i,j) = height
2581  xnsum(i,j) = xnsum(i,j) + 1
2582  ENDIF
2583  ENDDO
2584  ENDDO
2585  if ( i .lt. 5 .and. j .ge. jm-5 ) then
2586  print *,' I,J,ORO(I,J),XNSUM(I,J),ZMAX(I,J):',
2587  1 i,j,oro(i,j),xnsum(i,j),zmax(i,j)
2588  endif
2589  ENDDO
2590  ENDDO
2591 !
2592 C.... make ELVMAX ORO from MAKEMT sub
2593 C
2594 ! --- this will make work1 array take on oro's values on return
2595  DO j=1,jm
2596  DO i=1,numi(j)
2597 
2598  oro1(i,j) = oro(i,j)
2599  elvmax(i,j) = zmax(i,j)
2600  ENDDO
2601  ENDDO
2602 C........
2603 C The MAX elev peak (no averaging)
2604 C........
2605 ! DO J=1,JM
2606 ! DO I=1,numi(j)
2607 ! DO II1 = 1, IEN(I,J) - IST(I,J) + 1
2608 ! I1 = IST(I,J) + II1 - 1
2609 ! IF(I1.LE.0.) I1 = I1 + IMN
2610 ! IF(I1.GT.IMN) I1 = I1 - IMN
2611 ! DO J1=JST(J),JEN(J)
2612 ! if ( ELVMAX(I,J) .lt. ZMAX(I1,J1))
2613 ! 1 ELVMAX(I,J) = ZMAX(I1,J1)
2614 ! ENDDO
2615 ! ENDDO
2616 ! ENDDO
2617 ! ENDDO
2618 C
2619 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
2620 C IN A GRID BOX
2621  DO j=1,jm
2622  DO i=1,numi(j)
2623  xnsum1(i,j) = 0.0
2624  xnsum2(i,j) = 0.0
2625  xnsum3(i,j) = 0.0
2626  xnsum4(i,j) = 0.0
2627  ENDDO
2628  ENDDO
2629 ! --- loop
2630  DO j=1,jm1
2631  DO i=1,numi(j)
2632  hc = 1116.2 - 0.878 * var(i,j)
2633 ! print *,' I,J,HC,VAR:',I,J,HC,VAR(I,J)
2634  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2635  i1 = ist(i,j) + ii1 - 1
2636 ! IF (I1.LE.0.) print *,' I1 less than 0',I1,II1,IST(I,J),IEN(I,J)
2637 ! if ( J .lt. 3 .or. J .gt. JM-2 ) then
2638 ! IF(I1 .GT. IMN)print *,' I1 > IMN',J,I1,II1,IMN,IST(I,J),IEN(I,J)
2639 ! endif
2640  IF(i1.GT.imn) i1 = i1 - imn
2641  DO j1=jst(j),jen(j)
2642  IF(float(zavg(i1,j1)) .GT. hc)
2643  1 xnsum1(i,j) = xnsum1(i,j) + 1
2644  xnsum2(i,j) = xnsum2(i,j) + 1
2645  ENDDO
2646  ENDDO
2647 C
2648  inci = nint((ien(i,j)-ist(i,j)) * 0.5)
2649  isttt = min(max(ist(i,j)-inci,1),imn)
2650  ieddd = min(max(ien(i,j)-inci,1),imn)
2651 C
2652  incj = nint((jen(j)-jst(j)) * 0.5)
2653  jsttt = min(max(jst(j)-incj,1),jmn)
2654  jeddd = min(max(jen(j)-incj,1),jmn)
2655 ! if ( J .lt. 3 .or. J .gt. JM-3 ) then
2656 ! if(I .lt. 3 .or. I .gt. IM-3) then
2657 ! print *,' INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD:',
2658 ! 1 I,J,INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD
2659 ! endif
2660 ! endif
2661 C
2662  DO i1=isttt,ieddd
2663  DO j1=jsttt,jeddd
2664  IF(float(zavg(i1,j1)) .GT. hc)
2665  1 xnsum3(i,j) = xnsum3(i,j) + 1
2666  xnsum4(i,j) = xnsum4(i,j) + 1
2667  ENDDO
2668  ENDDO
2669 cx print*,' i j hc var ',i,j,hc,var(i,j)
2670 cx print*,'xnsum12 ',xnsum1(i,j),xnsum2(i,j)
2671 cx print*,'xnsum34 ',xnsum3(i,j),xnsum4(i,j)
2672  ENDDO
2673  ENDDO
2674 C
2675 C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
2676 C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
2677 C (KWD = 1 2 3 4)
2678 C ( WD = W S SW NW)
2679 C
2680 C
2681  DO kwd = 1, 4
2682  DO j=1,jm
2683  DO i=1,numi(j)
2684  oa4(i,j,kwd) = 0.0
2685  ENDDO
2686  ENDDO
2687  ENDDO
2688 C
2689  DO j=1,jm-2
2690  DO i=1,numi(j)
2691  ii = i + 1
2692  IF (ii .GT. numi(j)) ii = ii - numi(j)
2693  xnpu = xnsum(i,j) + xnsum(i,j+1)
2694  xnpd = xnsum(ii,j) + xnsum(ii,j+1)
2695  IF (xnpd .NE. xnpu) oa4(ii,j+1,1) = 1. - xnpd / max(xnpu , 1.)
2696  ol(ii,j+1,1) = (xnsum3(i,j+1)+xnsum3(ii,j+1))/
2697  1 (xnsum4(i,j+1)+xnsum4(ii,j+1))
2698 ! if ( I .lt. 20 .and. J .ge. JM-19 ) then
2699 ! PRINT*,' MAKEOA: I J IST IEN ',I,j,IST(I,J),IEN(I,J)
2700 ! PRINT*,' HC VAR ',HC,VAR(i,j)
2701 ! PRINT*,' MAKEOA: XNSUM(I,J)=',XNSUM(I,J),XNPU, XNPD
2702 ! PRINT*,' MAKEOA: XNSUM3(I,J+1),XNSUM3(II,J+1)',
2703 ! 1 XNSUM3(I,J+1),XNSUM3(II,J+1)
2704 ! PRINT*,' MAKEOA: II, OA4(II,J+1,1), OL(II,J+1,1):',
2705 ! 1 II, OA4(II,J+1,1), OL(II,J+1,1)
2706 ! endif
2707  ENDDO
2708  ENDDO
2709  DO j=1,jm-2
2710  DO i=1,numi(j)
2711  ii = i + 1
2712  IF (ii .GT. numi(j)) ii = ii - numi(j)
2713  xnpu = xnsum(i,j+1) + xnsum(ii,j+1)
2714  xnpd = xnsum(i,j) + xnsum(ii,j)
2715  IF (xnpd .NE. xnpu) oa4(ii,j+1,2) = 1. - xnpd / max(xnpu , 1.)
2716  ol(ii,j+1,2) = (xnsum3(ii,j)+xnsum3(ii,j+1))/
2717  1 (xnsum4(ii,j)+xnsum4(ii,j+1))
2718  ENDDO
2719  ENDDO
2720  DO j=1,jm-2
2721  DO i=1,numi(j)
2722  ii = i + 1
2723  IF (ii .GT. numi(j)) ii = ii - numi(j)
2724  xnpu = xnsum(i,j+1) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2725  xnpd = xnsum(ii,j) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2726  IF (xnpd .NE. xnpu) oa4(ii,j+1,3) = 1. - xnpd / max(xnpu , 1.)
2727  ol(ii,j+1,3) = (xnsum1(ii,j)+xnsum1(i,j+1))/
2728  1 (xnsum2(ii,j)+xnsum2(i,j+1))
2729  ENDDO
2730  ENDDO
2731  DO j=1,jm-2
2732  DO i=1,numi(j)
2733  ii = i + 1
2734  IF (ii .GT. numi(j)) ii = ii - numi(j)
2735  xnpu = xnsum(i,j) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2736  xnpd = xnsum(ii,j+1) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2737  IF (xnpd .NE. xnpu) oa4(ii,j+1,4) = 1. - xnpd / max(xnpu , 1.)
2738  ol(ii,j+1,4) = (xnsum1(i,j)+xnsum1(ii,j+1))/
2739  1 (xnsum2(i,j)+xnsum2(ii,j+1))
2740  ENDDO
2741  ENDDO
2742 C
2743  DO kwd = 1, 4
2744  DO i=1,numi(j)
2745  ol(i,1,kwd) = ol(i,2,kwd)
2746  ol(i,jm,kwd) = ol(i,jm-1,kwd)
2747  ENDDO
2748  ENDDO
2749 C
2750  DO kwd=1,4
2751  DO j=1,jm
2752  DO i=1,numi(j)
2753  t = oa4(i,j,kwd)
2754  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
2755  ENDDO
2756  ENDDO
2757  ENDDO
2758 C
2759  ns0 = 0
2760  ns1 = 0
2761  ns2 = 0
2762  ns3 = 0
2763  ns4 = 0
2764  ns5 = 0
2765  ns6 = 0
2766  DO kwd=1,4
2767  DO j=1,jm
2768  DO i=1,numi(j)
2769  t = abs( oa4(i,j,kwd) )
2770  IF(t .EQ. 0.) THEN
2771  ioa4(i,j,kwd) = 0
2772  ns0 = ns0 + 1
2773  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
2774  ioa4(i,j,kwd) = 1
2775  ns1 = ns1 + 1
2776  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
2777  ioa4(i,j,kwd) = 2
2778  ns2 = ns2 + 1
2779  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
2780  ioa4(i,j,kwd) = 3
2781  ns3 = ns3 + 1
2782  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
2783  ioa4(i,j,kwd) = 4
2784  ns4 = ns4 + 1
2785  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
2786  ioa4(i,j,kwd) = 5
2787  ns5 = ns5 + 1
2788  ELSE IF(t .GT. 10000.) THEN
2789  ioa4(i,j,kwd) = 6
2790  ns6 = ns6 + 1
2791  ENDIF
2792  ENDDO
2793  ENDDO
2794  ENDDO
2795 C
2796  WRITE(6,*) "! MAKEOA EXIT"
2797 C
2798  RETURN
2799  END SUBROUTINE makeoa
2800 
2810  function get_lon_angle(dx,lat, DEGRAD)
2811  implicit none
2812  real dx, lat, degrad
2813 
2814  real get_lon_angle
2815  real, parameter :: radius = 6371200
2816 
2817  get_lon_angle = 2*asin( sin(dx/radius*0.5)/cos(lat) )*degrad
2818 
2819  end function get_lon_angle
2820 
2829  function get_lat_angle(dy, DEGRAD)
2830  implicit none
2831  real dy, degrad
2832 
2833  real get_lat_angle
2834  real, parameter :: radius = 6371200
2835 
2836  get_lat_angle = dy/radius*degrad
2837 
2838  end function get_lat_angle
2839 
2876  SUBROUTINE makeoa2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2877  1 ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
2878  2 IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,dx,dy,
2879  3 is_south_pole,is_north_pole )
2880  implicit none
2881  real, parameter :: MISSING_VALUE = -9999.
2882  real, parameter :: d2r = 3.14159265358979/180.
2883  real, PARAMETER :: R2D=180./3.14159265358979
2884  integer im,jm,imn,jmn
2885  real GLAT(JMN)
2886  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2887  real ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
2888  real OA4(IM,JM,4)
2889  integer IOA4(IM,JM,4)
2890  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
2891  real lon_t(IM,JM), lat_t(IM,JM)
2892  real dx(IM,JM), dy(IM,JM)
2893  logical is_south_pole(IM,JM), is_north_pole(IM,JM)
2894  real XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
2895  real XNSUM3(IM,JM),XNSUM4(IM,JM)
2896  real VAR(IM,JM),OL(IM,JM,4)
2897  integer i,j,ilist(IMN),numx,i1,j1,ii1
2898  integer KWD
2899  real LONO(4),LATO(4),LONI,LATI
2900  real DELXN,HC,HEIGHT,XNPU,XNPD,T
2901  integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
2902  logical inside_a_polygon
2903  real lon,lat,dlon,dlat,dlat_old
2904  real lon1,lat1,lon2,lat2
2905  real xnsum11,xnsum12,xnsum21,xnsum22
2906  real HC_11, HC_12, HC_21, HC_22
2907  real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
2908  real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
2909  real get_lon_angle, get_lat_angle, get_xnsum
2910  integer jst, jen
2911 C
2912 C---- GLOBAL XLAT AND XLON ( DEGREE )
2913 C
2914  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2915 C
2916  DO j=1,jmn
2917  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2918  ENDDO
2919  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
2920 C
2921 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2922 C
2923 C
2924  DO j=1,jm
2925  DO i=1,im
2926  xnsum(i,j) = 0.0
2927  elvmax(i,j) = oro(i,j)
2928  zmax(i,j) = 0.0
2929 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
2930 C IN A GRID BOX
2931  xnsum1(i,j) = 0.0
2932  xnsum2(i,j) = 0.0
2933  xnsum3(i,j) = 0.0
2934  xnsum4(i,j) = 0.0
2935  oro1(i,j) = oro(i,j)
2936  elvmax(i,j) = zmax(i,j)
2937  ENDDO
2938  ENDDO
2939 
2940 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
2941 ! --- to JM or to JM1
2942 !$omp parallel do
2943 !$omp* private (j,i,hc,lono,lato,jst,jen,ilist,numx,j1,ii1,i1,loni,
2944 !$omp* lati,height)
2945  DO j=1,jm
2946 ! print*, "J=", J
2947  DO i=1,im
2948  hc = 1116.2 - 0.878 * var(i,j)
2949  lono(1) = lon_c(i,j)
2950  lono(2) = lon_c(i+1,j)
2951  lono(3) = lon_c(i+1,j+1)
2952  lono(4) = lon_c(i,j+1)
2953  lato(1) = lat_c(i,j)
2954  lato(2) = lat_c(i+1,j)
2955  lato(3) = lat_c(i+1,j+1)
2956  lato(4) = lat_c(i,j+1)
2957  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
2958  do j1 = jst, jen; do ii1 = 1, numx
2959  i1 = ilist(ii1)
2960  loni = i1*delxn
2961  lati = -90 + j1*delxn
2962  if(inside_a_polygon(loni*d2r,lati*d2r,4,
2963  & lono*d2r,lato*d2r))then
2964 
2965  height = float(zavg(i1,j1))
2966  IF(height.LT.-990.) height = 0.0
2967  IF ( height .gt. oro(i,j) ) then
2968  if ( height .gt. zmax(i,j) )zmax(i,j) = height
2969  ENDIF
2970  endif
2971  ENDDO ; ENDDO
2972  ENDDO
2973  ENDDO
2974 !$omp end parallel do
2975 C
2976 ! --- this will make work1 array take on oro's values on return
2977 ! --- this will make work1 array take on oro's values on return
2978  DO j=1,jm
2979  DO i=1,im
2980 
2981  oro1(i,j) = oro(i,j)
2982  elvmax(i,j) = zmax(i,j)
2983  ENDDO
2984  ENDDO
2985 
2986  DO kwd = 1, 4
2987  DO j=1,jm
2988  DO i=1,im
2989  oa4(i,j,kwd) = 0.0
2990  ol(i,j,kwd) = 0.0
2991  ENDDO
2992  ENDDO
2993  ENDDO
2994  !
2995 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
2996 C
2997 C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
2998 C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
2999 C (KWD = 1 2 3 4)
3000 C ( WD = W S SW NW)
3001 C
3002 C
3003 !$omp parallel do
3004 !$omp* private (j,i,lon,lat,kwd,dlon,dlat,lon1,lon2,lat1,lat2,
3005 !$omp* xnsum11,xnsum12,xnsum21,xnsum22,xnpu,xnpd,
3006 !$omp* xnsum1_11,xnsum2_11,hc_11, xnsum1_12,xnsum2_12,
3007 !$omp* hc_12,xnsum1_21,xnsum2_21,hc_21, xnsum1_22,
3008 !$omp* xnsum2_22,hc_22)
3009  DO j=1,jm
3010 ! print*, "j = ", j
3011  DO i=1,im
3012  lon = lon_t(i,j)
3013  lat = lat_t(i,j)
3014  !--- for around north pole, oa and ol are all 0
3015 
3016  if(is_north_pole(i,j)) then
3017  print*, "set oa1 = 0 and ol=0 at i,j=", i,j
3018  do kwd = 1, 4
3019  oa4(i,j,kwd) = 0.
3020  ol(i,j,kwd) = 0.
3021  enddo
3022  else if(is_south_pole(i,j)) then
3023  print*, "set oa1 = 0 and ol=1 at i,j=", i,j
3024  do kwd = 1, 4
3025  oa4(i,j,kwd) = 0.
3026  ol(i,j,kwd) = 1.
3027  enddo
3028  else
3029 
3030  !--- for each point, find a lat-lon grid box with same dx and dy as the cubic grid box
3031  dlon = get_lon_angle(dx(i,j), lat*d2r, r2d )
3032  dlat = get_lat_angle(dy(i,j), r2d)
3033  !--- adjust dlat if the points are close to pole.
3034  if( lat-dlat*0.5<-90.) then
3035  print*, "at i,j =", i,j, lat, dlat, lat-dlat*0.5
3036  print*, "FATAL ERROR: lat-dlat*0.5<-90."
3037  call errexit(4)
3038  endif
3039  if( lat+dlat*2 > 90.) then
3040  dlat_old = dlat
3041  dlat = (90-lat)*0.5
3042  print*, "at i,j=",i,j," adjust dlat from ",
3043  & dlat_old, " to ", dlat
3044  endif
3045  !--- lower left
3046  lon1 = lon-dlon*1.5
3047  lon2 = lon-dlon*0.5
3048  lat1 = lat-dlat*0.5
3049  lat2 = lat+dlat*0.5
3050 
3051  if(lat1<-90 .or. lat2>90) then
3052  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3053  endif
3054  xnsum11 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3055  & zavg,zslm,delxn)
3056 
3057  !--- upper left
3058  lon1 = lon-dlon*1.5
3059  lon2 = lon-dlon*0.5
3060  lat1 = lat+dlat*0.5
3061  lat2 = lat+dlat*1.5
3062  if(lat1<-90 .or. lat2>90) then
3063  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3064  endif
3065  xnsum12 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3066  & zavg,zslm,delxn)
3067 
3068  !--- lower right
3069  lon1 = lon-dlon*0.5
3070  lon2 = lon+dlon*0.5
3071  lat1 = lat-dlat*0.5
3072  lat2 = lat+dlat*0.5
3073  if(lat1<-90 .or. lat2>90) then
3074  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3075  endif
3076  xnsum21 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3077  & zavg,zslm,delxn)
3078 
3079  !--- upper right
3080  lon1 = lon-dlon*0.5
3081  lon2 = lon+dlon*0.5
3082  lat1 = lat+dlat*0.5
3083  lat2 = lat+dlat*1.5
3084  if(lat1<-90 .or. lat2>90) then
3085  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3086  endif
3087 
3088  xnsum22 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3089  & zavg,zslm,delxn)
3090 
3091  xnpu = xnsum11 + xnsum12
3092  xnpd = xnsum21 + xnsum22
3093  IF (xnpd .NE. xnpu) oa4(i,j,1) = 1. - xnpd / max(xnpu , 1.)
3094 
3095  xnpu = xnsum11 + xnsum21
3096  xnpd = xnsum12 + xnsum22
3097  IF (xnpd .NE. xnpu) oa4(i,j,2) = 1. - xnpd / max(xnpu , 1.)
3098 
3099  xnpu = xnsum11 + (xnsum12+xnsum21)*0.5
3100  xnpd = xnsum22 + (xnsum12+xnsum21)*0.5
3101  IF (xnpd .NE. xnpu) oa4(i,j,3) = 1. - xnpd / max(xnpu , 1.)
3102 
3103  xnpu = xnsum12 + (xnsum11+xnsum22)*0.5
3104  xnpd = xnsum21 + (xnsum11+xnsum22)*0.5
3105  IF (xnpd .NE. xnpu) oa4(i,j,4) = 1. - xnpd / max(xnpu , 1.)
3106 
3107 
3108  !--- calculate OL3 and OL4
3109  !--- lower left
3110  lon1 = lon-dlon*1.5
3111  lon2 = lon-dlon*0.5
3112  lat1 = lat-dlat*0.5
3113  lat2 = lat+dlat*0.5
3114  if(lat1<-90 .or. lat2>90) then
3115  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3116  endif
3117  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3118  & zavg,delxn, xnsum1_11, xnsum2_11, hc_11)
3119 
3120  !--- upper left
3121  lon1 = lon-dlon*1.5
3122  lon2 = lon-dlon*0.5
3123  lat1 = lat+dlat*0.5
3124  lat2 = lat+dlat*1.5
3125  if(lat1<-90 .or. lat2>90) then
3126  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3127  endif
3128  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3129  & zavg,delxn, xnsum1_12, xnsum2_12, hc_12)
3130 
3131  !--- lower right
3132  lon1 = lon-dlon*0.5
3133  lon2 = lon+dlon*0.5
3134  lat1 = lat-dlat*0.5
3135  lat2 = lat+dlat*0.5
3136  if(lat1<-90 .or. lat2>90) then
3137  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3138  endif
3139  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3140  & zavg,delxn, xnsum1_21, xnsum2_21, hc_21)
3141 
3142  !--- upper right
3143  lon1 = lon-dlon*0.5
3144  lon2 = lon+dlon*0.5
3145  lat1 = lat+dlat*0.5
3146  lat2 = lat+dlat*1.5
3147  if(lat1<-90 .or. lat2>90) then
3148  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3149  endif
3150  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3151  & zavg,delxn, xnsum1_22, xnsum2_22, hc_22)
3152 
3153  ol(i,j,3) = (xnsum1_22+xnsum1_11)/(xnsum2_22+xnsum2_11)
3154  ol(i,j,4) = (xnsum1_12+xnsum1_21)/(xnsum2_12+xnsum2_21)
3155 
3156  !--- calculate OL1 and OL2
3157  !--- lower left
3158  lon1 = lon-dlon*2.0
3159  lon2 = lon-dlon
3160  lat1 = lat
3161  lat2 = lat+dlat
3162  if(lat1<-90 .or. lat2>90) then
3163  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3164  endif
3165  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3166  & zavg,delxn, xnsum1_11, xnsum2_11, hc_11)
3167 
3168  !--- upper left
3169  lon1 = lon-dlon*2.0
3170  lon2 = lon-dlon
3171  lat1 = lat+dlat
3172  lat2 = lat+dlat*2.0
3173  if(lat1<-90 .or. lat2>90) then
3174  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3175  endif
3176 
3177  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3178  & zavg,delxn, xnsum1_12, xnsum2_12, hc_12)
3179 
3180  !--- lower right
3181  lon1 = lon-dlon
3182  lon2 = lon
3183  lat1 = lat
3184  lat2 = lat+dlat
3185  if(lat1<-90 .or. lat2>90) then
3186  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3187  endif
3188  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3189  & zavg,delxn, xnsum1_21, xnsum2_21, hc_21)
3190 
3191  !--- upper right
3192  lon1 = lon-dlon
3193  lon2 = lon
3194  lat1 = lat+dlat
3195  lat2 = lat+dlat*2.0
3196  if(lat1<-90 .or. lat2>90) then
3197  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3198  endif
3199 
3200  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3201  & zavg,delxn, xnsum1_22, xnsum2_22, hc_22)
3202 
3203  ol(i,j,1) = (xnsum1_11+xnsum1_21)/(xnsum2_11+xnsum2_21)
3204  ol(i,j,2) = (xnsum1_21+xnsum1_22)/(xnsum2_21+xnsum2_22)
3205  ENDIF
3206  ENDDO
3207  ENDDO
3208 !$omp end parallel do
3209  DO kwd=1,4
3210  DO j=1,jm
3211  DO i=1,im
3212  t = oa4(i,j,kwd)
3213  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3214  ENDDO
3215  ENDDO
3216  ENDDO
3217 C
3218  ns0 = 0
3219  ns1 = 0
3220  ns2 = 0
3221  ns3 = 0
3222  ns4 = 0
3223  ns5 = 0
3224  ns6 = 0
3225  DO kwd=1,4
3226  DO j=1,jm
3227  DO i=1,im
3228  t = abs( oa4(i,j,kwd) )
3229  IF(t .EQ. 0.) THEN
3230  ioa4(i,j,kwd) = 0
3231  ns0 = ns0 + 1
3232  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
3233  ioa4(i,j,kwd) = 1
3234  ns1 = ns1 + 1
3235  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
3236  ioa4(i,j,kwd) = 2
3237  ns2 = ns2 + 1
3238  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
3239  ioa4(i,j,kwd) = 3
3240  ns3 = ns3 + 1
3241  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
3242  ioa4(i,j,kwd) = 4
3243  ns4 = ns4 + 1
3244  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
3245  ioa4(i,j,kwd) = 5
3246  ns5 = ns5 + 1
3247  ELSE IF(t .GT. 10000.) THEN
3248  ioa4(i,j,kwd) = 6
3249  ns6 = ns6 + 1
3250  ENDIF
3251  ENDDO
3252  ENDDO
3253  ENDDO
3254 C
3255  WRITE(6,*) "! MAKEOA2 EXIT"
3256 C
3257  RETURN
3258 
3259  END SUBROUTINE makeoa2
3260 
3269  function spherical_distance(theta1,phi1,theta2,phi2)
3270 
3271  real, intent(in) :: theta1, phi1, theta2, phi2
3272  real :: spherical_distance, dot
3273 
3274  if(theta1 == theta2 .and. phi1 == phi2) then
3275  spherical_distance = 0.0
3276  return
3277  endif
3278 
3279  dot = cos(phi1)*cos(phi2)*cos(theta1-theta2) + sin(phi1)*sin(phi2)
3280  if(dot > 1. ) dot = 1.
3281  if(dot < -1.) dot = -1.
3282  spherical_distance = acos(dot)
3283 
3284  return
3285 
3286  end function spherical_distance
3287 
3304  subroutine get_mismatch_index(im_in, jm_in, geolon_in,geolat_in,
3305  & bitmap_in,num_out, lon_out,lat_out, iindx, jindx )
3306  integer, intent(in) :: im_in, jm_in, num_out
3307  real, intent(in) :: geolon_in(im_in,jm_in)
3308  real, intent(in) :: geolat_in(im_in,jm_in)
3309  logical*1, intent(in) :: bitmap_in(im_in,jm_in)
3310  real, intent(in) :: lon_out(num_out), lat_out(num_out)
3311  integer, intent(out):: iindx(num_out), jindx(num_out)
3312  real, parameter :: MAX_DIST = 1.e+20
3313  integer, parameter :: NUMNBR = 20
3314  integer :: i_c,j_c,jstart,jend
3315  real :: lon,lat
3316 
3317  print*, "im_in,jm_in = ", im_in, jm_in
3318  print*, "num_out = ", num_out
3319  print*, "max and min of lon_in is", minval(geolon_in),
3320  & maxval(geolon_in)
3321  print*, "max and min of lat_in is", minval(geolat_in),
3322  & maxval(geolat_in)
3323  print*, "max and min of lon_out is", minval(lon_out),
3324  & maxval(lon_out)
3325  print*, "max and min of lat_out is", minval(lat_out),
3326  & maxval(lat_out)
3327  print*, "count(bitmap_in)= ", count(bitmap_in), max_dist
3328 
3329  do n = 1, num_out
3330  ! print*, "n = ", n
3331  lon = lon_out(n)
3332  lat = lat_out(n)
3333  !--- find the j-index for the nearest point
3334  i_c = 0; j_c = 0
3335  do j = 1, jm_in-1
3336  if(lat .LE. geolat_in(1,j) .and.
3337  & lat .GE. geolat_in(1,j+1)) then
3338  j_c = j
3339  endif
3340  enddo
3341  if(lat > geolat_in(1,1)) j_c = 1
3342  if(lat < geolat_in(1,jm_in)) j_c = jm_in
3343  ! print*, "lat =", lat, geolat_in(1,jm_in), geolat_in(1,jm_in-1)
3344  ! The input is Gaussian grid.
3345  jstart = max(j_c-numnbr,1)
3346  jend = min(j_c+numnbr,jm_in)
3347  dist = max_dist
3348  iindx(n) = 0
3349  jindx(n) = 0
3350  ! print*, "jstart, jend =", jstart, jend
3351  do j = jstart, jend; do i = 1,im_in
3352  if(bitmap_in(i,j) ) then
3353  ! print*, "bitmap_in is true"
3354  d = spherical_distance(lon_out(n),lat_out(n),
3355  & geolon_in(i,j), geolat_in(i,j))
3356  if( dist > d ) then
3357  iindx(n) = i; jindx(n) = j
3358  dist = d
3359  endif
3360  endif
3361  enddo; enddo
3362  if(iindx(n) ==0) then
3363  print*, "lon,lat=", lon,lat
3364  print*, "jstart, jend=", jstart, jend, dist
3365  print*, "FATAL ERROR in get mismatch_index: "
3366  print*, "did not find nearest points."
3367  call errexit(4)
3368  endif
3369  enddo
3370 
3371  end subroutine get_mismatch_index
3372 
3386  subroutine interpolate_mismatch(im_in, jm_in, data_in,
3387  & num_out, data_out, iindx, jindx)
3388  integer, intent(in) :: im_in, jm_in, num_out
3389  real, intent(in) :: data_in(im_in,jm_in)
3390  real, intent(out):: data_out(num_out)
3391  integer, intent(in) :: iindx(num_out), jindx(num_out)
3392 
3393  do n = 1, num_out
3394  data_out(n) = data_in(iindx(n),jindx(n))
3395  enddo
3396 
3397  end subroutine interpolate_mismatch
3398 
3440  SUBROUTINE makeoa3(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
3441  1 ORO,SLM,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
3442  2 IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,
3443  3 IMI,JMI,OA_IN,OL_IN,
3444  4 slm_in,lon_in,lat_in)
3445 
3446 ! Required when using iplib v4.0 or higher.
3447 #ifdef IP_V4
3448  use ipolates_mod
3449 #endif
3450 
3451  implicit none
3452  real, parameter :: MISSING_VALUE = -9999.
3453  real, parameter :: d2r = 3.14159265358979/180.
3454  real, PARAMETER :: R2D=180./3.14159265358979
3455  integer im,jm,imn,jmn,imi,jmi
3456  real GLAT(JMN)
3457  INTEGER ZAVG(IMN,JMN)
3458  real SLM(IM,JM)
3459  real ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
3460  real OA4(IM,JM,4)
3461  integer IOA4(IM,JM,4)
3462  real OA_IN(IMI,JMI,4), OL_IN(IMI,JMI,4)
3463  real slm_in(IMI,JMI)
3464  real lon_in(IMI,JMI), lat_in(IMI,JMI)
3465  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
3466  real lon_t(IM,JM), lat_t(IM,JM)
3467  real XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
3468  real XNSUM3(IM,JM),XNSUM4(IM,JM)
3469  real VAR(IM,JM),OL(IM,JM,4)
3470  integer i,j,ilist(IMN),numx,i1,j1,ii1
3471  integer KWD
3472  real LONO(4),LATO(4),LONI,LATI
3473  real DELXN,HC,HEIGHT,T
3474  integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
3475  logical inside_a_polygon
3476  integer jst, jen
3477  integer int_opt, ipopt(20), kgds_input(200), kgds_output(200)
3478  integer count_land_output
3479  integer ij, ijmdl_output, iret, num_mismatch_land, num
3480  integer ibo(1), ibi(1)
3481  logical*1, allocatable :: bitmap_input(:,:)
3482  logical*1, allocatable :: bitmap_output(:,:)
3483  integer, allocatable :: ijsav_land_output(:)
3484  real, allocatable :: lats_land_output(:)
3485  real, allocatable :: lons_land_output(:)
3486  real, allocatable :: output_data_land(:,:)
3487  real, allocatable :: lons_mismatch_output(:)
3488  real, allocatable :: lats_mismatch_output(:)
3489  real, allocatable :: data_mismatch_output(:)
3490  integer, allocatable :: iindx(:), jindx(:)
3491 C
3492 C---- GLOBAL XLAT AND XLON ( DEGREE )
3493 C
3494  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
3495 C
3496  ijmdl_output = im*jm
3497 
3498  DO j=1,jmn
3499  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3500  ENDDO
3501  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
3502 C
3503 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
3504 C
3505 C
3506  DO j=1,jm
3507  DO i=1,im
3508  xnsum(i,j) = 0.0
3509  elvmax(i,j) = oro(i,j)
3510  zmax(i,j) = 0.0
3511 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
3512 C IN A GRID BOX
3513  xnsum1(i,j) = 0.0
3514  xnsum2(i,j) = 0.0
3515  xnsum3(i,j) = 0.0
3516  xnsum4(i,j) = 0.0
3517  oro1(i,j) = oro(i,j)
3518  elvmax(i,j) = zmax(i,j)
3519  ENDDO
3520  ENDDO
3521 
3522 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3523 ! --- to JM or to JM1
3524  DO j=1,jm
3525 ! print*, "J=", J
3526  DO i=1,im
3527  hc = 1116.2 - 0.878 * var(i,j)
3528  lono(1) = lon_c(i,j)
3529  lono(2) = lon_c(i+1,j)
3530  lono(3) = lon_c(i+1,j+1)
3531  lono(4) = lon_c(i,j+1)
3532  lato(1) = lat_c(i,j)
3533  lato(2) = lat_c(i+1,j)
3534  lato(3) = lat_c(i+1,j+1)
3535  lato(4) = lat_c(i,j+1)
3536  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3537  do j1 = jst, jen; do ii1 = 1, numx
3538  i1 = ilist(ii1)
3539  loni = i1*delxn
3540  lati = -90 + j1*delxn
3541  if(inside_a_polygon(loni*d2r,lati*d2r,4,
3542  & lono*d2r,lato*d2r))then
3543 
3544  height = float(zavg(i1,j1))
3545  IF(height.LT.-990.) height = 0.0
3546  IF ( height .gt. oro(i,j) ) then
3547  if ( height .gt. zmax(i,j) )zmax(i,j) = height
3548  ENDIF
3549  endif
3550  ENDDO ; ENDDO
3551  ENDDO
3552  ENDDO
3553 
3554 C
3555 ! --- this will make work1 array take on oro's values on return
3556 ! --- this will make work1 array take on oro's values on return
3557  DO j=1,jm
3558  DO i=1,im
3559 
3560  oro1(i,j) = oro(i,j)
3561  elvmax(i,j) = zmax(i,j)
3562  ENDDO
3563  ENDDO
3564 
3565  DO kwd = 1, 4
3566  DO j=1,jm
3567  DO i=1,im
3568  oa4(i,j,kwd) = 0.0
3569  ol(i,j,kwd) = 0.0
3570  ENDDO
3571  ENDDO
3572  ENDDO
3573 
3574  !--- use the nearest point to do remapping.
3575  int_opt = 2
3576  ipopt=0
3577  kgds_input = 0
3578  kgds_input(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
3579  kgds_input(2) = imi ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
3580  kgds_input(3) = jmi ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
3581  kgds_input(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
3582  kgds_input(5) = 0 ! OCT 14-16 - LON OF ORIGIN
3583  kgds_input(6) = 128 ! OCT 17 - RESOLUTION FLAG
3584  kgds_input(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
3585  kgds_input(8) = nint(-360000./imi) ! OCT 21-23 - LON OF EXTREME POINT
3586  kgds_input(9) = nint((360.0 / float(imi))*1000.0)
3587  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
3588  kgds_input(10) = jmi /2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
3589  kgds_input(12) = 255 ! OCT 29 - RESERVED
3590  kgds_input(20) = 255 ! OCT 5 - NOT USED, SET TO 255
3591 
3592 
3593  kgds_output = -1
3594 ! KGDS_OUTPUT(1) = IDRT ! OCT 6 - TYPE OF GRID (GAUSSIAN)
3595  kgds_output(2) = im ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
3596  kgds_output(3) = jm ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
3597  kgds_output(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
3598  kgds_output(5) = 0 ! OCT 14-16 - LON OF ORIGIN
3599  kgds_output(6) = 128 ! OCT 17 - RESOLUTION FLAG
3600  kgds_output(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
3601  kgds_output(8) = nint(-360000./im) ! OCT 21-23 - LON OF EXTREME POINT
3602  kgds_output(9) = nint((360.0 / float(im))*1000.0)
3603  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
3604  kgds_output(10) = jm /2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
3605  kgds_output(12) = 255 ! OCT 29 - RESERVED
3606  kgds_output(20) = 255 ! OCT 5 - NOT USED, SET TO 255
3607 
3608  count_land_output=0
3609  print*, "sum(slm) = ", sum(slm)
3610  do ij = 1, ijmdl_output
3611  j = (ij-1)/im + 1
3612  i = mod(ij-1,im) + 1
3613  if (slm(i,j) > 0.0) then
3614  count_land_output=count_land_output+1
3615  endif
3616  enddo
3617  allocate(bitmap_input(imi,jmi))
3618  bitmap_input=.false.
3619  print*, "number of land input=", sum(slm_in)
3620  where(slm_in > 0.0) bitmap_input=.true.
3621  print*, "count(bitmap_input)", count(bitmap_input)
3622 
3623  allocate(bitmap_output(count_land_output,1))
3624  allocate(output_data_land(count_land_output,1))
3625  allocate(ijsav_land_output(count_land_output))
3626  allocate(lats_land_output(count_land_output))
3627  allocate(lons_land_output(count_land_output))
3628 
3629 
3630 
3631  count_land_output=0
3632  do ij = 1, ijmdl_output
3633  j = (ij-1)/im + 1
3634  i = mod(ij-1,im) + 1
3635  if (slm(i,j) > 0.0) then
3636  count_land_output=count_land_output+1
3637  ijsav_land_output(count_land_output)=ij
3638  lats_land_output(count_land_output)=lat_t(i,j)
3639  lons_land_output(count_land_output)=lon_t(i,j)
3640  endif
3641  enddo
3642 
3643  oa4 = 0.0
3644  ol = 0.0
3645  ibi = 1
3646 
3647  do kwd=1,4
3648  bitmap_output = .false.
3649  output_data_land = 0.0
3650  call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3651  & (imi*jmi), count_land_output,
3652  & 1, ibi, bitmap_input, oa_in(:,:,kwd),
3653  & count_land_output, lats_land_output,
3654  & lons_land_output, ibo,
3655  & bitmap_output, output_data_land, iret)
3656  if (iret /= 0) then
3657  print*,'- FATAL ERROR IN IPOLATES ',iret
3658  call errexit(4)
3659  endif
3660 
3661  num_mismatch_land = 0
3662  do ij = 1, count_land_output
3663  if (bitmap_output(ij,1)) then
3664  j = (ijsav_land_output(ij)-1)/im + 1
3665  i = mod(ijsav_land_output(ij)-1,im) + 1
3666  oa4(i,j,kwd)=output_data_land(ij,1)
3667  else ! default value
3668  num_mismatch_land = num_mismatch_land + 1
3669  endif
3670  enddo
3671  print*, "num_mismatch_land = ", num_mismatch_land
3672 
3673  if(.not. allocated(data_mismatch_output)) then
3674  allocate(lons_mismatch_output(num_mismatch_land))
3675  allocate(lats_mismatch_output(num_mismatch_land))
3676  allocate(data_mismatch_output(num_mismatch_land))
3677  allocate(iindx(num_mismatch_land))
3678  allocate(jindx(num_mismatch_land))
3679 
3680  num = 0
3681  do ij = 1, count_land_output
3682  if (.not. bitmap_output(ij,1)) then
3683  num = num+1
3684  lons_mismatch_output(num) = lons_land_output(ij)
3685  lats_mismatch_output(num) = lats_land_output(ij)
3686  endif
3687  enddo
3688 
3689  ! For those land points that with bitmap_output=.false. use
3690  ! the nearest land points to interpolate.
3691  print*,"before get_mismatch_index", count(bitmap_input)
3692  call get_mismatch_index(imi,jmi,lon_in*d2r,
3693  & lat_in*d2r,bitmap_input,num_mismatch_land,
3694  & lons_mismatch_output*d2r,lats_mismatch_output*d2r,
3695  & iindx, jindx )
3696  endif
3697 
3698  data_mismatch_output = 0
3699  call interpolate_mismatch(imi,jmi,oa_in(:,:,kwd),
3700  & num_mismatch_land,data_mismatch_output,iindx,jindx)
3701 
3702  num = 0
3703  do ij = 1, count_land_output
3704  if (.not. bitmap_output(ij,1)) then
3705  num = num+1
3706  j = (ijsav_land_output(ij)-1)/im + 1
3707  i = mod(ijsav_land_output(ij)-1,im) + 1
3708  oa4(i,j,kwd) = data_mismatch_output(num)
3709  if(i==372 .and. j== 611) then
3710  print*, "ij=",ij, num, oa4(i,j,kwd),iindx(num),jindx(num)
3711  endif
3712  endif
3713  enddo
3714 
3715 
3716  enddo
3717 
3718  !OL
3719  do kwd=1,4
3720  bitmap_output = .false.
3721  output_data_land = 0.0
3722  call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3723  & (imi*jmi), count_land_output,
3724  & 1, ibi, bitmap_input, ol_in(:,:,kwd),
3725  & count_land_output, lats_land_output,
3726  & lons_land_output, ibo,
3727  & bitmap_output, output_data_land, iret)
3728  if (iret /= 0) then
3729  print*,'- FATAL ERROR IN IPOLATES ',iret
3730  call errexit(4)
3731  endif
3732 
3733  num_mismatch_land = 0
3734  do ij = 1, count_land_output
3735  if (bitmap_output(ij,1)) then
3736  j = (ijsav_land_output(ij)-1)/im + 1
3737  i = mod(ijsav_land_output(ij)-1,im) + 1
3738  ol(i,j,kwd)=output_data_land(ij,1)
3739  else ! default value
3740  num_mismatch_land = num_mismatch_land + 1
3741  endif
3742  enddo
3743  print*, "num_mismatch_land = ", num_mismatch_land
3744 
3745  data_mismatch_output = 0
3746  call interpolate_mismatch(imi,jmi,ol_in(:,:,kwd),
3747  & num_mismatch_land,data_mismatch_output,iindx,jindx)
3748 
3749  num = 0
3750  do ij = 1, count_land_output
3751  if (.not. bitmap_output(ij,1)) then
3752  num = num+1
3753  j = (ijsav_land_output(ij)-1)/im + 1
3754  i = mod(ijsav_land_output(ij)-1,im) + 1
3755  ol(i,j,kwd) = data_mismatch_output(num)
3756  if(i==372 .and. j== 611) then
3757  print*, "ij=",ij, num, ol(i,j,kwd),iindx(num),jindx(num)
3758  endif
3759  endif
3760  enddo
3761 
3762 
3763  enddo
3764 
3765  deallocate(lons_mismatch_output,lats_mismatch_output)
3766  deallocate(data_mismatch_output,iindx,jindx)
3767  deallocate(bitmap_output,output_data_land)
3768  deallocate(ijsav_land_output,lats_land_output)
3769  deallocate(lons_land_output)
3770 
3771  DO kwd=1,4
3772  DO j=1,jm
3773  DO i=1,im
3774  t = oa4(i,j,kwd)
3775  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3776  ENDDO
3777  ENDDO
3778  ENDDO
3779 C
3780  ns0 = 0
3781  ns1 = 0
3782  ns2 = 0
3783  ns3 = 0
3784  ns4 = 0
3785  ns5 = 0
3786  ns6 = 0
3787  DO kwd=1,4
3788  DO j=1,jm
3789  DO i=1,im
3790  t = abs( oa4(i,j,kwd) )
3791  IF(t .EQ. 0.) THEN
3792  ioa4(i,j,kwd) = 0
3793  ns0 = ns0 + 1
3794  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
3795  ioa4(i,j,kwd) = 1
3796  ns1 = ns1 + 1
3797  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
3798  ioa4(i,j,kwd) = 2
3799  ns2 = ns2 + 1
3800  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
3801  ioa4(i,j,kwd) = 3
3802  ns3 = ns3 + 1
3803  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
3804  ioa4(i,j,kwd) = 4
3805  ns4 = ns4 + 1
3806  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
3807  ioa4(i,j,kwd) = 5
3808  ns5 = ns5 + 1
3809  ELSE IF(t .GT. 10000.) THEN
3810  ioa4(i,j,kwd) = 6
3811  ns6 = ns6 + 1
3812  ENDIF
3813  ENDDO
3814  ENDDO
3815  ENDDO
3816 C
3817  WRITE(6,*) "! MAKEOA3 EXIT"
3818 C
3819  RETURN
3820  END SUBROUTINE makeoa3
3821 
3830  SUBROUTINE minmxj(IM,JM,A,title)
3831  implicit none
3832 
3833  real A(IM,JM),rmin,rmax
3834  integer i,j,IM,JM
3835  character*8 title
3836 
3837  rmin=1.e+10
3838  rmax=-rmin
3839 csela....................................................
3840 csela if(rmin.eq.1.e+10)return
3841 csela....................................................
3842  DO j=1,jm
3843  DO i=1,im
3844  if(a(i,j).ge.rmax)rmax=a(i,j)
3845  if(a(i,j).le.rmin)rmin=a(i,j)
3846  ENDDO
3847  ENDDO
3848  write(6,150)rmin,rmax,title
3849 150 format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
3850 C
3851  RETURN
3852  END
3853 
3865  SUBROUTINE mnmxja(IM,JM,A,imax,jmax,title)
3866  implicit none
3867 
3868  real A(IM,JM),rmin,rmax
3869  integer i,j,IM,JM,imax,jmax
3870  character*8 title
3871 
3872  rmin=1.e+10
3873  rmax=-rmin
3874 csela....................................................
3875 csela if(rmin.eq.1.e+10)return
3876 csela....................................................
3877  DO j=1,jm
3878  DO i=1,im
3879  if(a(i,j).ge.rmax)then
3880  rmax=a(i,j)
3881  imax=i
3882  jmax=j
3883  endif
3884  if(a(i,j).le.rmin)rmin=a(i,j)
3885  ENDDO
3886  ENDDO
3887  write(6,150)rmin,rmax,title
3888 150 format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
3889 C
3890  RETURN
3891  END
3892 
3897  subroutine read_g(glob)
3898  implicit none
3899 
3900  include 'netcdf.inc'
3901 
3902  integer*2, intent(out) :: glob(360*120,180*120)
3903 
3904  integer :: ncid, error, id_var, fsize
3905 
3906  fsize=65536
3907 
3908  error=nf__open("./topography.gmted2010.30s.nc",
3909  & nf_nowrite,fsize,ncid)
3910  call netcdf_err(error, 'Open file topography.gmted2010.30s.nc' )
3911  error=nf_inq_varid(ncid, 'topo', id_var)
3912  call netcdf_err(error, 'Inquire varid of topo')
3913  error=nf_get_var_int2(ncid, id_var, glob)
3914  call netcdf_err(error, 'Read topo')
3915  error = nf_close(ncid)
3916 
3917  print*,' '
3918  call maxmin (glob,360*120*180*120,'global0')
3919 
3920  return
3921  end subroutine read_g
3922 
3930  subroutine maxmin(ia,len,tile)
3931 ccmr
3932  implicit none
3933 ccmr
3934  integer*2 ia(len)
3935  character*7 tile
3936  integer iaamax, iaamin, len, m, ja, kount
3937  integer(8) sum2,std,mean,isum
3938  integer i_count_notset,kount_9
3939 ! --- missing is -9999
3940 c
3941  isum = 0
3942  sum2 = 0
3943  kount = 0
3944  kount_9 = 0
3945  iaamax = -9999999
3946 ccmr iaamin = 1
3947  iaamin = 9999999
3948  i_count_notset=0
3949  do 10 m=1,len
3950  ja=ia(m)
3951 ccmr if ( ja .lt. 0 ) print *,' ja < 0:',ja
3952 ccmr if ( ja .eq. -9999 ) goto 10
3953  if ( ja .eq. -9999 ) then
3954  kount_9=kount_9+1
3955  goto 10
3956  endif
3957  if ( ja .eq. -12345 ) i_count_notset=i_count_notset+1
3958 ccmr if ( ja .eq. 0 ) goto 11
3959  iaamax = max0( iaamax, ja )
3960  iaamin = min0( iaamin, ja )
3961 ! iaamax = max0( iaamax, ia(m,j) )
3962 ! iaamin = min0( iaamin, ia(m,j) )
3963  11 continue
3964  kount = kount + 1
3965  isum = isum + ja
3966 ccmr sum2 = sum2 + ifix( float(ja) * float(ja) )
3967  sum2 = sum2 + ja*ja
3968  10 continue
3969 !
3970  mean = isum/kount
3971  std = ifix(sqrt(float((sum2/(kount))-mean**2)))
3972  print*,tile,' max=',iaamax,' min=',iaamin,' sum=',isum,
3973  & ' i_count_notset=',i_count_notset
3974  print*,tile,' mean=',mean,' std.dev=',std,
3975  & ' ko9s=',kount,kount_9,kount+kount_9
3976  return
3977  end
3978 
3988  SUBROUTINE minmaxj(IM,JM,A,title)
3989  implicit none
3990 
3991  real(kind=4) A(IM,JM),rmin,rmax,undef
3992  integer i,j,IM,JM,imax,jmax,imin,jmin,iundef
3993  character*8 title,chara
3994  data chara/' '/
3995  chara=title
3996  rmin=1.e+10
3997  rmax=-rmin
3998  imax=0
3999  imin=0
4000  jmax=0
4001  jmin=0
4002  iundef=0
4003  undef=-9999.
4004 csela....................................................
4005 csela if(rmin.eq.1.e+10)return
4006 csela....................................................
4007  DO j=1,jm
4008  DO i=1,im
4009  if(a(i,j).ge.rmax)then
4010  rmax=a(i,j)
4011  imax=i
4012  jmax=j
4013  endif
4014  if(a(i,j).le.rmin)then
4015  if ( a(i,j) .eq. undef ) then
4016  iundef = iundef + 1
4017  else
4018  rmin=a(i,j)
4019  imin=i
4020  jmin=j
4021  endif
4022  endif
4023  ENDDO
4024  ENDDO
4025  write(6,150)chara,rmin,imin,jmin,rmax,imax,jmax,iundef
4026 150 format(1x,a8,2x,'rmin=',e13.4,2i6,2x,'rmax=',e13.4,3i6)
4027 C
4028  RETURN
4029  END
4030 
4040  subroutine latlon2xyz(siz,lon, lat, x, y, z)
4041  implicit none
4042  integer, intent(in) :: siz
4043  real, intent(in) :: lon(siz), lat(siz)
4044  real, intent(out) :: x(siz), y(siz), z(siz)
4045 
4046  integer n
4047 
4048  do n = 1, siz
4049  x(n) = cos(lat(n))*cos(lon(n))
4050  y(n) = cos(lat(n))*sin(lon(n))
4051  z(n) = sin(lat(n))
4052  enddo
4053  end
4054 
4062  FUNCTION spherical_angle(v1, v2, v3)
4063  implicit none
4064  real, parameter :: epsln30 = 1.e-30
4065  real, parameter :: pi=3.1415926535897931
4066  real v1(3), v2(3), v3(3)
4067  real spherical_angle
4068 
4069  real px, py, pz, qx, qy, qz, ddd;
4070 
4071  ! vector product between v1 and v2
4072  px = v1(2)*v2(3) - v1(3)*v2(2)
4073  py = v1(3)*v2(1) - v1(1)*v2(3)
4074  pz = v1(1)*v2(2) - v1(2)*v2(1)
4075  ! vector product between v1 and v3
4076  qx = v1(2)*v3(3) - v1(3)*v3(2);
4077  qy = v1(3)*v3(1) - v1(1)*v3(3);
4078  qz = v1(1)*v3(2) - v1(2)*v3(1);
4079 
4080  ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
4081  if ( ddd <= 0.0 ) then
4082  spherical_angle = 0.
4083  else
4084  ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
4085  if( abs(ddd-1) < epsln30 ) ddd = 1;
4086  if( abs(ddd+1) < epsln30 ) ddd = -1;
4087  if ( ddd>1. .or. ddd<-1. ) then
4088  !FIX to correctly handle co-linear points (angle near pi or 0) */
4089  if (ddd < 0.) then
4090  spherical_angle = pi
4091  else
4092  spherical_angle = 0.
4093  endif
4094  else
4095  spherical_angle = acos( ddd )
4096  endif
4097  endif
4098 
4099  return
4100  END
4101 
4112  FUNCTION inside_a_polygon(lon1, lat1, npts, lon2, lat2)
4113  implicit none
4114  real, parameter :: epsln10 = 1.e-10
4115  real, parameter :: epsln8 = 1.e-8
4116  real, parameter :: pi=3.1415926535897931
4117  real, parameter :: range_check_criteria=0.05
4118  real :: anglesum, angle, spherical_angle
4119  integer i, ip1
4120  real lon1, lat1
4121  integer npts
4122  real lon2(npts), lat2(npts)
4123  real x2(npts), y2(npts), z2(npts)
4124  real lon1_1d(1), lat1_1d(1)
4125  real x1(1), y1(1), z1(1)
4126  real pnt0(3),pnt1(3),pnt2(3)
4127  logical inside_a_polygon
4128  real max_x2,min_x2,max_y2,min_y2,max_z2,min_z2
4129  !first convert to cartesian grid */
4130  call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
4131  lon1_1d(1) = lon1
4132  lat1_1d(1) = lat1
4133  call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
4134  inside_a_polygon = .false.
4135  max_x2 = maxval(x2)
4136  if( x1(1) > max_x2+range_check_criteria ) return
4137  min_x2 = minval(x2)
4138  if( x1(1)+range_check_criteria < min_x2 ) return
4139  max_y2 = maxval(y2)
4140  if( y1(1) > max_y2+range_check_criteria ) return
4141  min_y2 = minval(y2)
4142  if( y1(1)+range_check_criteria < min_y2 ) return
4143  max_z2 = maxval(z2)
4144  if( z1(1) > max_z2+range_check_criteria ) return
4145  min_z2 = minval(z2)
4146  if( z1(1)+range_check_criteria < min_z2 ) return
4147 
4148  pnt0(1) = x1(1)
4149  pnt0(2) = y1(1)
4150  pnt0(3) = z1(1)
4151 
4152  anglesum = 0;
4153  do i = 1, npts
4154  if(abs(x1(1)-x2(i)) < epsln10 .and.
4155  & abs(y1(1)-y2(i)) < epsln10 .and.
4156  & abs(z1(1)-z2(i)) < epsln10 ) then ! same as the corner point
4157  inside_a_polygon = .true.
4158  return
4159  endif
4160  ip1 = i+1
4161  if(ip1>npts) ip1 = 1
4162  pnt1(1) = x2(i)
4163  pnt1(2) = y2(i)
4164  pnt1(3) = z2(i)
4165  pnt2(1) = x2(ip1)
4166  pnt2(2) = y2(ip1)
4167  pnt2(3) = z2(ip1)
4168 
4169  angle = spherical_angle(pnt0, pnt2, pnt1);
4170 ! anglesum = anglesum + spherical_angle(pnt0, pnt2, pnt1);
4171  anglesum = anglesum + angle
4172  enddo
4173 
4174  if(abs(anglesum-2*pi) < epsln8) then
4175  inside_a_polygon = .true.
4176  else
4177  inside_a_polygon = .false.
4178  endif
4179 
4180  return
4181 
4182  end
4183 
4207  function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
4208  & glat,zavg,zslm,delxn)
4209  implicit none
4210 
4211  real get_xnsum
4212  real, intent(in) :: lon1,lat1,lon2,lat2,delxn
4213  integer, intent(in) :: imn,jmn
4214  real, intent(in) :: glat(jmn)
4215  integer, intent(in) :: zavg(imn,jmn),zslm(imn,jmn)
4216  integer i, j, ist, ien, jst, jen, i1
4217  real oro, height
4218  real xland,xwatr,xl1,xs1,slm,xnsum
4219  !---figure out ist,ien,jst,jen
4220  do j = 1, jmn
4221  if( glat(j) .GT. lat1 ) then
4222  jst = j
4223  exit
4224  endif
4225  enddo
4226  do j = 1, jmn
4227  if( glat(j) .GT. lat2 ) then
4228  jen = j
4229  exit
4230  endif
4231  enddo
4232 
4233 
4234  ist = lon1/delxn + 1
4235  ien = lon2/delxn
4236  if(ist .le.0) ist = ist + imn
4237  if(ien < ist) ien = ien + imn
4238 
4239  !--- compute average oro
4240  oro = 0.0
4241  xnsum = 0
4242  xland = 0
4243  xwatr = 0
4244  xl1 = 0
4245  xs1 = 0
4246  do j = jst,jen
4247  do i1 = 1, ien - ist + 1
4248  i = ist + i1 -1
4249  if( i .LE. 0) i = i + imn
4250  if( i .GT. imn) i = i - imn
4251  xland = xland + float(zslm(i,j))
4252  xwatr = xwatr + float(1-zslm(i,j))
4253  xnsum = xnsum + 1.
4254  height = float(zavg(i,j))
4255  IF(height.LT.-990.) height = 0.0
4256  xl1 = xl1 + height * float(zslm(i,j))
4257  xs1 = xs1 + height * float(1-zslm(i,j))
4258  enddo
4259  enddo
4260  if( xnsum > 1.) THEN
4261  slm = float(nint(xland/xnsum))
4262  IF(slm.NE.0.) THEN
4263  oro= xl1 / xland
4264  ELSE
4265  oro = xs1 / xwatr
4266  ENDIF
4267  ENDIF
4268 
4269  get_xnsum = 0
4270  do j = jst, jen
4271  do i1= 1, ien-ist+1
4272  i = ist + i1 -1
4273  if( i .LE. 0) i = i + imn
4274  if( i .GT. imn) i = i - imn
4275  height = float(zavg(i,j))
4276  IF(height.LT.-990.) height = 0.0
4277  IF ( height .gt. oro ) get_xnsum = get_xnsum + 1
4278  enddo
4279  enddo
4280 
4281  end function get_xnsum
4282 
4310  subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
4311  & glat,zavg,delxn,xnsum1,xnsum2,HC)
4312  implicit none
4313 
4314  real, intent(out) :: xnsum1,xnsum2,HC
4315  real lon1,lat1,lon2,lat2,delxn
4316  integer IMN,JMN
4317  real glat(JMN)
4318  integer zavg(IMN,JMN)
4319  integer i, j, ist, ien, jst, jen, i1
4320  real HEIGHT, var
4321  real XW1,XW2,xnsum
4322  !---figure out ist,ien,jst,jen
4323  do j = 1, jmn
4324  if( glat(j) .GT. lat1 ) then
4325  jst = j
4326  exit
4327  endif
4328  enddo
4329  do j = 1, jmn
4330  if( glat(j) .GT. lat2 ) then
4331  jen = j
4332  exit
4333  endif
4334  enddo
4335 
4336 
4337  ist = lon1/delxn + 1
4338  ien = lon2/delxn
4339  if(ist .le.0) ist = ist + imn
4340  if(ien < ist) ien = ien + imn
4341 
4342  !--- compute average oro
4343  xnsum = 0
4344  xw1 = 0
4345  xw2 = 0
4346  do j = jst,jen
4347  do i1 = 1, ien - ist + 1
4348  i = ist + i1 -1
4349  if( i .LE. 0) i = i + imn
4350  if( i .GT. imn) i = i - imn
4351  xnsum = xnsum + 1.
4352  height = float(zavg(i,j))
4353  IF(height.LT.-990.) height = 0.0
4354  xw1 = xw1 + height
4355  xw2 = xw2 + height ** 2
4356  enddo
4357  enddo
4358  var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
4359  hc = 1116.2 - 0.878 * var
4360  xnsum1 = 0
4361  xnsum2 = 0
4362  do j = jst, jen
4363  do i1= 1, ien-ist+1
4364  i = ist + i1 -1
4365  if( i .LE. 0) i = i + imn
4366  if( i .GT. imn) i = i - imn
4367  height = float(zavg(i,j))
4368  IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4369  xnsum2 = xnsum2 + 1
4370  enddo
4371  enddo
4372 
4373  end subroutine get_xnsum2
4374 
4403  subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
4404  & glat,zavg,delxn,xnsum1,xnsum2,HC)
4405  implicit none
4406 
4407  real, intent(out) :: xnsum1,xnsum2
4408  real lon1,lat1,lon2,lat2,delxn
4409  integer IMN,JMN
4410  real glat(JMN)
4411  integer zavg(IMN,JMN)
4412  integer i, j, ist, ien, jst, jen, i1
4413  real HEIGHT, HC
4414  !---figure out ist,ien,jst,jen
4415  ! if lat1 or lat 2 is 90 degree. set jst = JMN
4416  jst = jmn
4417  jen = jmn
4418  do j = 1, jmn
4419  if( glat(j) .GT. lat1 ) then
4420  jst = j
4421  exit
4422  endif
4423  enddo
4424  do j = 1, jmn
4425  if( glat(j) .GT. lat2 ) then
4426  jen = j
4427  exit
4428  endif
4429  enddo
4430 
4431 
4432  ist = lon1/delxn + 1
4433  ien = lon2/delxn
4434  if(ist .le.0) ist = ist + imn
4435  if(ien < ist) ien = ien + imn
4436 
4437  xnsum1 = 0
4438  xnsum2 = 0
4439  do j = jst, jen
4440  do i1= 1, ien-ist+1
4441  i = ist + i1 -1
4442  if( i .LE. 0) i = i + imn
4443  if( i .GT. imn) i = i - imn
4444  height = float(zavg(i,j))
4445  IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4446  xnsum2 = xnsum2 + 1
4447  enddo
4448  enddo
4449 
4450  end subroutine get_xnsum3
4455  real function timef()
4456  character(8) :: date
4457  character(10) :: time
4458  character(5) :: zone
4459  integer,dimension(8) :: values
4460  integer :: total
4461  real :: elapsed
4462  call date_and_time(date,time,zone,values)
4463  total=(3600*values(5))+(60*values(6))
4464  * +values(7)
4465  elapsed=float(total) + (1.0e-3*float(values(8)))
4466  timef=elapsed
4467  return
4468  end
write_mask_netcdf
subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geolat)
Write the land mask file.
Definition: netcdf_io.F90:245
read_mask
subroutine read_mask(merge_file, slm, land_frac, lake_frac, im, jm)
Read the land mask file.
Definition: netcdf_io.F90:336
makemt
subroutine makemt(ZAVG, ZSLM, ORO, SLM, VAR, VAR4, GLAT, IST, IEN, JST, JEN, IM, JM, IMN, JMN, XLAT, numi)
Create the orography, land-mask, standard deviation of orography and the convexity on a model gaussia...
Definition: mtnlm7_oclsm.F:1321
get_index
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...
Definition: mtnlm7_oclsm.F:1482
timef
real function timef()
Get the date/time for the system clock.
Definition: mtnlm7_oclsm.F:4456
netcdf_err
subroutine netcdf_err(err, string)
Check NetCDF error code and output the error message.
Definition: netcdf_io.F90:219
get_lat_angle
real function get_lat_angle(dy, DEGRAD)
Convert the &#39;y&#39; direction distance of a cubed-sphere grid point to the corresponding distance in lati...
Definition: mtnlm7_oclsm.F:2830
mnmxja
subroutine mnmxja(IM, JM, A, imax, jmax, title)
Print out the maximum and minimum values of an array.
Definition: mtnlm7_oclsm.F:3866
makepc2
subroutine makepc2(ZAVG, ZSLM, THETA, GAMMA, SIGMA, GLAT, IM, JM, IMN, JMN, lon_c, lat_c, SLM)
Make the principle coordinates - slope of orography, anisotropy, angle of mountain range with respect...
Definition: mtnlm7_oclsm.F:2216
tersub
subroutine tersub(IMN, JMN, IM, JM, NM, NR, NF0, NF1, NW, EFAC, OUTGRID, INPUTOROG, MASK_ONLY, MERGE_FILE)
Driver routine to compute terrain.
Definition: mtnlm7_oclsm.F:167
spherical_angle
real function spherical_angle(v1, v2, v3)
Compute spherical angle.
Definition: mtnlm7_oclsm.F:4063
minmaxj
subroutine minmaxj(IM, JM, A, title)
Print out the maximum and minimum values of an array and their i/j location.
Definition: mtnlm7_oclsm.F:3989
lake_frac
program lake_frac
This program computes lake fraction and depth numbers for FV3 cubed sphere grid cells, from a high resolution lat/lon data set.
Definition: lakefrac.F90:21
write_netcdf
subroutine write_netcdf(im, jm, slm, land_frac, oro, orf, hprime, ntiles, tile, geolon, geolat, lon, lat)
Write out orography file in netcdf format.
Definition: netcdf_io.F90:22
makepc
subroutine makepc(ZAVG, ZSLM, THETA, GAMMA, SIGMA, GLAT, IST, IEN, JST, JEN, IM, JM, IMN, JMN, XLAT, numi)
Make the principle coordinates - slope of orography, anisotropy, angle of mountain range with respect...
Definition: mtnlm7_oclsm.F:1941
get_xnsum3
subroutine get_xnsum3(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
Definition: mtnlm7_oclsm.F:4405
inside_a_polygon
logical function inside_a_polygon(lon1, lat1, npts, lon2, lat2)
Check if a point is inside a polygon.
Definition: mtnlm7_oclsm.F:4113
spherical_distance
real function spherical_distance(theta1, phi1, theta2, phi2)
Compute a great circle distance between two points.
Definition: mtnlm7_oclsm.F:3270
makeoa3
subroutine makeoa3(ZAVG, VAR, GLAT, OA4, OL, IOA4, ELVMAX, ORO, SLM, oro1, XNSUM, XNSUM1, XNSUM2, XNSUM3, XNSUM4, IM, JM, IMN, JMN, lon_c, lat_c, lon_t, lat_t, IMI, JMI, OA_IN, OL_IN, slm_in, lon_in, lat_in)
Create orographic asymmetry and orographic length scale on the model grid.
Definition: mtnlm7_oclsm.F:3445
get_mismatch_index
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...
Definition: mtnlm7_oclsm.F:3306
make_mask
subroutine make_mask(zslm, SLM, land_frac, GLAT, IM, JM, IMN, JMN, lon_c, lat_c)
Create the land-mask, land fraction.
Definition: mtnlm7_oclsm.F:1571
maxmin
subroutine maxmin(ia, len, tile)
Print the maximum, mininum, mean and standard deviation of an array.
Definition: mtnlm7_oclsm.F:3931
makemt2
subroutine makemt2(ZAVG, ZSLM, ORO, SLM, VAR, VAR4, GLAT, IM, JM, IMN, JMN, lon_c, lat_c, lake_frac, land_frac)
Create the orography, land-mask, land fraction, standard deviation of orography and the convexity on ...
Definition: mtnlm7_oclsm.F:1706
interpolate_mismatch
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. ...
Definition: mtnlm7_oclsm.F:3388
get_lon_angle
real function get_lon_angle(dx, lat, DEGRAD)
Convert the &#39;x&#39; direction distance of a cubed-sphere grid point to the corresponding distance in long...
Definition: mtnlm7_oclsm.F:2811
get_xnsum2
subroutine get_xnsum2(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
Definition: mtnlm7_oclsm.F:4312
makeoa2
subroutine makeoa2(ZAVG, zslm, VAR, GLAT, OA4, OL, IOA4, ELVMAX, ORO, oro1, XNSUM, XNSUM1, XNSUM2, XNSUM3, XNSUM4, IM, JM, IMN, JMN, lon_c, lat_c, lon_t, lat_t, dx, dy, is_south_pole, is_north_pole)
Create orographic asymmetry and orographic length scale on the model grid.
Definition: mtnlm7_oclsm.F:2880
latlon2xyz
subroutine latlon2xyz(siz, lon, lat, x, y, z)
Convert from latitude and longitude to x,y,z coordinates.
Definition: mtnlm7_oclsm.F:4041
minmxj
subroutine minmxj(IM, JM, A, title)
Print out the maximum and minimum values of an array.
Definition: mtnlm7_oclsm.F:3831
makeoa
subroutine makeoa(ZAVG, VAR, GLAT, OA4, OL, IOA4, ELVMAX, ORO, oro1, XNSUM, XNSUM1, XNSUM2, XNSUM3, XNSUM4, IST, IEN, JST, JEN, IM, JM, IMN, JMN, XLAT, numi)
Create orographic asymmetry and orographic length scale on the model grid.
Definition: mtnlm7_oclsm.F:2485
get_xnsum
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 ...
Definition: mtnlm7_oclsm.F:4209
read_g
subroutine read_g(glob)
Read input global 30-arc second orography data.
Definition: mtnlm7_oclsm.F:3898