orog_mask_tools  1.3.0
 All Data Structures Files Functions Variables
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 navy
6 C> 10-minute 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> - UNIT235 - GTOPO 30" AVR for ZAVG elevation
53 C> - UNIT10 - 30" UMD land (lake) cover mask see MSKSRC switch
54 C> - XUNIT11 - GTOPO AVR
55 C> - XUNIT12 - GTOPO STD DEV
56 C> - XUNIT13 - GTOPO MAX
57 C> - UNIT14 - GTOPO SLM (10' NAVY if switched to get lakes
58 C> - UNIT15 - GICE Grumbine 30" RAMP Antarctica orog IMNx3616
59 C> - UNIT25 - Ocean land-sea mask on gaussian grid
60 C>
61 C> OUTPUT FILES:
62 C> - UNIT51 - SEA-LAND MASK (IM,JM)
63 C> - UNIT52 - GRIDDED OROGRAPHY (IM,JM)
64 C> - UNIT54 - SPECTRAL OROGRAPHY ((NM+1)*((NR+1)*NM+2))
65 C> - UNIT55 - UNFILTERED GRIDDED OROGRAPHY (IM,JM)
66 C> - UNIT57 - GRIB GRIDDED OROGRAPHY (IM,JM)
67 C>
68 C> SUBPROGRAMS CALLED:
69 C> - UNIQUE:
70 C> - TERSUB - MAIN SUBPROGRAM
71 C> - SPLAT - COMPUTE GAUSSIAN LATITUDES OR EQUALLY-SPACED LATITUDES
72 C> - LIBRARY:
73 C> - SPTEZ - SPHERICAL TRANSFORM
74 C> - GBYTES - UNPACK BITS
75 C>
76 C> @return 0 for success, error code otherwise.
77  include 'netcdf.inc'
78  logical fexist, opened
79  integer fsize, ncid, error, id_dim, nx, ny
80  character(len=256) :: outgrid = "none"
81  character(len=256) :: inputorog = "none"
82  logical :: do_oa = .true. ! create oa and ol data.
83  logical :: grid_from_file = .true.
84  integer :: mtnres,im,jm,nm,nr,nf0,nf1,efac,blat,nw
85  fsize=65536
86  READ(5,*) mtnres,im,jm,nm,nr,nf0,nf1,efac,blat
87  READ(5,*) outgrid
88  READ(5,*) inputorog
89 ! MTNRES=1
90 ! IM=48
91 ! JM=48
92 ! NM=46
93 ! NF0=0
94 ! NF1=0
95 ! efac=0
96 ! blat=0
97 ! NR=0
98 ! OUTGRID = "C48_grid.tile1.nc"
99 ! INPUTOROG = "oro.288x144.nc"
100  print*, "INPUTOROG=", trim(inputorog)
101  print*, "IM,JM=", im, jm
102 ! --- MTNRES defines the input (highest) elev resolution
103 ! --- =1 is topo30 30" in units of 1/2 minute.
104 ! so MTNRES for old values must be *2.
105 ! =16 is now Song Yu's 8' orog the old ops standard
106 ! --- other possibilities are =8 for 4' and =4 for 2' see
107 ! HJ for T1000 test. Must set to 1 for now.
108  mtnres=1
109  print*, mtnres,im,jm,nm,nr,nf0,nf1,efac,blat
110  nw=(nm+1)*((nr+1)*nm+2)
111  imn = 360*120/mtnres
112  jmn = 180*120/mtnres
113  print *, ' Starting terr12 mtnlm7_slm30.f IMN,JMN:',imn,jmn
114 
115 ! --- read the grid resolution if the OUTGRID exists.
116  if( trim(outgrid) .NE. "none" ) then
117  inquire(file=trim(outgrid), exist=fexist)
118  if(.not. fexist) then
119  print*, "file "//trim(outgrid)//" does not exist"
120  CALL errexit(4)
121  endif
122  do ncid = 103, 512
123  inquire( ncid,opened=opened )
124  if( .NOT.opened )exit
125  end do
126 
127  print*, "outgrid=", trim(outgrid)
128  error=nf__open(trim(outgrid),nf_nowrite,fsize,ncid)
129  call netcdf_err(error, 'Open file '//trim(outgrid) )
130  error=nf_inq_dimid(ncid, 'nx', id_dim)
131  call netcdf_err(error, 'inquire dimension nx from file '//
132  & trim(outgrid) )
133  error=nf_inq_dimlen(ncid,id_dim,nx)
134  call netcdf_err(error, 'inquire dimension nx length '//
135  & 'from file '//trim(outgrid) )
136 
137  error=nf_inq_dimid(ncid, 'ny', id_dim)
138  call netcdf_err(error, 'inquire dimension ny from file '//
139  & trim(outgrid) )
140  error=nf_inq_dimlen(ncid,id_dim,ny)
141  call netcdf_err(error, 'inquire dimension ny length '//
142  & 'from file '//trim(outgrid) )
143  print*, "nx = ", nx
144  if(im .ne. nx/2) then
145  print*, "IM=",im, " /= grid file nx/2=",nx/2
146  print*, "Set IM = ", nx/2
147  im = nx/2
148  endif
149  if(jm .ne. ny/2) then
150  print*, "JM=",jm, " /= grid file ny/2=",ny/2
151  print*, "Set JM = ", ny/2
152  jm = ny/2
153  endif
154  error=nf_close(ncid)
155  call netcdf_err(error, 'close file '//trim(outgrid) )
156 
157  endif
158 
159 
160  CALL tersub(imn,jmn,im,jm,nm,nr,nf0,nf1,nw,efac,blat,
161  & outgrid,inputorog)
162  stop
163  END
164 
185  SUBROUTINE tersub(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
186  & outgrid,inputorog)
187 !jaa use ipfort
188  implicit none
189  include 'netcdf.inc'
190 C
191  integer :: imn,jmn,im,jm,nw
192  character(len=*), intent(in) :: outgrid
193  character(len=*), intent(in) :: inputorog
194  integer :: nr,nf0,nf1
195  real, parameter :: missing_value=-9999.
196  real, PARAMETER :: pi=3.1415926535897931
197  integer :: efac, blat
198  integer, PARAMETER :: nmt=14
199  INTEGER zslmx(2700,1350)
200  integer nm
201  logical latlongrid
202  INTEGER,allocatable:: zavg(:,:),zslm(:,:)
203  REAL(4),allocatable:: gice(:,:),oclsm(:,:)
204  real :: degrad
205  integer*1,allocatable:: umd(:,:)
206  integer*1 i3save
207  integer*2 glob(imn,jmn), i2save
208  logical grid_from_file
209  INTEGER kpds(200),kgds(200), zsave1,zsave2,itopo,kount
210  INTEGER kount2,islmx,jslmx,oldslm,msksrc,mskocn,notocn
211  REAL cosclt(jm),wgtclt(jm),rclt(jm),xlat(jm),diffx(jm/2)
212  REAL xlon(im)
213  LOGICAL is_south_pole(im,jm), is_north_pole(im,jm)
214  REAL geolon(im,jm),geolat(im,jm)
215  REAL, allocatable :: tmpvar(:,:)
216  REAL geolon_c(im+1,jm+1),geolat_c(im+1,jm+1)
217  REAL dx(im,jm),dy(im,jm)
218  REAL slm(im,jm),oro(im,jm),var(im,jm),ors(nw),orf(im,jm)
219  REAL land_frac(im,jm)
220  REAL theta(im,jm),gamma(im,jm),sigma(im,jm),elvmax(im,jm)
221  REAL wz4(im,jm),var4(im,jm),oa(im,jm,4),ol(im,jm,4),slmi(im,jm)
222  integer ist(im,jm),ien(im,jm),jst(jm),jen(jm)
223  integer iwork(im,jm,4)
224  real work1(im,jm),work2(im,jm),work3(im,jm),work4(im,jm)
225  real work5(im,jm),work6(im,jm),glat(jmn)
226  LOGICAL spectr, revlat, filter
227  logical fexist
228  real hprime(im,jm,14)
229  real oaa(4),ola(4),sumdif,avedif,alon,alat
230  real, allocatable :: oa_in(:,:,:), ol_in(:,:,:)
231  real, allocatable :: slm_in(:,:), lon_in(:,:), lat_in(:,:)
232  integer numi(jm),ios,iosg,latg2,istat
233  integer maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
234  integer lonsperlat(jm/2),itest,jtest
235  integer i, j, nx, ny, ncid, js, jn, iw, ie, k
236  integer it,jt,i1,error,id_dim,id_var,nx_in,ny_in
237  integer i_south_pole,j_south_pole,i_north_pole,j_north_pole
238  real maxlat, minlat
239  logical opened
240  logical lb(im*jm)
241  integer fsize,wgta,in,inw,ine,is,isw,ise,m,n,imt,iret
242  complex ffj(im/2+1)
243  real dlat,phi,delxn,rs,rn,slma,oroa,vara,var4a,xn,xs,fff
244  real www
245  real :: timef,tbeg,tend,tbeg1
246  logical :: output_binary
247  output_binary = .false.
248  tbeg1=timef()
249  tbeg=timef()
250  fsize = 65536
251 
252  allocate (zavg(imn,jmn))
253  allocate (zslm(imn,jmn))
254  allocate (gice(imn+1,3601))
255  allocate (umd(imn,jmn))
256  allocate (oclsm(im,jm))
257 
258 !
259 ! SET CONSTANTS AND ZERO FIELDS
260 !
261  degrad = 180./pi
262  spectr = nm .GT. 0 ! if NM <=0 grid is assumed lat/lon
263  filter = .true. ! Spectr Filter defaults true and set by NF1 & NF0
264 ! MSKSRC = 0 ! MSKSRC=0 navy 10 lake msk, 1 UMD 30, -1 no lakes
265  msksrc = 1 ! MSKSRC=0 navy 10 lake msk, 1 UMD 30, -1 no lakes
266  revlat = blat .LT. 0 ! Reverse latitude/longitude for output
267  itopo = 1 ! topo 30" read, otherwise tiles (opt offline)
268  mskocn = 1 ! Ocean land sea mask =1, =0 if not present
269  notocn = 1 ! =1 Ocean lsm input reverse: Ocean=1, land=0
270 ! --- The LSM Gaussian file from the ocean model sometimes arrives with
271 ! --- 0=Ocean and 1=Land or it arrives with 1=Ocean and 0=land without
272 ! --- metadata to distinguish its disposition. The AI below mitigates this.
273 
274  print *,' In TERSUB, ITOPO=',itopo
275  if (mskocn .eq. 1)then
276  print *,' Ocean Model LSM Present and '
277  print *, ' Overrides OCEAN POINTS in LSM: mskocn=',mskocn
278  if (notocn .eq. 1)then
279  print *,' Ocean LSM Reversed: NOTOCN=',notocn
280  endif
281  endif
282 C
283 C --- old S-Y. files
284 C- OPEN(UNIT=11,FORM='FORMATTED',ERR=900) ! average
285 C- OPEN(UNIT=12,FORM='FORMATTED',ERR=900) ! Std Dev
286 C- OPEN(UNIT=13,FORM='FORMATTED',ERR=900) ! maximum
287 C- OPEN(UNIT=14,FORM='FORMATTED',ERR=900) ! sea-land-lake-mask
288 C
289 ! --- READ(11,11) ZAVG
290 ! --- READ(12,11) ZVAR
291 ! --- READ(13,11) ZMAX
292 ! --- 11 FORMAT(20I4)
293 !
294 ! --- MSKSRC 0 navy 10' lake mask, =1 for 30" UMD lake mask,
295 ! --- MSKSRC internally set if above fails at -1 for no lakes
296 ! ---
297  IF (msksrc .eq. 0 ) then
298  READ(14,12,iostat=ios) zslmx
299  12 FORMAT(80i1)
300  if (ios.ne.0) then
301  msksrc=-1
302  print *,' navy10 lake mask rd fail -- ios,MSKSRC:',ios,msksrc
303  endif
304  ELSE
305  print *,' Attempt to open/read UMD 30" slmsk MSKSRC=',msksrc
306 ! --- not 0 so MSKSRC=1 and attempt to open/read UMD 30" slmsk
307 ! open(10,file=
308 ! &"/scratch2/portfolios/NCEPDEV/global/noscrub/Jordan.Alpert/wx23ja
309 ! &/terrain30/landcover30.fixed",
310 ! & recl=43200*21600, access='direct',iostat=istat)
311  open(10,file="landcover30.fixed",
312  & recl=43200*21600, access='direct',iostat=istat)
313 
314  IF (istat.ne.0) then
315  msksrc=-1
316  print *,' UMD lake mask open failed -- ios,MSKSRC:',istat,msksrc
317  ELSE
318 !
319  read(10, rec=1,iostat=istat) umd
320  print *,' UMD lake mask opened OK -- ios,MSKSRC:',istat,msksrc
321 !
322  ENDIF
323 ! --------------
324  IF (istat.ne.0) then
325 ! --- When UMD read fails attempt to read navy 10'
326  print *,' UMD lake mask rd err -- trying navy 10',istat
327  msksrc=0
328  print *,' ***** MSKSRC set to 0 MSKSRC=',msksrc
329  if (msksrc .eq. 0 ) then
330  READ(14,12,iostat=ios) zslmx
331  if (ios.ne.0) then
332  msksrc=-1
333  print *,' navy10 lake mask rd fail - ios,MSKSRC:',ios,msksrc
334  endif
335  endif
336  ELSE
337  print *,' UMD lake, UMD(50,50)=',umd(50,50),msksrc
338  ENDIF
339 ! --------------
340 ! --- good UMD land cover read and MSKSRC=1
341  ENDIF
342 C
343 C- READ_G for global 30" terrain
344 C
345  print *,' About to call read_g, ITOPO=',itopo
346  if ( itopo .ne. 0 ) call read_g(glob,itopo)
347 ! --- transpose even though glob 30" is from S to N and NCEP std is N to S
348  do j=1,jmn/2
349  do i=1,imn
350  jt=jmn - j + 1
351  i2save = glob(i,j)
352  glob(i,j)=glob(i,jt)
353  glob(i,jt) = i2save
354  enddo
355  enddo
356 ! --- transpose glob as USGS 30" is from dateline and NCEP std is 0
357  do j=1,jmn
358  do i=1,imn/2
359  it=imn/2 + i
360  i2save = glob(i,j)
361  glob(i,j)=glob(it,j)
362  glob(it,j) = i2save
363  enddo
364  enddo
365  print *,' After read_g, glob(500,500)=',glob(500,500)
366 !
367 
368 ! --- IMN,JMN
369  print*, ' IM, JM, NM, NR, NF0, NF1, EFAC, BLAT'
370  print*, im,jm,nm,nr,nf0,nf1,efac,blat
371  print *,' imn,jmn,glob(imn,jmn)=',imn,jmn,glob(imn,jmn)
372  print *,' UBOUND ZAVG=',ubound(zavg)
373  print *,' UBOUND glob=',ubound(glob)
374  print *,' UBOUND ZSLM=',ubound(zslm)
375  print *,' UBOUND GICE=',ubound(gice)
376  print *,' UBOUND OCLSM=',ubound(oclsm)
377 !
378 ! --- 0 is ocean and 1 is land for slm
379 !
380 C
381 ! --- ZSLM initialize with all land 1, ocean 0
382 ! ZSLM=1
383  do j=1,jmn
384  do i=1,imn
385  zslm(i,j)=1
386  enddo
387  enddo
388 
389  selectcase(msksrc)
390 C---- 30" sea land mask. 0 are water (lake or ocean)
391  CASE(1)
392 ! --- transpose even though glob 30" is from S to N and NCEP std is N to S
393  do j=1,jmn/2
394  do i=1,imn
395  jt=jmn - j + 1
396  i3save = umd(i,j)
397  umd(i,j)=umd(i,jt)
398  umd(i,jt) = i3save
399  enddo
400  enddo
401 ! --- transpose UMD as USGS 30" is from dateline and NCEP std is 0
402  do j=1,jmn
403  do i=1,imn/2
404  it=imn/2 + i
405  i3save = umd(i,j)
406  umd(i,j)=umd(it,j)
407  umd(it,j) = i3save
408  enddo
409  enddo
410 ! --- UMD slmsk with 30" lakes and set ZAVG from glob
411  do j=1,jmn
412  do i=1,imn
413  if ( umd(i,j) .eq. 0 ) zslm(i,j) = 0
414  zavg(i,j) = glob(i,j)
415  enddo
416  enddo
417 ! --- Global land in slm plus lakes on 30" grid and elev set over globe
418 ! ---
419 ! --- When navy 10' mask is set MSKSRC=0
420  CASE(0)
421 ! --- MSKSRC 0 navy 10' lake mask, =1 for 30" UMD lake mask, -1 no lakes
422  print *,' NAVY 10 (8) slmsk for lakes, MSKSRC=',msksrc
423  kount = 0
424  kount2 = 0
425  do j=1,jmn
426  oldslm = zslm(imn,j)
427  do i=1,imn
428  i1 = i + 1
429 ! --- slmsk with 10' lakes and set ZAVG from 30" glob
430  zavg(i,j) = glob(i,j)
431  if ( glob(i,j) .eq. -9999 ) then
432  zslm(i,j) = 0
433  kount = kount + 1
434  endif
435  islmx=(i-1)/16 + 1
436  jslmx=(j-1)/16 + 1
437  if ( zslmx(islmx,jslmx) .eq. 0 ) then
438  if ( j .gt. 8 .and. j .lt. jmn-8 ) then
439  if (i1 .gt. imn ) i1 = i1 - imn
440 ! -----
441  if(zslm(i,j).eq.1 .and. oldslm .eq. 1 .and. zslm(i1,j).eq.1)then
442  if (i .ne. 1) oldslm = zslm(i,j)
443  zslm(i,j) = 0
444  kount2 = kount2 + 1
445  endif
446 ! -----
447  endif
448  endif
449  enddo
450  enddo
451 ! ---
452  CASE(-1)
453  print *,' ***** set ZAVG and slm from 30" glob, MSKSRC=',msksrc
454  kount = 0
455  kount2 = 0
456  do j=1,jmn
457  do i=1,imn
458  i1 = i + 1
459 ! --- UMD slmsk with 10' lakes and set ZAVG from 30" glob
460  zavg(i,j) = glob(i,j)
461  if ( glob(i,j) .eq. -9999 ) then
462  zslm(i,j) = 0
463  kount = kount + 1
464  endif
465  enddo
466  enddo
467  END SELECT
468 ! ---
469 ! --- Fixing an error in the topo 30" data set at pole (-9999).
470  do i=1,imn
471  zslm(i,1)=0
472  zslm(i,jmn)=1
473  enddo
474 !
475 ! print *,' kount1,2,ZAVG(1,1),ZAVG(imn,jmn),ZAVG(500,500)',
476 ! & kount,kount2,ZAVG(1,1),ZAVG(imn,jmn),ZAVG(500,500)
477 ! --- The center of pixel (1,1) is 89.9958333N/179.9958333W with dx/dy
478 ! --- spacing of 1/120 degrees.
479 !
480 ! READ REDUCED GRID EXTENTS IF GIVEN
481 !
482  read(20,*,iostat=ios) latg2,lonsperlat
483  if(ios.ne.0.or.2*latg2.ne.jm) then
484  do j=1,jm
485  numi(j)=im
486  enddo
487  print *,ios,latg2,'COMPUTE TERRAIN ON A FULL GAUSSIAN GRID'
488  else
489  do j=1,jm/2
490  numi(j)=lonsperlat(j)
491  enddo
492  do j=jm/2+1,jm
493  numi(j)=lonsperlat(jm+1-j)
494  enddo
495  print *,ios,latg2,'COMPUTE TERRAIN ON A REDUCED GAUSSIAN GRID',
496  & numi
497 C print *,ios,latg2,'COMPUTE TERRAIN ON A REDUCED GAUSSIAN GRID'
498  endif
499 ! print *,ios,latg2,'TERRAIN ON GAUSSIAN GRID',numi
500 
501 !
502 ! This code assumes that lat runs from north to south for gg!
503 !
504  print *,' SPECTR=',spectr,' REVLAT=',revlat,' ** with GICE-07 **'
505  IF (spectr) THEN
506  CALL splat(4,jm,cosclt,wgtclt)
507  DO j=1,jm/2
508  rclt(j) = acos(cosclt(j))
509  ENDDO
510  DO j = 1,jm/2
511  phi = rclt(j) * degrad
512  xlat(j) = 90. - phi
513  xlat(jm-j+1) = phi - 90.
514  ENDDO
515  ELSE
516  CALL splat(0,jm,cosclt,wgtclt)
517  DO j=1,jm
518  rclt(j) = acos(cosclt(j))
519  xlat(j) = 90.0 - rclt(j) * degrad
520  ENDDO
521  ENDIF
522 !
523 c print *,' cosclt=',cosclt
524 ! print *,' RCLT(1)=',RCLT(1)
525  sumdif = 0.
526  DO j = jm/2,2,-1
527  diffx(j) = xlat(j) - xlat(j-1)
528  sumdif = sumdif + diffx(j)
529  ENDDO
530  avedif=sumdif/(float(jm/2))
531 ! print *,' XLAT= avedif: ',avedif
532 ! write (6,107) (xlat(J)-xlat(j-1),J=JM,2,-1)
533  print *,' XLAT='
534  write (6,106) (xlat(j),j=jm,1,-1)
535  106 format( 10(f7.3,1x))
536  107 format( 10(f9.5,1x))
537 C
538  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
539 C
540  DO j=1,jmn
541  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
542  ENDDO
543  print *,
544  & ' Before GICE ZAVG(1,2)=',zavg(1,2),zslm(1,2)
545  print *,
546  & ' Before GICE ZAVG(1,12)=',zavg(1,12),zslm(1,12)
547  print *,
548  & ' Before GICE ZAVG(1,52)=',zavg(1,52),zslm(1,52)
549  print *,
550  & ' Before GICE ZAVG(1,112)=',zavg(1,jmn-112),zslm(1,112)
551 ! GICE: Grumbine 30" Antarctica orog IMNx3616 from S to N & wraped E-W.
552 ! NB: Zfields are S to N and W-E!
553  iosg = 0
554  READ(15,iostat=iosg) gice
555  if(iosg .ne. 0 ) then
556  print *,' *** Err on reading GICE record, iosg=',iosg
557  print *,' exec continues but NO GICE correction done '
558 ! stop
559  else
560  print *,' GICE 30" Antarctica RAMP orog 43200x3616 read OK'
561  print *,' Processing! '
562  print *,' Processing! '
563  print *,' Processing! '
564  do j = 1, 3601
565  do i = 1, imn
566  zsave1 = zavg(i,j)
567  zsave2 = zslm(i,j)
568  if( gice(i,j) .ne. -99. .and. gice(i,j) .ne. -1.0 ) then
569  if ( gice(i,j) .gt. 0.) then
570  zavg(i,j) = int( gice(i,j) + 0.5 )
571 !! --- for GICE values less than or equal to 0 (0, -1, or -99) then
572 !! --- radar-sat (RAMP) values are not valid and revert back to old orog
573  zslm(i,j) = 1
574  endif
575  endif
576 !jaa ALON = float(i-1) * 360./float(IMN)
577 !jaa ALAT = glat(j)
578 ! if( ZAVG(i,j) .ne. zsave1 .and. i .lt. 3 )
579 ! & print *,' antarctica change to ZAVG(i=',i,'j=',j,')=',
580 ! & ZAVG(i,j),ZSLM(i,j),' from originally:',zsave1,zsave2
581 ! &write(6,151)i,j,ZAVG(i,j),ZSLM(i,j),zsave1,zsave2,ALAT,ALON
582 ! 151 format(1x,'antarctica ZAVG(i=',i3,' j=',i3,')=',i5,i3,
583 ! &' orig:',i5,i3,' Lat=',f8.3,f9.3,'E')
584 !jaa if( ZAVG(i,j) .ne. zsave1 ) then
585 !jaa if ( i .le. 1201 .and. i .gt. 1200 )then
586 !jaa write(6,152)i,j,ZAVG(i,j),ZSLM(i,j),zsave1,zsave2,ALAT,ALON,
587 !jaa & GICE(i,j)
588 !jaa endif
589 !jaa endif
590  152 format(1x,' ZAVG(i=',i4,' j=',i4,')=',i5,i3,
591  &' orig:',i5,i4,' Lat=',f7.3,f8.2,'E',' GICE=',f8.1)
592  enddo
593  enddo
594  endif
595 ! print *,
596 ! & ' After GICE ZAVG(1,2)=',ZAVG(1,2),ZSLM(1,2)
597 ! print *,
598 ! & ' After GICE ZAVG(1,12)=',ZAVG(1,12),ZSLM(1,12)
599 ! print *,
600 ! & ' After GICE ZAVG(1,52)=',ZAVG(1,52),ZSLM(1,52)
601 ! print *,
602 ! & ' After GICE ZAVG(1,112)=',ZAVG(1,112),ZSLM(1,112)
603 !C
604 C COMPUTE MOUNTAIN DATA : ORO SLM VAR (Std Dev) OC
605 C
606 ! --- The coupled ocean model is already on a Guasian grid if (IM,JM)
607 ! --- Attempt to Open the file if mskocn=1
608  istat=0
609  if (mskocn .eq. 1) then
610 ! open(25,form='unformatted',iostat=istat)
611 ! open(25,form='binary',iostat=istat)
612 ! --- open to fort.25 with link to file in script
613  open(25,form='formatted',iostat=istat)
614  if (istat.ne.0) then
615  mskocn = 0
616  print *,' Ocean lsm file Open failure: mskocn,istat=',mskocn,istat
617  else
618  mskocn = 1
619  print *,' Ocean lsm file Opened OK: mskocn,istat=',mskocn,istat
620  endif
621 ! --- Read it in
622  ios=0
623  oclsm=0.
624 ! read(25,iostat=ios)OCLSM
625  read(25,*,iostat=ios)oclsm
626  if (ios.ne.0) then
627  mskocn = 0
628 ! --- did not properly read Gaussian grid ocean land-sea mask, but
629 ! continue using ZSLMX
630  print *,' Rd fail: Ocean lsm - continue, mskocn,ios=',mskocn,ios
631  else
632  mskocn = 1
633  print *,' Rd OK: ocean lsm: mskocn,ios=',mskocn,ios
634 ! --- LSM initialized to ocean mask especially for case where Ocean
635 ! --- changed by ocean model to land to cope with its problems
636 ! --- remember, that lake mask is in zslm to be assigned in MAKEMT.
637  if ( mskocn .eq. 1 ) then
638  DO j = 1,jm
639  DO i = 1,numi(j)
640  if ( notocn .eq. 0 ) then
641  slmi(i,j) = float(nint(oclsm(i,j)))
642  else
643  if ( nint(oclsm(i,j)) .eq. 0) then
644  slmi(i,j) = 1
645  else
646  slmi(i,j) = 0
647  endif
648  endif
649  enddo
650  enddo
651  print *,' OCLSM',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
652  print *,' SLMI:',slmi(1,1),slmi(50,50),slmi(75,75),slmi(im,jm)
653 ! --- Diag
654 ! WRITE(27,iostat=ios) REAL(SLMI,4)
655 ! print *,' write SLMI/OCLSM diag input:',ios
656  endif
657  endif
658 
659  else
660  print *,' Not using Ocean model land sea mask'
661  endif
662 
663  if (mskocn .eq. 1)then
664  print *,' LSM:',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
665  endif
666 
667 !--- reading grid file.
668  grid_from_file = .false.
669  is_south_pole = .false.
670  is_north_pole = .false.
671  i_south_pole = 0
672  j_south_pole = 0
673  i_north_pole = 0
674  j_north_pole = 0
675  if( trim(outgrid) .NE. "none" ) then
676  grid_from_file = .true.
677  inquire(file=trim(outgrid), exist=fexist)
678  if(.not. fexist) then
679  print*, "file "//trim(outgrid)//" does not exist"
680  CALL errexit(4)
681  endif
682  do ncid = 103, 512
683  inquire( ncid,opened=opened )
684  if( .NOT.opened )exit
685  end do
686 
687  print*, "outgrid=", trim(outgrid)
688  error=nf__open(trim(outgrid),nf_nowrite,fsize,ncid)
689  call netcdf_err(error, 'Open file '//trim(outgrid) )
690  error=nf_inq_dimid(ncid, 'nx', id_dim)
691  call netcdf_err(error, 'inquire dimension nx from file '//
692  & trim(outgrid) )
693  nx = 2*im
694  ny = 2*jm
695 ! error=nf_inq_dimlen(ncid,id_dim,nx)
696 ! print*, "nx = ", nx, id_dim
697 ! call netcdf_err(error, 'inquire dimension nx length '//
698 ! & 'from file '//trim(OUTGRID) )
699 ! error=nf_inq_dimid(ncid, 'ny', id_dim)
700 ! call netcdf_err(error, 'inquire dimension ny from file '//
701 ! & trim(OUTGRID) )
702 ! error=nf_inq_dimlen(ncid,id_dim,ny)
703 ! call netcdf_err(error, 'inquire dimension ny length '//
704 ! & 'from file '//trim(OUTGRID) )
705 ! IM should equal nx/2 and JM should equal ny/2
706 ! if(IM .ne. nx/2) then
707 ! print*, "IM=",IM, " /= grid file nx/2=",nx/2
708 ! CALL ERREXIT(4)
709 ! endif
710 ! if(JM .ne. ny/2) then
711 ! print*, "JM=",JM, " /= grid file ny/2=",ny/2
712 ! CALL ERREXIT(4)
713 ! endif
714  print*, "Read the grid from file "//trim(outgrid)
715 
716  allocate(tmpvar(nx+1,ny+1))
717 
718  error=nf_inq_varid(ncid, 'x', id_var)
719  call netcdf_err(error, 'inquire varid of x from file '
720  & //trim(outgrid) )
721  error=nf_get_var_double(ncid, id_var, tmpvar)
722  call netcdf_err(error, 'inquire data of x from file '
723  & //trim(outgrid) )
724  !--- adjust lontitude to be between 0 and 360.
725  do j = 1,ny+1; do i = 1,nx+1
726  if(tmpvar(i,j) .NE. missing_value) then
727  if(tmpvar(i,j) .GT. 360) tmpvar(i,j) = tmpvar(i,j) - 360
728  if(tmpvar(i,j) .LT. 0) tmpvar(i,j) = tmpvar(i,j) + 360
729  endif
730  enddo; enddo
731 
732  geolon(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
733  geolon_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
734 
735  error=nf_inq_varid(ncid, 'y', id_var)
736  call netcdf_err(error, 'inquire varid of y from file '
737  & //trim(outgrid) )
738  error=nf_get_var_double(ncid, id_var, tmpvar)
739  call netcdf_err(error, 'inquire data of y from file '
740  & //trim(outgrid) )
741  geolat(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
742  geolat_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
743 
744  !--- figure out pole location.
745  maxlat = -90
746  minlat = 90
747  i_north_pole = 0
748  j_north_pole = 0
749  i_south_pole = 0
750  j_south_pole = 0
751  do j = 1, ny+1; do i = 1, nx+1
752  if( tmpvar(i,j) > maxlat ) then
753  i_north_pole=i
754  j_north_pole=j
755  maxlat = tmpvar(i,j)
756  endif
757  if( tmpvar(i,j) < minlat ) then
758  i_south_pole=i
759  j_south_pole=j
760  minlat = tmpvar(i,j)
761  endif
762  enddo ; enddo
763  !--- only when maxlat is close to 90. the point is north pole
764  if(maxlat < 89.9 ) then
765  i_north_pole = 0
766  j_north_pole = 0
767  endif
768  if(minlat > -89.9 ) then
769  i_south_pole = 0
770  j_south_pole = 0
771  endif
772  print*, "minlat=", minlat, "maxlat=", maxlat
773  print*, "north pole supergrid index is ",
774  & i_north_pole, j_north_pole
775  print*, "south pole supergrid index is ",
776  & i_south_pole, j_south_pole
777  deallocate(tmpvar)
778 
779  if(i_south_pole >0 .and. j_south_pole > 0) then
780  if(mod(i_south_pole,2)==0) then ! stretched grid
781  do j = 1, jm; do i = 1, im
782  if(i==i_south_pole/2 .and. (j==j_south_pole/2
783  & .or. j==j_south_pole/2+1) ) then
784  is_south_pole(i,j) = .true.
785  print*, "south pole at i,j=", i, j
786  endif
787  enddo; enddo
788  else
789  do j = 1, jm; do i = 1, im
790  if((i==i_south_pole/2 .or. i==i_south_pole/2+1)
791  & .and. (j==j_south_pole/2 .or.
792  & j==j_south_pole/2+1) ) then
793  is_south_pole(i,j) = .true.
794  print*, "south pole at i,j=", i, j
795  endif
796  enddo; enddo
797  endif
798  endif
799  if(i_north_pole >0 .and. j_north_pole > 0) then
800  if(mod(i_north_pole,2)==0) then ! stretched grid
801  do j = 1, jm; do i = 1, im
802  if(i==i_north_pole/2 .and. (j==j_north_pole/2 .or.
803  & j==j_north_pole/2+1) ) then
804  is_north_pole(i,j) = .true.
805  print*, "north pole at i,j=", i, j
806  endif
807  enddo; enddo
808  else
809  do j = 1, jm; do i = 1, im
810  if((i==i_north_pole/2 .or. i==i_north_pole/2+1)
811  & .and. (j==j_north_pole/2 .or.
812  & j==j_north_pole/2+1) ) then
813  is_north_pole(i,j) = .true.
814  print*, "north pole at i,j=", i, j
815  endif
816  enddo; enddo
817  endif
818  endif
819 
820 
821  allocate(tmpvar(nx,ny))
822  error=nf_inq_varid(ncid, 'area', id_var)
823  call netcdf_err(error, 'inquire varid of area from file '
824  & //trim(outgrid) )
825  error=nf_get_var_double(ncid, id_var, tmpvar)
826  call netcdf_err(error, 'inquire data of area from file '
827  & //trim(outgrid) )
828 
829  do j = 1, jm
830  do i = 1, im
831  dx(i,j) = sqrt(tmpvar(2*i-1,2*j-1)+tmpvar(2*i,2*j-1)
832  & +tmpvar(2*i-1,2*j )+tmpvar(2*i,2*j ))
833  dy(i,j) = dx(i,j)
834  enddo
835  enddo
836 ! allocate(tmpvar(nx,ny+1))
837 
838 ! error=nf_inq_varid(ncid, 'dx', id_var)
839 ! call netcdf_err(error, 'inquire varid of dx from file '
840 ! & //trim(OUTGRID) )
841 ! error=nf_get_var_double(ncid, id_var, tmpvar)
842 ! call netcdf_err(error, 'inquire data of dx from file '
843 ! & //trim(OUTGRID) )
844 ! dx(1:IM,1:JM+1) = tmpvar(2:nx:2,1:ny+1:2)
845 ! deallocate(tmpvar)
846 
847 ! allocate(tmpvar(nx+1,ny))
848 ! error=nf_inq_varid(ncid, 'dy', id_var)
849 ! call netcdf_err(error, 'inquire varid of dy from file '
850 ! & //trim(OUTGRID) )
851 ! error=nf_get_var_double(ncid, id_var, tmpvar)
852 ! call netcdf_err(error, 'inquire data of dy from file '
853 ! & //trim(OUTGRID) )
854 ! dy(1:IM+1,1:JM) = tmpvar(1:nx+1:2,2:ny:2)
855 ! deallocate(tmpvar)
856  endif
857  tend=timef()
858  write(6,*)' Timer 1 time= ',tend-tbeg
859  !
860  if(grid_from_file) then
861  tbeg=timef()
862  CALL makemt2(zavg,zslm,oro,slm,land_frac,var,var4,glat,
863  & im,jm,imn,jmn,geolon_c,geolat_c)
864  tend=timef()
865  write(6,*)' MAKEMT2 time= ',tend-tbeg
866  else
867  CALL makemt(zavg,zslm,oro,slm,var,var4,glat,
868  & ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
869  endif
870 
871  call minmxj(im,jm,oro,' ORO')
872  call minmxj(im,jm,slm,' SLM')
873  call minmxj(im,jm,var,' VAR')
874  call minmxj(im,jm,var4,' VAR4')
875 C --- check for nands in above
876 ! call nanc(ORO,IM*JM,"MAKEMT_ORO")
877 ! call nanc(SLM,IM*JM,"MAKEMT_SLM")
878 ! call nanc(VAR,IM*JM,"MAKEMT_VAR")
879 ! call nanc(VAR4,IM*JM,"MAKEMT_VAR4")
880 !
881 C check antarctic pole
882 ! DO J = 1,JM
883 ! DO I = 1,numi(j)
884 ! if ( i .le. 100 .and. i .ge. 1 )then
885 ! if (j .ge. JM-1 )then
886 ! if (height .eq. 0.) print *,'I,J,SLM:',I,J,SLM(I,J)
887 ! write(6,153)i,j,ORO(i,j),HEIGHT,SLM(i,j)
888 ! endif
889 ! endif
890 ! ENDDO
891 ! ENDDO
892 C
893 C === Compute mtn principal coord HTENSR: THETA,GAMMA,SIGMA
894 C
895  if(grid_from_file) then
896  tbeg=timef()
897  CALL makepc2(zavg,zslm,theta,gamma,sigma,glat,
898  1 im,jm,imn,jmn,geolon_c,geolat_c)
899  tend=timef()
900  write(6,*)' MAKEPC2 time= ',tend-tbeg
901  else
902  CALL makepc(zavg,zslm,theta,gamma,sigma,glat,
903  1 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
904  endif
905 
906  call minmxj(im,jm,theta,' THETA')
907  call minmxj(im,jm,gamma,' GAMMA')
908  call minmxj(im,jm,sigma,' SIGMA')
909 
910 C --- check for nands in above
911 ! call nanc(THETA,IM*JM,"MAKEPC_THETA")
912 ! call nanc(GAMMA,IM*JM,"MAKEPC_GAMMA")
913 ! call nanc(SIGMA,IM*JM,"MAKEPC_SIGMA")
914 C
915 C COMPUTE MOUNTAIN DATA : OA OL
916 C
917  call minmxj(im,jm,oro,' ORO')
918  print*, "inputorog=", trim(inputorog)
919  if(grid_from_file) then
920  if(trim(inputorog) == "none") then
921  print*, "calling MAKEOA2 to compute OA, OL"
922  tbeg=timef()
923  CALL makeoa2(zavg,zslm,var,glat,oa,ol,iwork,elvmax,oro,
924  1 work1,work2,work3,work4,work5,work6,
925  2 im,jm,imn,jmn,geolon_c,geolat_c,
926  3 geolon,geolat,dx,dy,is_south_pole,is_north_pole)
927  tend=timef()
928  write(6,*)' MAKEOA2 time= ',tend-tbeg
929  else
930  !-- read the data from INPUTOROG file.
931  error=nf__open(trim(inputorog),nf_nowrite,fsize,ncid)
932  call netcdf_err(error, 'Open file '//trim(inputorog) )
933  error=nf_inq_dimid(ncid, 'lon', id_dim)
934  call netcdf_err(error, 'inquire dimension lon from file '//
935  & trim(inputorog) )
936  error=nf_inq_dimlen(ncid,id_dim,nx_in)
937  call netcdf_err(error, 'inquire dimension lon length '//
938  & 'from file '//trim(inputorog) )
939  error=nf_inq_dimid(ncid, 'lat', id_dim)
940  call netcdf_err(error, 'inquire dimension lat from file '//
941  & trim(inputorog) )
942  error=nf_inq_dimlen(ncid,id_dim,ny_in)
943  call netcdf_err(error, 'inquire dimension lat length '//
944  & 'from file '//trim(inputorog) )
945 
946  print*, "extrapolate OA, OL from Gaussian grid with nx=",
947  & nx_in, ", ny=", ny_in
948  allocate(oa_in(nx_in,ny_in,4), ol_in(nx_in,ny_in,4))
949  allocate(slm_in(nx_in,ny_in) )
950  allocate(lon_in(nx_in,ny_in), lat_in(nx_in,ny_in) )
951 
952  error=nf_inq_varid(ncid, 'oa1', id_var)
953  call netcdf_err(error, 'inquire varid of oa1 from file '
954  & //trim(inputorog) )
955  error=nf_get_var_double(ncid, id_var, oa_in(:,:,1))
956  call netcdf_err(error, 'inquire data of oa1 from file '
957  & //trim(inputorog) )
958  error=nf_inq_varid(ncid, 'oa2', id_var)
959  call netcdf_err(error, 'inquire varid of oa2 from file '
960  & //trim(inputorog) )
961  error=nf_get_var_double(ncid, id_var, oa_in(:,:,2))
962  call netcdf_err(error, 'inquire data of oa2 from file '
963  & //trim(inputorog) )
964  error=nf_inq_varid(ncid, 'oa3', id_var)
965  call netcdf_err(error, 'inquire varid of oa3 from file '
966  & //trim(inputorog) )
967  error=nf_get_var_double(ncid, id_var, oa_in(:,:,3))
968  call netcdf_err(error, 'inquire data of oa3 from file '
969  & //trim(inputorog) )
970  error=nf_inq_varid(ncid, 'oa4', id_var)
971  call netcdf_err(error, 'inquire varid of oa4 from file '
972  & //trim(inputorog) )
973  error=nf_get_var_double(ncid, id_var, oa_in(:,:,4))
974  call netcdf_err(error, 'inquire data of oa4 from file '
975  & //trim(inputorog) )
976 
977  error=nf_inq_varid(ncid, 'ol1', id_var)
978  call netcdf_err(error, 'inquire varid of ol1 from file '
979  & //trim(inputorog) )
980  error=nf_get_var_double(ncid, id_var, ol_in(:,:,1))
981  call netcdf_err(error, 'inquire data of ol1 from file '
982  & //trim(inputorog) )
983  error=nf_inq_varid(ncid, 'ol2', id_var)
984  call netcdf_err(error, 'inquire varid of ol2 from file '
985  & //trim(inputorog) )
986  error=nf_get_var_double(ncid, id_var, ol_in(:,:,2))
987  call netcdf_err(error, 'inquire data of ol2 from file '
988  & //trim(inputorog) )
989  error=nf_inq_varid(ncid, 'ol3', id_var)
990  call netcdf_err(error, 'inquire varid of ol3 from file '
991  & //trim(inputorog) )
992  error=nf_get_var_double(ncid, id_var, ol_in(:,:,3))
993  call netcdf_err(error, 'inquire data of ol3 from file '
994  & //trim(inputorog) )
995  error=nf_inq_varid(ncid, 'ol4', id_var)
996  call netcdf_err(error, 'inquire varid of ol4 from file '
997  & //trim(inputorog) )
998  error=nf_get_var_double(ncid, id_var, ol_in(:,:,4))
999  call netcdf_err(error, 'inquire data of ol4 from file '
1000  & //trim(inputorog) )
1001 
1002  error=nf_inq_varid(ncid, 'slmsk', id_var)
1003  call netcdf_err(error, 'inquire varid of slmsk from file '
1004  & //trim(inputorog) )
1005  error=nf_get_var_double(ncid, id_var, slm_in)
1006  call netcdf_err(error, 'inquire data of slmsk from file '
1007  & //trim(inputorog) )
1008 
1009  error=nf_inq_varid(ncid, 'geolon', id_var)
1010  call netcdf_err(error, 'inquire varid of geolon from file '
1011  & //trim(inputorog) )
1012  error=nf_get_var_double(ncid, id_var, lon_in)
1013  call netcdf_err(error, 'inquire data of geolon from file '
1014  & //trim(inputorog) )
1015 
1016  error=nf_inq_varid(ncid, 'geolat', id_var)
1017  call netcdf_err(error, 'inquire varid of geolat from file '
1018  & //trim(inputorog) )
1019  error=nf_get_var_double(ncid, id_var, lat_in)
1020  call netcdf_err(error, 'inquire data of geolat from file '
1021  & //trim(inputorog) )
1022 
1023  ! set slmsk=2 to be ocean (0)
1024  do j=1,ny_in; do i=1,nx_in
1025  if(slm_in(i,j) == 2) slm_in(i,j) = 0
1026  enddo; enddo
1027 
1028  error=nf_close(ncid)
1029  call netcdf_err(error, 'close file '
1030  & //trim(inputorog) )
1031 
1032  print*, "calling MAKEOA3 to compute OA, OL"
1033  CALL makeoa3(zavg,zslm,var,glat,oa,ol,iwork,elvmax,oro,slm,
1034  1 work1,work2,work3,work4,work5,work6,
1035  2 im,jm,imn,jmn,geolon_c,geolat_c,
1036  3 geolon,geolat,is_south_pole,is_north_pole,nx_in,ny_in,
1037  4 oa_in,ol_in,slm_in,lon_in,lat_in)
1038  endif
1039  else
1040  CALL makeoa(zavg,var,glat,oa,ol,iwork,elvmax,oro,
1041  1 work1,work2,work3,work4,
1042  2 work5,work6,
1043  3 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
1044  endif
1045  tbeg=timef()
1046  call minmxj(im,jm,oa,' OA')
1047  call minmxj(im,jm,ol,' OL')
1048  call minmxj(im,jm,elvmax,' ELVMAX')
1049  call minmxj(im,jm,oro,' ORO')
1050 C --- check for nands in above
1051 ! call nanc(OA(1,1,1), IM*JM,"MAKEOA_OA(1,1,1)")
1052 ! call nanc(OA(1,1,2), IM*JM,"MAKEOA_OA(1,1,2)")
1053 ! call nanc(OA(1,1,3), IM*JM,"MAKEOA_OA(1,1,3)")
1054 ! call nanc(OA(1,1,4), IM*JM,"MAKEOA_OA(1,1,4)")
1055 ! call nanc(OL(1,1,1), IM*JM,"MAKEOA_OL(1,1,1)")
1056 ! call nanc(OL(1,1,2), IM*JM,"MAKEOA_OL(1,1,2)")
1057 ! call nanc(OL(1,1,3), IM*JM,"MAKEOA_OL(1,1,3)")
1058 ! call nanc(OL(1,1,4), IM*JM,"MAKEOA_OL(1,1,4)")
1059 ! call nanc(ELVMAX, IM*JM,"MAKEPC_ELVMAX")
1060 
1061  maxc3 = 0
1062  maxc4 = 0
1063  maxc5 = 0
1064  maxc6 = 0
1065  maxc7 = 0
1066  maxc8 = 0
1067  DO j = 1,jm
1068  DO i = 1,numi(j)
1069  if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
1070  if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
1071  if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
1072  if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
1073  if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
1074  if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
1075  ENDDO
1076  ENDDO
1077  print *,' MAXC3:',maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
1078 !
1079 c itest=151
1080 c jtest=56
1081 C
1082  print *,' ===> Replacing ELVMAX with ELVMAX-ORO <=== '
1083  print *,' ===> if ELVMAX<=ORO replace with proxy <=== '
1084  print *,' ===> the sum of mean orog (ORO) and std dev <=== '
1085  DO j = 1,jm
1086  DO i = 1,numi(j)
1087  if (elvmax(i,j) .lt. oro(i,j) ) then
1088 C--- subtracting off ORO leaves std dev (this should never happen)
1089  elvmax(i,j) = max( 3. * var(i,j),0.)
1090  else
1091  elvmax(i,j) = max( elvmax(i,j) - oro(i,j),0.)
1092  endif
1093  ENDDO
1094  ENDDO
1095  maxc3 = 0
1096  maxc4 = 0
1097  maxc5 = 0
1098  maxc6 = 0
1099  maxc7 = 0
1100  maxc8 = 0
1101  DO j = 1,jm
1102  DO i = 1,numi(j)
1103  if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
1104  if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
1105  if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
1106  if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
1107  if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
1108  if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
1109  ENDDO
1110  ENDDO
1111  print *,' after MAXC 3-6 km:',maxc3,maxc4,maxc5,maxc6
1112 c
1113  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1114 ! if (JM .gt. 0) stop
1115 C
1116 C ZERO OVER OCEAN
1117 C
1118  print *,' Testing at point (itest,jtest)=',itest,jtest
1119  print *,' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
1120  print *,' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
1121  DO j = 1,jm
1122  DO i = 1,numi(j)
1123  IF(slm(i,j).EQ.0.) THEN
1124 C VAR(I,J) = 0.
1125  var4(i,j) = 0.
1126  oa(i,j,1) = 0.
1127  oa(i,j,2) = 0.
1128  oa(i,j,3) = 0.
1129  oa(i,j,4) = 0.
1130  ol(i,j,1) = 0.
1131  ol(i,j,2) = 0.
1132  ol(i,j,3) = 0.
1133  ol(i,j,4) = 0.
1134 C THETA(I,J) =0.
1135 C GAMMA(I,J) =0.
1136 C SIGMA(I,J) =0.
1137 C ELVMAX(I,J)=0.
1138 ! --- the sub-grid scale parameters for mtn blocking and gwd retain
1139 ! --- properties even if over ocean but there is elevation within the
1140 ! --- gaussian grid box.
1141  ENDIF
1142  ENDDO
1143  ENDDO
1144 C
1145 ! --- if mskocn=1 ocean land sea mask given, =0 if not present
1146 ! --- OCLSM is real(*4) array with fractional values possible
1147 ! --- 0 is ocean and 1 is land for slm
1148 ! --- Step 1: Only change SLM after GFS SLM is applied
1149 ! --- SLM is only field that will be altered by OCLSM
1150 ! --- Ocean land sea mask ocean points made ocean in atm model
1151 ! --- Land and Lakes and all other atm elv moments remain unchanged.
1152  if ( mskocn .eq. 1 ) then
1153 
1154  DO j = 1,jm
1155  DO i = 1,numi(j)
1156  if (abs(oro(i,j)) .lt. 1. ) then
1157  slm(i,j) = slmi(i,j)
1158  else
1159  if ( slmi(i,j) .eq. 1. .and. slm(i,j) .eq. 1) slm(i,j) = 1
1160  if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
1161  if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 1) slm(i,j) = 0
1162  if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
1163  endif
1164  enddo
1165  enddo
1166  endif
1167  print *,' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
1168  print *,' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
1169 
1170 C REMOVE ISOLATED POINTS
1171  DO j=2,jm-1
1172  jn=j-1
1173  js=j+1
1174  rn=REAL(numi(jn))/REAL(numi(j))
1175  rs=REAL(numi(js))/REAL(numi(j))
1176  DO i=1,numi(j)
1177  iw=mod(i+im-2,im)+1
1178  ie=mod(i,im)+1
1179  slma=slm(iw,j)+slm(ie,j)
1180  oroa=oro(iw,j)+oro(ie,j)
1181  vara=var(iw,j)+var(ie,j)
1182  var4a=var4(iw,j)+var4(ie,j)
1183  DO k=1,4
1184  oaa(k)=oa(iw,j,k)+oa(ie,j,k)
1185 ! --- (*j*) fix typo:
1186  ola(k)=ol(iw,j,k)+ol(ie,j,k)
1187  ENDDO
1188  wgta=2
1189  xn=rn*(i-1)+1
1190  IF(abs(xn-nint(xn)).LT.1.e-2) THEN
1191  in=mod(nint(xn)-1,numi(jn))+1
1192  inw=mod(in+numi(jn)-2,numi(jn))+1
1193  ine=mod(in,numi(jn))+1
1194  slma=slma+slm(inw,jn)+slm(in,jn)+slm(ine,jn)
1195  oroa=oroa+oro(inw,jn)+oro(in,jn)+oro(ine,jn)
1196  vara=vara+var(inw,jn)+var(in,jn)+var(ine,jn)
1197  var4a=var4a+var4(inw,jn)+var4(in,jn)+var4(ine,jn)
1198  DO k=1,4
1199  oaa(k)=oaa(k)+oa(inw,jn,k)+oa(in,jn,k)+oa(ine,jn,k)
1200  ola(k)=ola(k)+ol(inw,jn,k)+ol(in,jn,k)+ol(ine,jn,k)
1201  ENDDO
1202  wgta=wgta+3
1203  ELSE
1204  inw=int(xn)
1205  ine=mod(inw,numi(jn))+1
1206  slma=slma+slm(inw,jn)+slm(ine,jn)
1207  oroa=oroa+oro(inw,jn)+oro(ine,jn)
1208  vara=vara+var(inw,jn)+var(ine,jn)
1209  var4a=var4a+var4(inw,jn)+var4(ine,jn)
1210  DO k=1,4
1211  oaa(k)=oaa(k)+oa(inw,jn,k)+oa(ine,jn,k)
1212  ola(k)=ola(k)+ol(inw,jn,k)+ol(ine,jn,k)
1213  ENDDO
1214  wgta=wgta+2
1215  ENDIF
1216  xs=rs*(i-1)+1
1217  IF(abs(xs-nint(xs)).LT.1.e-2) THEN
1218  is=mod(nint(xs)-1,numi(js))+1
1219  isw=mod(is+numi(js)-2,numi(js))+1
1220  ise=mod(is,numi(js))+1
1221  slma=slma+slm(isw,js)+slm(is,js)+slm(ise,js)
1222  oroa=oroa+oro(isw,js)+oro(is,js)+oro(ise,js)
1223  vara=vara+var(isw,js)+var(is,js)+var(ise,js)
1224  var4a=var4a+var4(isw,js)+var4(is,js)+var4(ise,js)
1225  DO k=1,4
1226  oaa(k)=oaa(k)+oa(isw,js,k)+oa(is,js,k)+oa(ise,js,k)
1227  ola(k)=ola(k)+ol(isw,js,k)+ol(is,js,k)+ol(ise,js,k)
1228  ENDDO
1229  wgta=wgta+3
1230  ELSE
1231  isw=int(xs)
1232  ise=mod(isw,numi(js))+1
1233  slma=slma+slm(isw,js)+slm(ise,js)
1234  oroa=oroa+oro(isw,js)+oro(ise,js)
1235  vara=vara+var(isw,js)+var(ise,js)
1236  var4a=var4a+var4(isw,js)+var4(ise,js)
1237  DO k=1,4
1238  oaa(k)=oaa(k)+oa(isw,js,k)+oa(ise,js,k)
1239  ola(k)=ola(k)+ol(isw,js,k)+ol(ise,js,k)
1240  ENDDO
1241  wgta=wgta+2
1242  ENDIF
1243  oroa=oroa/wgta
1244  vara=vara/wgta
1245  var4a=var4a/wgta
1246  DO k=1,4
1247  oaa(k)=oaa(k)/wgta
1248  ola(k)=ola(k)/wgta
1249  ENDDO
1250  IF(slm(i,j).EQ.0..AND.slma.EQ.wgta) THEN
1251  print '("SEA ",2F8.0," MODIFIED TO LAND",2F8.0," AT ",2I8)',
1252  & oro(i,j),var(i,j),oroa,vara,i,j
1253  slm(i,j)=1.
1254  oro(i,j)=oroa
1255  var(i,j)=vara
1256  var4(i,j)=var4a
1257  DO k=1,4
1258  oa(i,j,k)=oaa(k)
1259  ol(i,j,k)=ola(k)
1260  ENDDO
1261  ELSEIF(slm(i,j).EQ.1..AND.slma.EQ.0.) THEN
1262  print '("LAND",2F8.0," MODIFIED TO SEA ",2F8.0," AT ",2I8)',
1263  & oro(i,j),var(i,j),oroa,vara,i,j
1264  slm(i,j)=0.
1265  oro(i,j)=oroa
1266  var(i,j)=vara
1267  var4(i,j)=var4a
1268  DO k=1,4
1269  oa(i,j,k)=oaa(k)
1270  ol(i,j,k)=ola(k)
1271  ENDDO
1272  ENDIF
1273  ENDDO
1274  ENDDO
1275 C--- print for testing after isolated points removed
1276  print *,' after isolated points removed'
1277  call minmxj(im,jm,oro,' ORO')
1278 C print *,' JM=',JM,' numi=',numi
1279  print *,' ORO(itest,jtest)=',oro(itest,jtest)
1280  print *,' VAR(itest,jtest)=',var(itest,jtest)
1281  print *,' VAR4(itest,jtest)=',var4(itest,jtest)
1282  print *,' OA(itest,jtest,1)=',oa(itest,jtest,1)
1283  print *,' OA(itest,jtest,2)=',oa(itest,jtest,2)
1284  print *,' OA(itest,jtest,3)=',oa(itest,jtest,3)
1285  print *,' OA(itest,jtest,4)=',oa(itest,jtest,4)
1286  print *,' OL(itest,jtest,1)=',ol(itest,jtest,1)
1287  print *,' OL(itest,jtest,2)=',ol(itest,jtest,2)
1288  print *,' OL(itest,jtest,3)=',ol(itest,jtest,3)
1289  print *,' OL(itest,jtest,4)=',ol(itest,jtest,4)
1290  print *,' Testing at point (itest,jtest)=',itest,jtest
1291  print *,' THETA(itest,jtest)=',theta(itest,jtest)
1292  print *,' GAMMA(itest,jtest)=',gamma(itest,jtest)
1293  print *,' SIGMA(itest,jtest)=',sigma(itest,jtest)
1294  print *,' ELVMAX(itest,jtest)=',elvmax(itest,jtest)
1295  print *,' EFAC=',efac
1296 C
1297  DO j=1,jm
1298  DO i=1,numi(j)
1299  oro(i,j) = oro(i,j) + efac*var(i,j)
1300  hprime(i,j,1) = var(i,j)
1301  hprime(i,j,2) = var4(i,j)
1302  hprime(i,j,3) = oa(i,j,1)
1303  hprime(i,j,4) = oa(i,j,2)
1304  hprime(i,j,5) = oa(i,j,3)
1305  hprime(i,j,6) = oa(i,j,4)
1306  hprime(i,j,7) = ol(i,j,1)
1307  hprime(i,j,8) = ol(i,j,2)
1308  hprime(i,j,9) = ol(i,j,3)
1309  hprime(i,j,10)= ol(i,j,4)
1310  hprime(i,j,11)= theta(i,j)
1311  hprime(i,j,12)= gamma(i,j)
1312  hprime(i,j,13)= sigma(i,j)
1313  hprime(i,j,14)= elvmax(i,j)
1314  ENDDO
1315  ENDDO
1316 !
1317  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1318 ! --- Quadratic filter applied by default.
1319 ! --- NF0 is normally set to an even value beyond the previous truncation,
1320 ! --- for example, for jcap=382, NF0=254+2
1321 ! --- NF1 is set as jcap+2 (and/or nearest even), eg., for t382, NF1=382+2=384
1322 ! --- if no filter is desired then NF1=NF0=0 and ORF=ORO
1323 ! --- if no filter but spectral to grid (with gibbs) then NF1=jcap+2, and NF1=jcap+1
1324 !
1325  IF ( nf1 - nf0 .eq. 0 ) filter=.false.
1326  print *,' NF1, NF0, FILTER=',nf1,nf0,filter
1327  IF (filter) THEN
1328 C SPECTRALLY TRUNCATE AND FILTER OROGRAPHY
1329  do j=1,jm
1330  if(numi(j).lt.im) then
1331  ffj=cmplx(0.,0.)
1332  call spfft1(numi(j),im/2+1,numi(j),1,ffj,oro(1,j),-1)
1333  call spfft1(im,im/2+1,im,1,ffj,oro(1,j),+1)
1334  endif
1335  enddo
1336  CALL sptez(nr,nm,4,im,jm,ors,oro,-1)
1337 !
1338  print *,' about to apply spectral filter '
1339  fff=1./(nf1-nf0)**2
1340  i=0
1341  DO m=0,nm
1342  DO n=m,nm+nr*m
1343  IF(n.GT.nf0) THEN
1344  www=max(1.-fff*(n-nf0)**2,0.)
1345  ors(i+1)=ors(i+1)*www
1346  ors(i+2)=ors(i+2)*www
1347  ENDIF
1348  i=i+2
1349  ENDDO
1350  ENDDO
1351 !
1352  CALL sptez(nr,nm,4,im,jm,ors,orf,+1)
1353  do j=1,jm
1354  if(numi(j).lt.im) then
1355  call spfft1(im,im/2+1,im,1,ffj,orf(1,j),-1)
1356  call spfft1(numi(j),im/2+1,numi(j),1,ffj,orf(1,j),+1)
1357  endif
1358  enddo
1359 
1360  ELSE
1361  IF (revlat) THEN
1362  CALL revers(im, jm, numi, slm, work1)
1363  CALL revers(im, jm, numi, oro, work1)
1364  DO imt=1,nmt
1365  CALL revers(im, jm, numi, hprime(1,1,imt), work1)
1366  ENDDO
1367  ENDIF
1368  ors=0.
1369  orf=oro
1370  ENDIF
1371  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1372  print *,' ELVMAX(',itest,jtest,')=',elvmax(itest,jtest)
1373  print *,' after spectral filter is applied'
1374  call minmxj(im,jm,oro,' ORO')
1375  call minmxj(im,jm,orf,' ORF')
1376 C
1377 C USE NEAREST NEIGHBOR INTERPOLATION TO FILL FULL GRIDS
1378  call rg2gg(im,jm,numi,slm)
1379  call rg2gg(im,jm,numi,oro)
1380  call rg2gg(im,jm,numi,orf)
1381 C --- not apply to new prin coord and ELVMAX (*j*)
1382  do imt=1,10
1383  call rg2gg(im,jm,numi,hprime(1,1,imt))
1384  enddo
1385 C
1386  print *,' after nearest neighbor interpolation applied '
1387  call minmxj(im,jm,oro,' ORO')
1388  call minmxj(im,jm,orf,' ORF')
1389  call mnmxja(im,jm,elvmax,itest,jtest,' ELVMAX')
1390  print *,' ORO,ORF(itest,jtest),itest,jtest:',
1391  & oro(itest,jtest),orf(itest,jtest),itest,jtest
1392  print *,' ELVMAX(',itest,jtest,')=',elvmax(itest,jtest)
1393 
1394 
1395 C check antarctic pole
1396  DO j = 1,jm
1397  DO i = 1,numi(j)
1398  if ( i .le. 21 .and. i .ge. 1 )then
1399  if (j .eq. jm )write(6,153)i,j,oro(i,j),elvmax(i,j),slm(i,j)
1400  153 format(1x,' ORO,ELVMAX(i=',i4,' j=',i4,')=',2e14.5,f5.1)
1401  endif
1402  ENDDO
1403  ENDDO
1404  tend=timef()
1405  write(6,*)' Timer 5 time= ',tend-tbeg
1406  if (output_binary) then
1407  tbeg=timef()
1408 C OUTPUT BINARY FIELDS
1409  print *,' OUTPUT BINARY FIELDS'
1410  WRITE(51) REAL(slm,4)
1411  WRITE(52) REAL(orf,4)
1412  WRITE(53) REAL(hprime,4)
1413  WRITE(54) REAL(ors,4)
1414  WRITE(55) REAL(oro,4)
1415  WRITE(66) REAL(theta,4)
1416  WRITE(67) REAL(gamma,4)
1417  WRITE(68) REAL(sigma,4)
1418 ! --- OCLSM is real(4) write only if ocean mask is present
1419  if ( mskocn .eq. 1 ) then
1420  ios=0
1421  WRITE(27,iostat=ios) oclsm
1422  print *,' write OCLSM input:',ios
1423 ! print *,' LSM:',OCLSM(1,1),OCLSM(50,50),OCLSM(75,75),OCLSM(IM,JM)
1424  endif
1425 C
1426  call minmxj(im,jm,oro,' ORO')
1427  print *,' IM=',im,' JM=',jm,' SPECTR=',spectr
1428 C--- Test binary file output:
1429  WRITE(71) REAL(slm,4)
1430  DO imt=1,nmt
1431  WRITE(71) REAL(HPRIME(:,:,IMT),4)
1432  print *,' HPRIME(',itest,jtest,imt,')=',hprime(itest,jtest,imt)
1433  ENDDO
1434  WRITE(71) REAL(oro,4)
1435  IF (spectr) THEN
1436  WRITE(71) REAL(ORF,4) ! smoothed spectral orography!
1437  ENDIF
1438 C OUTPUT GRIB FIELDS
1439  kpds=0
1440  kpds(1)=7
1441  kpds(2)=78
1442  kpds(3)=255
1443  kpds(4)=128
1444  kpds(5)=81
1445  kpds(6)=1
1446  kpds(8)=2004
1447  kpds(9)=1
1448  kpds(10)=1
1449  kpds(13)=4
1450  kpds(15)=1
1451  kpds(16)=51
1452  kpds(17)=1
1453  kpds(18)=1
1454  kpds(19)=1
1455  kpds(21)=20
1456  kpds(22)=0
1457  kgds=0
1458  kgds(1)=4
1459  kgds(2)=im
1460  kgds(3)=jm
1461  kgds(4)=90000-180000/pi*rclt(1)
1462  kgds(6)=128
1463  kgds(7)=180000/pi*rclt(1)-90000
1464  kgds(8)=-nint(360000./im)
1465  kgds(9)=nint(360000./im)
1466  kgds(10)=jm/2
1467  kgds(20)=255
1468 ! --- SLM
1469  CALL baopen(56,'fort.56',iret)
1470  if (iret .ne. 0) print *,' BAOPEN ERROR UNIT 56: IRET=',iret
1471  CALL putgb(56,im*jm,kpds,kgds,lb,slm,iret)
1472  print *,' SLM: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1473  if (iret .ne. 0) print *,' SLM PUTGB ERROR: UNIT 56: IRET=',iret
1474  print *,' SLM: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1475 ! --- OCLSM if present
1476 ! if ( mskocn .eq. 1 ) then
1477 ! CALL BAOPEN(27,'fort.27',IRET)
1478 ! if (iret .ne. 0) print *,' OCLSM BAOPEN ERROR UNIT 27:IRET=',IRET
1479 ! CALL PUTGB(27,IM*JM,KPDS,KGDS,LB,OCLSM,IRET)
1480 ! if (iret .ne. 0) print *,' OCLSM PUTGB ERROR: UNIT 27:IRET=',IRET
1481 ! print *,' OCLSM: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
1482 ! endif
1483 
1484  kpds(5)=8
1485  IF (spectr) THEN
1486  CALL baopen(57,'fort.57',iret)
1487  CALL putgb(57,im*jm,kpds,kgds,lb,orf,iret)
1488  print *,' ORF (ORO): putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1489  ENDIF
1490 C
1491 C === write out theta (angle of land to East) using #101 (wave dir)
1492 C === [radians] and since < 1 scale adjust kpds(22)
1493 C
1494  kpds(5)=101
1495  CALL baopen(58,'fort.58',iret)
1496  CALL putgb(58,im*jm,kpds,kgds,lb,theta,iret)
1497  print *,' THETA: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1498 C
1499 C === write out (land aspect ratio or anisotropy) using #102
1500 C === (as in wind wave hgt)
1501 C
1502  kpds(22)=2
1503  kpds(5)=102
1504  CALL baopen(60,'fort.60',iret)
1505  CALL putgb(60,im*jm,kpds,kgds,lb,sigma,iret)
1506  print *,' SIGMA: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1507 C
1508 C === write out (slope parameter sigma) using #9
1509 C === (as in std hgt)
1510 C
1511  kpds(22)=1
1512  kpds(5)=103
1513  CALL baopen(59,'fort.59',iret)
1514  CALL putgb(59,im*jm,kpds,kgds,lb,gamma,iret)
1515  print *,' GAMMA: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1516 C
1517  kpds(22)=1
1518  kpds(5)=9
1519  CALL baopen(61,'fort.61',iret)
1520  CALL putgb(61,im*jm,kpds,kgds,lb,hprime,iret)
1521  print *,' HPRIME: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1522 C
1523 C
1524  kpds(22)=0
1525  kpds(5)=8
1526  CALL baopen(62,'fort.62',iret)
1527  CALL putgb(62,im*jm,kpds,kgds,lb,elvmax,iret)
1528  print *,' ELVMAX: putgb-KPDS(22,5),iret:',kpds(22),kpds(5),iret
1529  endif ! output_binary
1530 C
1531  delxn = 360./im
1532  do i=1,im
1533  xlon(i) = delxn*(i-1)
1534  enddo
1535  IF(trim(outgrid) == "none") THEN
1536  do j=1,jm
1537  do i=1,im
1538  geolon(i,j) = xlon(i)
1539  geolat(i,j) = xlat(j)
1540  enddo
1541  enddo
1542  else
1543  do j = 1, jm
1544  xlat(j) = geolat(1,j)
1545  enddo
1546  do i = 1, im
1547  xlon(i) = geolon(i,1)
1548  enddo
1549  endif
1550  tend=timef()
1551  write(6,*)' Binary output time= ',tend-tbeg
1552  tbeg=timef()
1553  CALL write_netcdf(im,jm,slm,land_frac,oro,orf,hprime,1,1,
1554  1 geolon(1:im,1:jm),geolat(1:im,1:jm), xlon,xlat)
1555  tend=timef()
1556  write(6,*)' WRITE_NETCDF time= ',tend-tbeg
1557  print *,' wrote netcdf file out.oro.tile?.nc'
1558 
1559  print *,' ===== Deallocate Arrays and ENDING MTN VAR OROG program'
1560  deallocate (zavg)
1561  deallocate (zslm)
1562  deallocate (umd)
1563  deallocate (gice)
1564  tend=timef()
1565  write(6,*)' Total runtime time= ',tend-tbeg1
1566  RETURN
1567  END
1568 
1598  SUBROUTINE makemt(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
1599  1 glat,ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
1600  dimension glat(jmn),xlat(jm)
1601 ! REAL*4 OCLSM
1602  INTEGER zavg(imn,jmn),zslm(imn,jmn)
1603  dimension oro(im,jm),slm(im,jm),var(im,jm),var4(im,jm)
1604  dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
1605  INTEGER mskocn,isave
1606  LOGICAL flag, debug
1607 C==== DATA DEBUG/.TRUE./
1608  DATA debug/.false./
1609 C
1610 ! ---- OCLSM holds the ocean (im,jm) grid
1611 ! --- mskocn=1 Use ocean model sea land mask, OK and present,
1612 ! --- mskocn=0 dont use Ocean model sea land mask, not OK, not present
1613  print *,' _____ SUBROUTINE MAKEMT '
1614 C---- GLOBAL XLAT AND XLON ( DEGREE )
1615 C
1616  jm1 = jm - 1
1617  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1618 C
1619  DO j=1,jmn
1620  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1621  ENDDO
1622 C
1623 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1624 C
1625 C (*j*) for hard wired zero offset (lambda s =0) for terr05
1626  DO j=1,jm
1627  DO i=1,numi(j)
1628  im1 = numi(j) - 1
1629  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
1630  faclon = delx / delxn
1631  ist(i,j) = faclon * float(i-1) - faclon * 0.5 + 1
1632  ien(i,j) = faclon * float(i) - faclon * 0.5 + 1
1633 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
1634 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
1635 C
1636  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
1637  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
1638 !
1639 ! if ( I .lt. 10 .and. J .ge. JM-1 )
1640 ! 1 PRINT*,' MAKEMT: I j IST IEN ',I,j,IST(I,j),IEN(I,j)
1641  ENDDO
1642 ! if ( J .ge. JM-1 ) then
1643 ! print *,' *** FACLON=',FACLON, 'numi(j=',j,')=',numi(j)
1644 ! endif
1645  ENDDO
1646  print *,' DELX=',delx,' DELXN=',delxn
1647  DO j=1,jm-1
1648  flag=.true.
1649  DO j1=1,jmn
1650  xxlat = (xlat(j)+xlat(j+1))/2.
1651  IF(flag.AND.glat(j1).GT.xxlat) THEN
1652  jst(j) = j1
1653  jen(j+1) = j1 - 1
1654  flag = .false.
1655  ENDIF
1656  ENDDO
1657 CX PRINT*, ' J JST JEN ',J,JST(J),JEN(J),XLAT(J),GLAT(J1)
1658  ENDDO
1659  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
1660  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
1661 ! PRINT*, ' JM JST JEN=',JST(JM),JEN(JM),XLAT(JM),GLAT(JMN)
1662 C
1663 C...FIRST, AVERAGED HEIGHT
1664 C
1665  DO j=1,jm
1666  DO i=1,numi(j)
1667  oro(i,j) = 0.0
1668  var(i,j) = 0.0
1669  var4(i,j) = 0.0
1670  xnsum = 0.0
1671  xland = 0.0
1672  xwatr = 0.0
1673  xl1 = 0.0
1674  xs1 = 0.0
1675  xw1 = 0.0
1676  xw2 = 0.0
1677  xw4 = 0.0
1678  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
1679  i1 = ist(i,j) + ii1 - 1
1680  IF(i1.LE.0.) i1 = i1 + imn
1681  IF(i1.GT.imn) i1 = i1 - imn
1682 ! if ( i .le. 10 .and. i .ge. 1 ) then
1683 ! if (j .eq. JM )
1684 ! &print *,' J,JST,JEN,IST,IEN,I1=',
1685 ! &J,JST(j),JEN(J),IST(I,j),IEN(I,j),I1
1686 ! endif
1687  DO j1=jst(j),jen(j)
1688  xland = xland + float(zslm(i1,j1))
1689  xwatr = xwatr + float(1-zslm(i1,j1))
1690  xnsum = xnsum + 1.
1691  height = float(zavg(i1,j1))
1692 C.........
1693  IF(height.LT.-990.) height = 0.0
1694  xl1 = xl1 + height * float(zslm(i1,j1))
1695  xs1 = xs1 + height * float(1-zslm(i1,j1))
1696  xw1 = xw1 + height
1697  xw2 = xw2 + height ** 2
1698 C check antarctic pole
1699 ! if ( i .le. 10 .and. i .ge. 1 )then
1700 ! if (j .ge. JM-1 )then
1701 C=== degub testing
1702 ! print *," I,J,I1,J1,XL1,XS1,XW1,XW2:",I,J,I1,J1,XL1,XS1,XW1,XW2
1703 ! 153 format(1x,' ORO,ELVMAX(i=',i4,' j=',i4,')=',2E14.5,3f5.1)
1704 ! endif
1705 ! endif
1706  ENDDO
1707  ENDDO
1708  IF(xnsum.GT.1.) THEN
1709 ! --- SLM initialized with OCLSM calc from all land points except ....
1710 ! --- 0 is ocean and 1 is land for slm
1711 ! --- Step 1 is to only change SLM after GFS SLM is applied
1712 
1713  slm(i,j) = float(nint(xland/xnsum))
1714  IF(slm(i,j).NE.0.) THEN
1715  oro(i,j)= xl1 / xland
1716  ELSE
1717  oro(i,j)= xs1 / xwatr
1718  ENDIF
1719  var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1720  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
1721  i1 = ist(i,j) + ii1 - 1
1722  IF(i1.LE.0.) i1 = i1 + imn
1723  IF(i1.GT.imn) i1 = i1 - imn
1724  DO j1=jst(j),jen(j)
1725  height = float(zavg(i1,j1))
1726  IF(height.LT.-990.) height = 0.0
1727  xw4 = xw4 + (height-oro(i,j)) ** 4
1728  ENDDO
1729  ENDDO
1730  IF(var(i,j).GT.1.) THEN
1731 ! if ( I .lt. 20 .and. J .ge. JM-19 ) then
1732 ! print *,'I,J,XW4,XNSUM,VAR(I,J)',I,J,XW4,XNSUM,VAR(I,J)
1733 ! endif
1734  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
1735  ENDIF
1736  ENDIF
1737  ENDDO
1738  ENDDO
1739  WRITE(6,*) "! MAKEMT ORO SLM VAR VAR4 DONE"
1740 C
1741 
1742  RETURN
1743  END
1744 
1764  SUBROUTINE get_index(IMN,JMN,npts,lonO,latO,DELXN,
1765  & jst,jen,ilist,numx)
1766  implicit none
1767  integer, intent(in) :: imn,jmn
1768  integer :: npts
1769  real, intent(in) :: lono(npts), lato(npts)
1770  real, intent(in) :: delxn
1771  integer, intent(out) :: jst,jen
1772  integer, intent(out) :: ilist(imn)
1773  integer, intent(out) :: numx
1774  real minlat,maxlat,minlon,maxlon
1775  integer i2, ii, ist, ien
1776 
1777  minlat = minval(lato)
1778  maxlat = maxval(lato)
1779  minlon = minval(lono)
1780  maxlon = maxval(lono)
1781  ist = minlon/delxn+1
1782  ien = maxlon/delxn+1
1783  jst = (minlat+90)/delxn+1
1784  jen = (maxlat+90)/delxn
1785  !--- add a few points to both ends of j-direction
1786  jst = jst - 5
1787  if(jst<1) jst = 1
1788  jen = jen + 5
1789  if(jen>jmn) jen = jmn
1790 
1791  !--- when around the pole, just search through all the points.
1792  if((jst == 1 .OR. jen == jmn) .and.
1793  & (ien-ist+1 > imn/2) )then
1794  numx = imn
1795  do i2 = 1, imn
1796  ilist(i2) = i2
1797  enddo
1798  else if( ien-ist+1 > imn/2 ) then ! cross longitude = 0
1799  !--- find the minimum that greater than IMN/2
1800  !--- and maximum that less than IMN/2
1801  ist = 0
1802  ien = imn+1
1803  do i2 = 1, npts
1804  ii = lono(i2)/delxn+1
1805  if(ii <0 .or. ii>imn) print*,"ii=",ii,imn,lono(i2),delxn
1806  if( ii < imn/2 ) then
1807  ist = max(ist,ii)
1808  else if( ii > imn/2 ) then
1809  ien = min(ien,ii)
1810  endif
1811  enddo
1812  if(ist<1 .OR. ist>imn) then
1813  print*, .or."ist<1 ist>IMN"
1814  call abort()
1815  endif
1816  if(ien<1 .OR. ien>imn) then
1817  print*, .or."iend<1 iend>IMN"
1818  call abort()
1819  endif
1820 
1821  numx = imn - ien + 1
1822  do i2 = 1, numx
1823  ilist(i2) = ien + (i2-1)
1824  enddo
1825  do i2 = 1, ist
1826  ilist(numx+i2) = i2
1827  enddo
1828  numx = numx+ist
1829  else
1830  numx = ien-ist+1
1831  do i2 = 1, numx
1832  ilist(i2) = ist + (i2-1)
1833  enddo
1834  endif
1835 
1836  END
1837 
1858  SUBROUTINE makemt2(ZAVG,ZSLM,ORO,SLM,land_frac,VAR,VAR4,
1859  1 glat,im,jm,imn,jmn,lon_c,lat_c)
1860  implicit none
1861  real, parameter :: d2r = 3.14159265358979/180.
1862  integer, parameter :: maxsum=20000000
1863  real, dimension(:), allocatable :: hgt_1d
1864  integer im, jm, imn, jmn
1865  real glat(jmn), glon(imn)
1866  INTEGER zavg(imn,jmn),zslm(imn,jmn)
1867  real land_frac(im,jm)
1868  real oro(im,jm),slm(im,jm),var(im,jm),var4(im,jm)
1869  integer ist,ien,jst, jen
1870  real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
1871  INTEGER mskocn,isave
1872  LOGICAL flag, debug
1873  real lono(4),lato(4),loni,lati
1874  real height
1875  integer jm1,i,j,nsum,ii,jj,i1,numx,i2
1876  integer ilist(imn)
1877  real delxn,xnsum,xland,xwatr,xl1,xs1,xw1,xw2,xw4
1878 !jaa
1879  real :: xnsum_j,xland_j,xwatr_j
1880  logical inside_a_polygon
1881 C==== DATA DEBUG/.TRUE./
1882  DATA debug/.false./
1883 C
1884 ! ---- OCLSM holds the ocean (im,jm) grid
1885 ! --- mskocn=1 Use ocean model sea land mask, OK and present,
1886 ! --- mskocn=0 dont use Ocean model sea land mask, not OK, not present
1887  print *,' _____ SUBROUTINE MAKEMT2 '
1888  allocate(hgt_1d(maxsum))
1889 C---- GLOBAL XLAT AND XLON ( DEGREE )
1890 C
1891  jm1 = jm - 1
1892  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1893 C
1894  DO j=1,jmn
1895  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1896  ENDDO
1897  DO i=1,imn
1898  glon(i) = 0. + (i-1) * delxn + delxn * 0.5
1899  ENDDO
1900 
1901  land_frac(:,:) = 0.0
1902 C
1903 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1904 C
1905 C (*j*) for hard wired zero offset (lambda s =0) for terr05
1906 !$omp parallel do
1907 !$omp* private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,xw2,xw4,lono,
1908 !$omp* lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,height,
1909 !$omp* hgt_1d)
1910  DO j=1,jm
1911 ! print*, "J=", J
1912  DO i=1,im
1913  oro(i,j) = 0.0
1914  var(i,j) = 0.0
1915  var4(i,j) = 0.0
1916  xnsum = 0.0
1917  xland = 0.0
1918  xwatr = 0.0
1919  nsum = 0
1920  xl1 = 0.0
1921  xs1 = 0.0
1922  xw1 = 0.0
1923  xw2 = 0.0
1924  xw4 = 0.0
1925 
1926  lono(1) = lon_c(i,j)
1927  lono(2) = lon_c(i+1,j)
1928  lono(3) = lon_c(i+1,j+1)
1929  lono(4) = lon_c(i,j+1)
1930  lato(1) = lat_c(i,j)
1931  lato(2) = lat_c(i+1,j)
1932  lato(3) = lat_c(i+1,j+1)
1933  lato(4) = lat_c(i,j+1)
1934  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
1935  do jj = jst, jen; do i2 = 1, numx
1936  ii = ilist(i2)
1937  loni = ii*delxn
1938  lati = -90 + jj*delxn
1939  if(inside_a_polygon(loni*d2r,lati*d2r,4,
1940  & lono*d2r,lato*d2r))then
1941 
1942  xland = xland + float(zslm(ii,jj))
1943  xwatr = xwatr + float(1-zslm(ii,jj))
1944  xnsum = xnsum + 1.
1945  height = float(zavg(ii,jj))
1946  nsum = nsum+1
1947  if(nsum > maxsum) then
1948  print*, "nsum is greater than MAXSUM, increase MAXSUM"
1949  call abort()
1950  endif
1951  hgt_1d(nsum) = height
1952  IF(height.LT.-990.) height = 0.0
1953  xl1 = xl1 + height * float(zslm(ii,jj))
1954  xs1 = xs1 + height * float(1-zslm(ii,jj))
1955  xw1 = xw1 + height
1956  xw2 = xw2 + height ** 2
1957  endif
1958  enddo ; enddo
1959 
1960 
1961  IF(xnsum.GT.1.) THEN
1962 ! --- SLM initialized with OCLSM calc from all land points except ....
1963 ! --- 0 is ocean and 1 is land for slm
1964 ! --- Step 1 is to only change SLM after GFS SLM is applied
1965  land_frac(i,j) = xland/xnsum
1966  slm(i,j) = float(nint(xland/xnsum))
1967  IF(slm(i,j).NE.0.) THEN
1968  oro(i,j)= xl1 / xland
1969  ELSE
1970  oro(i,j)= xs1 / xwatr
1971  ENDIF
1972  var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1973  do i1 = 1, nsum
1974  xw4 = xw4 + (hgt_1d(i1) - oro(i,j)) ** 4
1975  enddo
1976 
1977  IF(var(i,j).GT.1.) THEN
1978  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
1979  ENDIF
1980  ENDIF
1981  ENDDO
1982  ENDDO
1983 !$omp end parallel do
1984  WRITE(6,*) "! MAKEMT2 ORO SLM VAR VAR4 DONE"
1985 C
1986  deallocate(hgt_1d)
1987  RETURN
1988  END
1989 
2018  SUBROUTINE makepc(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2019  1 glat,ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
2020 C
2021 C=== PC: principal coordinates of each Z avg orog box for L&M
2022 C
2023  parameter(rearth=6.3712e+6)
2024  dimension glat(jmn),xlat(jm),deltax(jmn)
2025  INTEGER zavg(imn,jmn),zslm(imn,jmn)
2026  dimension oro(im,jm),slm(im,jm),hl(im,jm),hk(im,jm)
2027  dimension hx2(im,jm),hy2(im,jm),hxy(im,jm),hlprim(im,jm)
2028  dimension theta(im,jm),gamma(im,jm),sigma2(im,jm),sigma(im,jm)
2029  dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
2030  LOGICAL flag, debug
2031 C=== DATA DEBUG/.TRUE./
2032  DATA debug/.false./
2033 C
2034  pi = 4.0 * atan(1.0)
2035  certh = pi * rearth
2036 C---- GLOBAL XLAT AND XLON ( DEGREE )
2037 C
2038  jm1 = jm - 1
2039  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2040  deltay = certh / float(jmn)
2041  print *, 'MAKEPC: DELTAY=',deltay
2042 C
2043  DO j=1,jmn
2044  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2045  deltax(j) = deltay * cos(glat(j)*pi/180.0)
2046  ENDDO
2047 C
2048 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2049 C
2050  DO j=1,jm
2051  DO i=1,numi(j)
2052 C IM1 = numi(j) - 1
2053  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
2054  faclon = delx / delxn
2055  ist(i,j) = faclon * float(i-1) - faclon * 0.5
2056  ist(i,j) = ist(i,j) + 1
2057  ien(i,j) = faclon * float(i) - faclon * 0.5
2058 C if (debug) then
2059 C if ( I .lt. 10 .and. J .lt. 10 )
2060 C 1 PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
2061 C endif
2062 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
2063 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
2064  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2065  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2066  if (debug) then
2067  if ( i .lt. 10 .and. j .lt. 10 )
2068  1 print*, ' I j IST IEN ',i,j,ist(i,j),ien(i,j)
2069  endif
2070  IF (ien(i,j) .LT. ist(i,j))
2071  1 print *,' MAKEPC: IEN < IST: I,J,IST(I,J),IEN(I,J)',
2072  2 i,j,ist(i,j),ien(i,j)
2073  ENDDO
2074  ENDDO
2075  DO j=1,jm-1
2076  flag=.true.
2077  DO j1=1,jmn
2078  xxlat = (xlat(j)+xlat(j+1))/2.
2079  IF(flag.AND.glat(j1).GT.xxlat) THEN
2080  jst(j) = j1
2081  jen(j+1) = j1 - 1
2082  flag = .false.
2083  ENDIF
2084  ENDDO
2085  ENDDO
2086  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2087  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2088  if (debug) then
2089  print*, ' IST,IEN(1,1-numi(1,JM))',ist(1,1),ien(1,1),
2090  1 ist(numi(jm),jm),ien(numi(jm),jm), numi(jm)
2091  print*, ' JST,JEN(1,JM) ',jst(1),jen(1),jst(jm),jen(jm)
2092  endif
2093 C
2094 C... DERIVITIVE TENSOR OF HEIGHT
2095 C
2096  DO j=1,jm
2097  DO i=1,numi(j)
2098  oro(i,j) = 0.0
2099  hx2(i,j) = 0.0
2100  hy2(i,j) = 0.0
2101  hxy(i,j) = 0.0
2102  xnsum = 0.0
2103  xland = 0.0
2104  xwatr = 0.0
2105  xl1 = 0.0
2106  xs1 = 0.0
2107  xfp = 0.0
2108  yfp = 0.0
2109  xfpyfp = 0.0
2110  xfp2 = 0.0
2111  yfp2 = 0.0
2112  hl(i,j) = 0.0
2113  hk(i,j) = 0.0
2114  hlprim(i,j) = 0.0
2115  theta(i,j) = 0.0
2116  gamma(i,j) = 0.
2117  sigma2(i,j) = 0.
2118  sigma(i,j) = 0.
2119 C
2120  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2121  i1 = ist(i,j) + ii1 - 1
2122  IF(i1.LE.0.) i1 = i1 + imn
2123  IF(i1.GT.imn) i1 = i1 - imn
2124 C
2125 C=== set the rest of the indexs for ave: 2pt staggered derivitive
2126 C
2127  i0 = i1 - 1
2128  if (i1 - 1 .le. 0 ) i0 = i0 + imn
2129  if (i1 - 1 .gt. imn) i0 = i0 - imn
2130 C
2131  ip1 = i1 + 1
2132  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2133  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2134 C
2135  DO j1=jst(j),jen(j)
2136  if (debug) then
2137  if ( i1 .eq. ist(i,j) .and. j1 .eq. jst(j) )
2138  1 print*, ' J, J1,IST,JST,DELTAX,GLAT ',
2139  2 j,j1,ist(i,j),jst(j),deltax(j1),glat(j1)
2140  if ( i1 .eq. ien(i,j) .and. j1 .eq. jen(j) )
2141  1 print*, ' J, J1,IEN,JEN,DELTAX,GLAT ',
2142  2 j,j1,ien(i,j),jen(j),deltax(j1),glat(j1)
2143  endif
2144  xland = xland + float(zslm(i1,j1))
2145  xwatr = xwatr + float(1-zslm(i1,j1))
2146  xnsum = xnsum + 1.
2147 C
2148  height = float(zavg(i1,j1))
2149  hi0 = float(zavg(i0,j1))
2150  hip1 = float(zavg(ip1,j1))
2151 C
2152  IF(height.LT.-990.) height = 0.0
2153  if(hi0 .lt. -990.) hi0 = 0.0
2154  if(hip1 .lt. -990.) hip1 = 0.0
2155 C........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
2156  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2157  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2158 C
2159 ! --- not at boundaries
2160 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
2161  if ( j1 .ne. jst(jm) .and. j1 .ne. jen(1) ) then
2162  hj0 = float(zavg(i1,j1-1))
2163  hjp1 = float(zavg(i1,j1+1))
2164  if(hj0 .lt. -990.) hj0 = 0.0
2165  if(hjp1 .lt. -990.) hjp1 = 0.0
2166 C....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
2167  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2168  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2169 C
2170 C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
2171 C === the NH pole: NB J1 goes from High at NP to Low toward SP
2172 C
2173 !RAB elseif ( J1 .eq. JST(1) ) then
2174  elseif ( j1 .eq. jst(jm) ) then
2175  ijax = i1 + imn/2
2176  if (ijax .le. 0 ) ijax = ijax + imn
2177  if (ijax .gt. imn) ijax = ijax - imn
2178 C..... at N pole we stay at the same latitude j1 but cross to opp side
2179  hijax = float(zavg(ijax,j1))
2180  hi1j1 = float(zavg(i1,j1))
2181  if(hijax .lt. -990.) hijax = 0.0
2182  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2183 C....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
2184  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2185  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2186  1 / deltay )**2
2187 C
2188 C === the SH pole: NB J1 goes from High at NP to Low toward SP
2189 C
2190 !RAB elseif ( J1 .eq. JEN(JM) ) then
2191  elseif ( j1 .eq. jen(1) ) then
2192  ijax = i1 + imn/2
2193  if (ijax .le. 0 ) ijax = ijax + imn
2194  if (ijax .gt. imn) ijax = ijax - imn
2195  hijax = float(zavg(ijax,j1))
2196  hi1j1 = float(zavg(i1,j1))
2197  if(hijax .lt. -990.) hijax = 0.0
2198  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2199  if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,hijax,hi1j1
2200 C..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
2201  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2202  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2203  1 / deltay )**2
2204  endif
2205 C
2206 C === The above does an average across the pole for the bndry in j.
2207 C23456789012345678901234567890123456789012345678901234567890123456789012......
2208 C
2209  xfpyfp = xfpyfp + xfp * yfp
2210  xl1 = xl1 + height * float(zslm(i1,j1))
2211  xs1 = xs1 + height * float(1-zslm(i1,j1))
2212 C
2213 C === average the HX2, HY2 and HXY
2214 C === This will be done over all land
2215 C
2216  ENDDO
2217  ENDDO
2218 C
2219 C === HTENSR
2220 C
2221  IF(xnsum.GT.1.) THEN
2222  slm(i,j) = float(nint(xland/xnsum))
2223  IF(slm(i,j).NE.0.) THEN
2224  oro(i,j)= xl1 / xland
2225  hx2(i,j) = xfp2 / xland
2226  hy2(i,j) = yfp2 / xland
2227  hxy(i,j) = xfpyfp / xland
2228  ELSE
2229  oro(i,j)= xs1 / xwatr
2230  ENDIF
2231 C=== degub testing
2232  if (debug) then
2233  print *," I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2234  1 xland,slm(i,j)
2235  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2236  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2237  ENDIF
2238 C
2239 C === make the principal axes, theta, and the degree of anisotropy,
2240 C === and sigma2, the slope parameter
2241 C
2242  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2243  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2244  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2245  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
2246 C
2247  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) * 180.0 / pi
2248 C === for testing print out in degrees
2249 C THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
2250  ENDIF
2251  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2252  if ( sigma2(i,j) .GE. 0. ) then
2253  sigma(i,j) = sqrt(sigma2(i,j) )
2254  if (sigma2(i,j) .ne. 0. .and.
2255  & hk(i,j) .GE. hlprim(i,j) )
2256  1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2257  else
2258  sigma(i,j)=0.
2259  endif
2260  ENDIF
2261  if (debug) then
2262  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2263  1 sigma(i,j),gamma(i,j)
2264  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2265  endif
2266  ENDDO
2267  ENDDO
2268  WRITE(6,*) "! MAKE Principal Coord DONE"
2269 C
2270  RETURN
2271  END
2272 
2292  SUBROUTINE makepc2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2293  1 glat,im,jm,imn,jmn,lon_c,lat_c)
2294 C
2295 C=== PC: principal coordinates of each Z avg orog box for L&M
2296 C
2297  implicit none
2298  real, parameter :: rearth=6.3712e+6
2299  real, parameter :: d2r = 3.14159265358979/180.
2300  integer :: im,jm,imn,jmn
2301  real :: glat(jmn),deltax(jmn)
2302  INTEGER zavg(imn,jmn),zslm(imn,jmn)
2303  real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
2304  real oro(im,jm),slm(im,jm),hl(im,jm),hk(im,jm)
2305  real hx2(im,jm),hy2(im,jm),hxy(im,jm),hlprim(im,jm)
2306  real theta(im,jm),gamma(im,jm),sigma2(im,jm),sigma(im,jm)
2307  real pi,certh,delxn,deltay,xnsum,xland,xwatr,xl1,xs1
2308  real xfp,yfp,xfpyfp,xfp2,yfp2,height
2309  real hi0,hip1,hj0,hjp1,hijax,hi1j1
2310  real lono(4),lato(4),loni,lati
2311  integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
2312  integer ilist(imn)
2313  logical inside_a_polygon
2314  LOGICAL flag, debug
2315 C=== DATA DEBUG/.TRUE./
2316  DATA debug/.false./
2317 C
2318  pi = 4.0 * atan(1.0)
2319  certh = pi * rearth
2320 C---- GLOBAL XLAT AND XLON ( DEGREE )
2321 C
2322  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2323  deltay = certh / float(jmn)
2324  print *, 'MAKEPC2: DELTAY=',deltay
2325 C
2326  DO j=1,jmn
2327  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2328  deltax(j) = deltay * cos(glat(j)*d2r)
2329  ENDDO
2330 C
2331 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2332 C
2333 
2334 C... DERIVITIVE TENSOR OF HEIGHT
2335 C
2336 !$omp parallel do
2337 !$omp* private (j,i,xnsum,xland,xwatr,xl1,xs1,xfp,yfp,xfpyfp,
2338 !$omp* xfp2,yfp2,lono,lato,jst,jen,ilist,numx,j1,i2,i1,
2339 !$omp* loni,lati,i0,ip1,height,hi0,hip1,hj0,hjp1,ijax,
2340 !$omp* hijax,hi1j1)
2341  DO j=1,jm
2342 ! print*, "J=", J
2343  DO i=1,im
2344  oro(i,j) = 0.0
2345  hx2(i,j) = 0.0
2346  hy2(i,j) = 0.0
2347  hxy(i,j) = 0.0
2348  xnsum = 0.0
2349  xland = 0.0
2350  xwatr = 0.0
2351  xl1 = 0.0
2352  xs1 = 0.0
2353  xfp = 0.0
2354  yfp = 0.0
2355  xfpyfp = 0.0
2356  xfp2 = 0.0
2357  yfp2 = 0.0
2358  hl(i,j) = 0.0
2359  hk(i,j) = 0.0
2360  hlprim(i,j) = 0.0
2361  theta(i,j) = 0.0
2362  gamma(i,j) = 0.
2363  sigma2(i,j) = 0.
2364  sigma(i,j) = 0.
2365 
2366  lono(1) = lon_c(i,j)
2367  lono(2) = lon_c(i+1,j)
2368  lono(3) = lon_c(i+1,j+1)
2369  lono(4) = lon_c(i,j+1)
2370  lato(1) = lat_c(i,j)
2371  lato(2) = lat_c(i+1,j)
2372  lato(3) = lat_c(i+1,j+1)
2373  lato(4) = lat_c(i,j+1)
2374  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
2375 
2376  do j1 = jst, jen; do i2 = 1, numx
2377  i1 = ilist(i2)
2378  loni = i1*delxn
2379  lati = -90 + j1*delxn
2380  if(inside_a_polygon(loni*d2r,lati*d2r,4,
2381  & lono*d2r,lato*d2r))then
2382 
2383 C=== set the rest of the indexs for ave: 2pt staggered derivitive
2384 C
2385  i0 = i1 - 1
2386  if (i1 - 1 .le. 0 ) i0 = i0 + imn
2387  if (i1 - 1 .gt. imn) i0 = i0 - imn
2388 C
2389  ip1 = i1 + 1
2390  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2391  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2392 
2393  xland = xland + float(zslm(i1,j1))
2394  xwatr = xwatr + float(1-zslm(i1,j1))
2395  xnsum = xnsum + 1.
2396 C
2397  height = float(zavg(i1,j1))
2398  hi0 = float(zavg(i0,j1))
2399  hip1 = float(zavg(ip1,j1))
2400 C
2401  IF(height.LT.-990.) height = 0.0
2402  if(hi0 .lt. -990.) hi0 = 0.0
2403  if(hip1 .lt. -990.) hip1 = 0.0
2404 C........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
2405  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2406  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2407 C
2408 ! --- not at boundaries
2409 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
2410  if ( j1 .ne. 1 .and. j1 .ne. jmn ) then
2411  hj0 = float(zavg(i1,j1-1))
2412  hjp1 = float(zavg(i1,j1+1))
2413  if(hj0 .lt. -990.) hj0 = 0.0
2414  if(hjp1 .lt. -990.) hjp1 = 0.0
2415 C....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
2416  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2417  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2418 C
2419 C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
2420 C === the NH pole: NB J1 goes from High at NP to Low toward SP
2421 C
2422 !RAB elseif ( J1 .eq. JST(1) ) then
2423  elseif ( j1 .eq. 1 ) then
2424  ijax = i1 + imn/2
2425  if (ijax .le. 0 ) ijax = ijax + imn
2426  if (ijax .gt. imn) ijax = ijax - imn
2427 C..... at N pole we stay at the same latitude j1 but cross to opp side
2428  hijax = float(zavg(ijax,j1))
2429  hi1j1 = float(zavg(i1,j1))
2430  if(hijax .lt. -990.) hijax = 0.0
2431  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2432 C....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
2433  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2434  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2435  1 / deltay )**2
2436 C
2437 C === the SH pole: NB J1 goes from High at NP to Low toward SP
2438 C
2439 !RAB elseif ( J1 .eq. JEN(JM) ) then
2440  elseif ( j1 .eq. jmn ) then
2441  ijax = i1 + imn/2
2442  if (ijax .le. 0 ) ijax = ijax + imn
2443  if (ijax .gt. imn) ijax = ijax - imn
2444  hijax = float(zavg(ijax,j1))
2445  hi1j1 = float(zavg(i1,j1))
2446  if(hijax .lt. -990.) hijax = 0.0
2447  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2448  if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,
2449  & hijax,hi1j1
2450 C..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
2451  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2452  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2453  1 / deltay )**2
2454  endif
2455 C
2456 C === The above does an average across the pole for the bndry in j.
2457 C23456789012345678901234567890123456789012345678901234567890123456789012......
2458 C
2459  xfpyfp = xfpyfp + xfp * yfp
2460  xl1 = xl1 + height * float(zslm(i1,j1))
2461  xs1 = xs1 + height * float(1-zslm(i1,j1))
2462  ENDIF
2463 C
2464 C === average the HX2, HY2 and HXY
2465 C === This will be done over all land
2466 C
2467  ENDDO
2468  ENDDO
2469 C
2470 C === HTENSR
2471 C
2472  IF(xnsum.GT.1.) THEN
2473  slm(i,j) = float(nint(xland/xnsum))
2474  IF(slm(i,j).NE.0.) THEN
2475  oro(i,j)= xl1 / xland
2476  hx2(i,j) = xfp2 / xland
2477  hy2(i,j) = yfp2 / xland
2478  hxy(i,j) = xfpyfp / xland
2479  ELSE
2480  oro(i,j)= xs1 / xwatr
2481  ENDIF
2482 C=== degub testing
2483  if (debug) then
2484  print *," I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2485  1 xland,slm(i,j)
2486  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2487  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2488  ENDIF
2489 C
2490 C === make the principal axes, theta, and the degree of anisotropy,
2491 C === and sigma2, the slope parameter
2492 C
2493  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2494  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2495  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2496  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
2497 C
2498  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) / d2r
2499 C === for testing print out in degrees
2500 C THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
2501  ENDIF
2502  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2503  if ( sigma2(i,j) .GE. 0. ) then
2504  sigma(i,j) = sqrt(sigma2(i,j) )
2505  if (sigma2(i,j) .ne. 0. .and.
2506  & hk(i,j) .GE. hlprim(i,j) )
2507  1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2508  else
2509  sigma(i,j)=0.
2510  endif
2511  ENDIF
2512  if (debug) then
2513  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2514  1 sigma(i,j),gamma(i,j)
2515  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2516  endif
2517  ENDDO
2518  ENDDO
2519 !$omp end parallel do
2520  WRITE(6,*) "! MAKE Principal Coord DONE"
2521 C
2522  RETURN
2523  END SUBROUTINE makepc2
2524 
2566  SUBROUTINE makeoa(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2567  1 oro,oro1,xnsum,xnsum1,xnsum2,xnsum3,xnsum4,
2568  2 ist,ien,jst,jen,im,jm,imn,jmn,xlat,numi)
2569  dimension glat(jmn),xlat(jm)
2570  INTEGER zavg(imn,jmn)
2571  dimension oro(im,jm),oro1(im,jm),elvmax(im,jm),zmax(im,jm)
2572  dimension oa4(im,jm,4),ioa4(im,jm,4)
2573  dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm)
2574  dimension xnsum(im,jm),xnsum1(im,jm),xnsum2(im,jm)
2575  dimension xnsum3(im,jm),xnsum4(im,jm)
2576  dimension var(im,jm),ol(im,jm,4),numi(jm)
2577  LOGICAL flag
2578 C
2579 C---- GLOBAL XLAT AND XLON ( DEGREE )
2580 C
2581 ! --- IM1 = IM - 1 removed (not used in this sub)
2582  jm1 = jm - 1
2583  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2584 C
2585  DO j=1,jmn
2586  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2587  ENDDO
2588  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
2589 C
2590 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2591 C
2592  DO j=1,jm
2593  DO i=1,numi(j)
2594  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
2595  faclon = delx / delxn
2596 C --- minus sign here in IST and IEN as in MAKEMT!
2597  ist(i,j) = faclon * float(i-1) - faclon * 0.5
2598  ist(i,j) = ist(i,j) + 1
2599  ien(i,j) = faclon * float(i) - faclon * 0.5
2600 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
2601 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
2602  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2603  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2604 cx PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
2605  if ( i .lt. 3 .and. j .lt. 3 )
2606  1print*,' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2607  if ( i .lt. 3 .and. j .ge. jm-1 )
2608  1print*,' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2609  ENDDO
2610  ENDDO
2611  print *,'MAKEOA: DELXN,DELX,FACLON',delxn,delx,faclon
2612  print *, ' ***** ready to start JST JEN section '
2613  DO j=1,jm-1
2614  flag=.true.
2615  DO j1=1,jmn
2616 ! --- XXLAT added as in MAKEMT and in next line as well
2617  xxlat = (xlat(j)+xlat(j+1))/2.
2618  IF(flag.AND.glat(j1).GT.xxlat) THEN
2619  jst(j) = j1
2620 ! --- JEN(J+1) = J1 - 1
2621  flag = .false.
2622  if ( j .eq. 1 )
2623  1print*,' MAKEOA: XX j JST JEN ',j,jst(j),jen(j)
2624  ENDIF
2625  ENDDO
2626  if ( j .lt. 3 )
2627  1print*,' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2628  if ( j .ge. jm-2 )
2629  1print*,' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2630 C FLAG=.TRUE.
2631 C DO J1=JST(J),JMN
2632 C IF(FLAG.AND.GLAT(J1).GT.XLAT(J)) THEN
2633 C JEN(J) = J1 - 1
2634 C FLAG = .FALSE.
2635 C ENDIF
2636 C ENDDO
2637  ENDDO
2638  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2639  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2640  print *,' ***** JST(1) JEN(1) ',jst(1),jen(1)
2641  print *,' ***** JST(JM) JEN(JM) ',jst(jm),jen(jm)
2642 C
2643  DO j=1,jm
2644  DO i=1,numi(j)
2645  xnsum(i,j) = 0.0
2646  elvmax(i,j) = oro(i,j)
2647  zmax(i,j) = 0.0
2648  ENDDO
2649  ENDDO
2650 !
2651 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
2652 ! --- to JM or to JM1
2653  DO j=1,jm
2654  DO i=1,numi(j)
2655  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2656  i1 = ist(i,j) + ii1 - 1
2657 ! --- next line as in makemt (I1 not II1) (*j*) 20070701
2658  IF(i1.LE.0.) i1 = i1 + imn
2659  IF (i1 .GT. imn) i1 = i1 - imn
2660  DO j1=jst(j),jen(j)
2661  height = float(zavg(i1,j1))
2662  IF(height.LT.-990.) height = 0.0
2663  IF ( height .gt. oro(i,j) ) then
2664  if ( height .gt. zmax(i,j) )zmax(i,j) = height
2665  xnsum(i,j) = xnsum(i,j) + 1
2666  ENDIF
2667  ENDDO
2668  ENDDO
2669  if ( i .lt. 5 .and. j .ge. jm-5 ) then
2670  print *,' I,J,ORO(I,J),XNSUM(I,J),ZMAX(I,J):',
2671  1 i,j,oro(i,j),xnsum(i,j),zmax(i,j)
2672  endif
2673  ENDDO
2674  ENDDO
2675 !
2676 C.... make ELVMAX ORO from MAKEMT sub
2677 C
2678 ! --- this will make work1 array take on oro's values on return
2679  DO j=1,jm
2680  DO i=1,numi(j)
2681 
2682  oro1(i,j) = oro(i,j)
2683  elvmax(i,j) = zmax(i,j)
2684  ENDDO
2685  ENDDO
2686 C........
2687 C The MAX elev peak (no averaging)
2688 C........
2689 ! DO J=1,JM
2690 ! DO I=1,numi(j)
2691 ! DO II1 = 1, IEN(I,J) - IST(I,J) + 1
2692 ! I1 = IST(I,J) + II1 - 1
2693 ! IF(I1.LE.0.) I1 = I1 + IMN
2694 ! IF(I1.GT.IMN) I1 = I1 - IMN
2695 ! DO J1=JST(J),JEN(J)
2696 ! if ( ELVMAX(I,J) .lt. ZMAX(I1,J1))
2697 ! 1 ELVMAX(I,J) = ZMAX(I1,J1)
2698 ! ENDDO
2699 ! ENDDO
2700 ! ENDDO
2701 ! ENDDO
2702 C
2703 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
2704 C IN A GRID BOX
2705  DO j=1,jm
2706  DO i=1,numi(j)
2707  xnsum1(i,j) = 0.0
2708  xnsum2(i,j) = 0.0
2709  xnsum3(i,j) = 0.0
2710  xnsum4(i,j) = 0.0
2711  ENDDO
2712  ENDDO
2713 ! --- loop
2714  DO j=1,jm1
2715  DO i=1,numi(j)
2716  hc = 1116.2 - 0.878 * var(i,j)
2717 ! print *,' I,J,HC,VAR:',I,J,HC,VAR(I,J)
2718  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2719  i1 = ist(i,j) + ii1 - 1
2720 ! IF (I1.LE.0.) print *,' I1 less than 0',I1,II1,IST(I,J),IEN(I,J)
2721 ! if ( J .lt. 3 .or. J .gt. JM-2 ) then
2722 ! IF(I1 .GT. IMN)print *,' I1 > IMN',J,I1,II1,IMN,IST(I,J),IEN(I,J)
2723 ! endif
2724  IF(i1.GT.imn) i1 = i1 - imn
2725  DO j1=jst(j),jen(j)
2726  IF(float(zavg(i1,j1)) .GT. hc)
2727  1 xnsum1(i,j) = xnsum1(i,j) + 1
2728  xnsum2(i,j) = xnsum2(i,j) + 1
2729  ENDDO
2730  ENDDO
2731 C
2732  inci = nint((ien(i,j)-ist(i,j)) * 0.5)
2733  isttt = min(max(ist(i,j)-inci,1),imn)
2734  ieddd = min(max(ien(i,j)-inci,1),imn)
2735 C
2736  incj = nint((jen(j)-jst(j)) * 0.5)
2737  jsttt = min(max(jst(j)-incj,1),jmn)
2738  jeddd = min(max(jen(j)-incj,1),jmn)
2739 ! if ( J .lt. 3 .or. J .gt. JM-3 ) then
2740 ! if(I .lt. 3 .or. I .gt. IM-3) then
2741 ! print *,' INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD:',
2742 ! 1 I,J,INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD
2743 ! endif
2744 ! endif
2745 C
2746  DO i1=isttt,ieddd
2747  DO j1=jsttt,jeddd
2748  IF(float(zavg(i1,j1)) .GT. hc)
2749  1 xnsum3(i,j) = xnsum3(i,j) + 1
2750  xnsum4(i,j) = xnsum4(i,j) + 1
2751  ENDDO
2752  ENDDO
2753 cx print*,' i j hc var ',i,j,hc,var(i,j)
2754 cx print*,'xnsum12 ',xnsum1(i,j),xnsum2(i,j)
2755 cx print*,'xnsum34 ',xnsum3(i,j),xnsum4(i,j)
2756  ENDDO
2757  ENDDO
2758 C
2759 C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
2760 C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
2761 C (KWD = 1 2 3 4)
2762 C ( WD = W S SW NW)
2763 C
2764 C
2765  DO kwd = 1, 4
2766  DO j=1,jm
2767  DO i=1,numi(j)
2768  oa4(i,j,kwd) = 0.0
2769  ENDDO
2770  ENDDO
2771  ENDDO
2772 C
2773  DO j=1,jm-2
2774  DO i=1,numi(j)
2775  ii = i + 1
2776  IF (ii .GT. numi(j)) ii = ii - numi(j)
2777  xnpu = xnsum(i,j) + xnsum(i,j+1)
2778  xnpd = xnsum(ii,j) + xnsum(ii,j+1)
2779  IF (xnpd .NE. xnpu) oa4(ii,j+1,1) = 1. - xnpd / max(xnpu , 1.)
2780  ol(ii,j+1,1) = (xnsum3(i,j+1)+xnsum3(ii,j+1))/
2781  1 (xnsum4(i,j+1)+xnsum4(ii,j+1))
2782 ! if ( I .lt. 20 .and. J .ge. JM-19 ) then
2783 ! PRINT*,' MAKEOA: I J IST IEN ',I,j,IST(I,J),IEN(I,J)
2784 ! PRINT*,' HC VAR ',HC,VAR(i,j)
2785 ! PRINT*,' MAKEOA: XNSUM(I,J)=',XNSUM(I,J),XNPU, XNPD
2786 ! PRINT*,' MAKEOA: XNSUM3(I,J+1),XNSUM3(II,J+1)',
2787 ! 1 XNSUM3(I,J+1),XNSUM3(II,J+1)
2788 ! PRINT*,' MAKEOA: II, OA4(II,J+1,1), OL(II,J+1,1):',
2789 ! 1 II, OA4(II,J+1,1), OL(II,J+1,1)
2790 ! endif
2791  ENDDO
2792  ENDDO
2793  DO j=1,jm-2
2794  DO i=1,numi(j)
2795  ii = i + 1
2796  IF (ii .GT. numi(j)) ii = ii - numi(j)
2797  xnpu = xnsum(i,j+1) + xnsum(ii,j+1)
2798  xnpd = xnsum(i,j) + xnsum(ii,j)
2799  IF (xnpd .NE. xnpu) oa4(ii,j+1,2) = 1. - xnpd / max(xnpu , 1.)
2800  ol(ii,j+1,2) = (xnsum3(ii,j)+xnsum3(ii,j+1))/
2801  1 (xnsum4(ii,j)+xnsum4(ii,j+1))
2802  ENDDO
2803  ENDDO
2804  DO j=1,jm-2
2805  DO i=1,numi(j)
2806  ii = i + 1
2807  IF (ii .GT. numi(j)) ii = ii - numi(j)
2808  xnpu = xnsum(i,j+1) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2809  xnpd = xnsum(ii,j) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2810  IF (xnpd .NE. xnpu) oa4(ii,j+1,3) = 1. - xnpd / max(xnpu , 1.)
2811  ol(ii,j+1,3) = (xnsum1(ii,j)+xnsum1(i,j+1))/
2812  1 (xnsum2(ii,j)+xnsum2(i,j+1))
2813  ENDDO
2814  ENDDO
2815  DO j=1,jm-2
2816  DO i=1,numi(j)
2817  ii = i + 1
2818  IF (ii .GT. numi(j)) ii = ii - numi(j)
2819  xnpu = xnsum(i,j) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2820  xnpd = xnsum(ii,j+1) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2821  IF (xnpd .NE. xnpu) oa4(ii,j+1,4) = 1. - xnpd / max(xnpu , 1.)
2822  ol(ii,j+1,4) = (xnsum1(i,j)+xnsum1(ii,j+1))/
2823  1 (xnsum2(i,j)+xnsum2(ii,j+1))
2824  ENDDO
2825  ENDDO
2826 C
2827  DO kwd = 1, 4
2828  DO i=1,numi(j)
2829  ol(i,1,kwd) = ol(i,2,kwd)
2830  ol(i,jm,kwd) = ol(i,jm-1,kwd)
2831  ENDDO
2832  ENDDO
2833 C
2834  DO kwd=1,4
2835  DO j=1,jm
2836  DO i=1,numi(j)
2837  t = oa4(i,j,kwd)
2838  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
2839  ENDDO
2840  ENDDO
2841  ENDDO
2842 C
2843  ns0 = 0
2844  ns1 = 0
2845  ns2 = 0
2846  ns3 = 0
2847  ns4 = 0
2848  ns5 = 0
2849  ns6 = 0
2850  DO kwd=1,4
2851  DO j=1,jm
2852  DO i=1,numi(j)
2853  t = abs( oa4(i,j,kwd) )
2854  IF(t .EQ. 0.) THEN
2855  ioa4(i,j,kwd) = 0
2856  ns0 = ns0 + 1
2857  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
2858  ioa4(i,j,kwd) = 1
2859  ns1 = ns1 + 1
2860  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
2861  ioa4(i,j,kwd) = 2
2862  ns2 = ns2 + 1
2863  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
2864  ioa4(i,j,kwd) = 3
2865  ns3 = ns3 + 1
2866  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
2867  ioa4(i,j,kwd) = 4
2868  ns4 = ns4 + 1
2869  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
2870  ioa4(i,j,kwd) = 5
2871  ns5 = ns5 + 1
2872  ELSE IF(t .GT. 10000.) THEN
2873  ioa4(i,j,kwd) = 6
2874  ns6 = ns6 + 1
2875  ENDIF
2876  ENDDO
2877  ENDDO
2878  ENDDO
2879 C
2880  WRITE(6,*) "! MAKEOA EXIT"
2881 C
2882  RETURN
2883  END SUBROUTINE makeoa
2884 
2894  function get_lon_angle(dx,lat, DEGRAD)
2895  implicit none
2896  real dx, lat, degrad
2897 
2898  real get_lon_angle
2899  real, parameter :: radius = 6371000
2900 
2901  get_lon_angle = 2*asin( sin(dx/radius*0.5)/cos(lat) )*degrad
2902 
2903  end function get_lon_angle
2904 
2913  function get_lat_angle(dy, DEGRAD)
2914  implicit none
2915  real dy, degrad
2916 
2917  real get_lat_angle
2918  real, parameter :: radius = 6371000
2919 
2920  get_lat_angle = dy/radius*degrad
2921 
2922  end function get_lat_angle
2923 
2960  SUBROUTINE makeoa2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2961  1 oro,oro1,xnsum,xnsum1,xnsum2,xnsum3,xnsum4,
2962  2 im,jm,imn,jmn,lon_c,lat_c,lon_t,lat_t,dx,dy,
2963  3 is_south_pole,is_north_pole )
2964  implicit none
2965  real, parameter :: missing_value = -9999.
2966  real, parameter :: d2r = 3.14159265358979/180.
2967  real, PARAMETER :: r2d=180./3.14159265358979
2968  integer im,jm,imn,jmn
2969  real glat(jmn)
2970  INTEGER zavg(imn,jmn),zslm(imn,jmn)
2971  real oro(im,jm),oro1(im,jm),elvmax(im,jm),zmax(im,jm)
2972  real oa4(im,jm,4)
2973  integer ioa4(im,jm,4)
2974  real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
2975  real lon_t(im,jm), lat_t(im,jm)
2976  real dx(im,jm), dy(im,jm)
2977  logical is_south_pole(im,jm), is_north_pole(im,jm)
2978  real xnsum(im,jm),xnsum1(im,jm),xnsum2(im,jm)
2979  real xnsum3(im,jm),xnsum4(im,jm)
2980  real var(im,jm),ol(im,jm,4)
2981  LOGICAL flag
2982  integer i,j,ilist(imn),numx,i1,j1,ii1
2983  integer kwd,ii,npts
2984  real lono(4),lato(4),loni,lati
2985  real delxn,hc,height,xnpu,xnpd,t
2986  integer ns0,ns1,ns2,ns3,ns4,ns5,ns6
2987  logical inside_a_polygon
2988  real lon,lat,dlon,dlat,dlat_old
2989  real lon1,lat1,lon2,lat2
2990  real xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
2991  real hc_11, hc_12, hc_21, hc_22
2992  real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
2993  real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
2994  real get_lon_angle, get_lat_angle, get_xnsum
2995  integer ist, ien, jst, jen
2996  real xland,xwatr,xl1,xs1,oroavg,slm
2997 C
2998 C---- GLOBAL XLAT AND XLON ( DEGREE )
2999 C
3000  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
3001 C
3002  DO j=1,jmn
3003  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3004  ENDDO
3005  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
3006 C
3007 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
3008 C
3009 C
3010  DO j=1,jm
3011  DO i=1,im
3012  xnsum(i,j) = 0.0
3013  elvmax(i,j) = oro(i,j)
3014  zmax(i,j) = 0.0
3015 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
3016 C IN A GRID BOX
3017  xnsum1(i,j) = 0.0
3018  xnsum2(i,j) = 0.0
3019  xnsum3(i,j) = 0.0
3020  xnsum4(i,j) = 0.0
3021  oro1(i,j) = oro(i,j)
3022  elvmax(i,j) = zmax(i,j)
3023  ENDDO
3024  ENDDO
3025 
3026 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3027 ! --- to JM or to JM1
3028 !$omp parallel do
3029 !$omp* private (j,i,hc,lono,lato,jst,jen,ilist,numx,j1,ii1,i1,loni,
3030 !$omp* lati,height)
3031  DO j=1,jm
3032 ! print*, "J=", J
3033  DO i=1,im
3034  hc = 1116.2 - 0.878 * var(i,j)
3035  lono(1) = lon_c(i,j)
3036  lono(2) = lon_c(i+1,j)
3037  lono(3) = lon_c(i+1,j+1)
3038  lono(4) = lon_c(i,j+1)
3039  lato(1) = lat_c(i,j)
3040  lato(2) = lat_c(i+1,j)
3041  lato(3) = lat_c(i+1,j+1)
3042  lato(4) = lat_c(i,j+1)
3043  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3044  do j1 = jst, jen; do ii1 = 1, numx
3045  i1 = ilist(ii1)
3046  loni = i1*delxn
3047  lati = -90 + j1*delxn
3048  if(inside_a_polygon(loni*d2r,lati*d2r,4,
3049  & lono*d2r,lato*d2r))then
3050 
3051  height = float(zavg(i1,j1))
3052  IF(height.LT.-990.) height = 0.0
3053  IF ( height .gt. oro(i,j) ) then
3054  if ( height .gt. zmax(i,j) )zmax(i,j) = height
3055  ENDIF
3056  endif
3057  ENDDO ; ENDDO
3058  ENDDO
3059  ENDDO
3060 !$omp end parallel do
3061 C
3062 ! --- this will make work1 array take on oro's values on return
3063 ! --- this will make work1 array take on oro's values on return
3064  DO j=1,jm
3065  DO i=1,im
3066 
3067  oro1(i,j) = oro(i,j)
3068  elvmax(i,j) = zmax(i,j)
3069  ENDDO
3070  ENDDO
3071 
3072  DO kwd = 1, 4
3073  DO j=1,jm
3074  DO i=1,im
3075  oa4(i,j,kwd) = 0.0
3076  ol(i,j,kwd) = 0.0
3077  ENDDO
3078  ENDDO
3079  ENDDO
3080  !
3081 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3082 C
3083 C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
3084 C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
3085 C (KWD = 1 2 3 4)
3086 C ( WD = W S SW NW)
3087 C
3088 C
3089 !$omp parallel do
3090 !$omp* private (j,i,lon,lat,kwd,dlon,dlat,lon1,lon2,lat1,lat2,
3091 !$omp* xnsum11,xnsum12,xnsum21,xnsum22,xnpu,xnpd,
3092 !$omp* xnsum1_11,xnsum2_11,hc_11, xnsum1_12,xnsum2_12,
3093 !$omp* hc_12,xnsum1_21,xnsum2_21,hc_21, xnsum1_22,
3094 !$omp* xnsum2_22,hc_22)
3095  DO j=1,jm
3096 ! print*, "j = ", j
3097  DO i=1,im
3098  lon = lon_t(i,j)
3099  lat = lat_t(i,j)
3100  !--- for around north pole, oa and ol are all 0
3101 
3102  if(is_north_pole(i,j)) then
3103  print*, "set oa1 = 0 and ol=0 at i,j=", i,j
3104  do kwd = 1, 4
3105  oa4(i,j,kwd) = 0.
3106  ol(i,j,kwd) = 0.
3107  enddo
3108  else if(is_south_pole(i,j)) then
3109  print*, "set oa1 = 0 and ol=1 at i,j=", i,j
3110  do kwd = 1, 4
3111  oa4(i,j,kwd) = 0.
3112  ol(i,j,kwd) = 1.
3113  enddo
3114  else
3115 
3116  !--- for each point, find a lat-lon grid box with same dx and dy as the cubic grid box
3117  dlon = get_lon_angle(dx(i,j), lat*d2r, r2d )
3118  dlat = get_lat_angle(dy(i,j), r2d)
3119  !--- adjust dlat if the points are close to pole.
3120  if( lat-dlat*0.5<-90.) then
3121  print*, "at i,j =", i,j, lat, dlat, lat-dlat*0.5
3122  print*, "ERROR: lat-dlat*0.5<-90."
3123  call errexit(4)
3124  endif
3125  if( lat+dlat*2 > 90.) then
3126  dlat_old = dlat
3127  dlat = (90-lat)*0.5
3128  print*, "at i,j=",i,j," adjust dlat from ",
3129  & dlat_old, " to ", dlat
3130  endif
3131  !--- lower left
3132  lon1 = lon-dlon*1.5
3133  lon2 = lon-dlon*0.5
3134  lat1 = lat-dlat*0.5
3135  lat2 = lat+dlat*0.5
3136 
3137  if(lat1<-90 .or. lat2>90) then
3138  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3139  endif
3140  xnsum11 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3141  & zavg,zslm,delxn)
3142 
3143  !--- upper left
3144  lon1 = lon-dlon*1.5
3145  lon2 = lon-dlon*0.5
3146  lat1 = lat+dlat*0.5
3147  lat2 = lat+dlat*1.5
3148  if(lat1<-90 .or. lat2>90) then
3149  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3150  endif
3151  xnsum12 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3152  & zavg,zslm,delxn)
3153 
3154  !--- lower right
3155  lon1 = lon-dlon*0.5
3156  lon2 = lon+dlon*0.5
3157  lat1 = lat-dlat*0.5
3158  lat2 = lat+dlat*0.5
3159  if(lat1<-90 .or. lat2>90) then
3160  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3161  endif
3162  xnsum21 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3163  & zavg,zslm,delxn)
3164 
3165  !--- upper right
3166  lon1 = lon-dlon*0.5
3167  lon2 = lon+dlon*0.5
3168  lat1 = lat+dlat*0.5
3169  lat2 = lat+dlat*1.5
3170  if(lat1<-90 .or. lat2>90) then
3171  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3172  endif
3173 
3174  xnsum22 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3175  & zavg,zslm,delxn)
3176 
3177  xnpu = xnsum11 + xnsum12
3178  xnpd = xnsum21 + xnsum22
3179  IF (xnpd .NE. xnpu) oa4(i,j,1) = 1. - xnpd / max(xnpu , 1.)
3180 
3181  xnpu = xnsum11 + xnsum21
3182  xnpd = xnsum12 + xnsum22
3183  IF (xnpd .NE. xnpu) oa4(i,j,2) = 1. - xnpd / max(xnpu , 1.)
3184 
3185  xnpu = xnsum11 + (xnsum12+xnsum21)*0.5
3186  xnpd = xnsum22 + (xnsum12+xnsum21)*0.5
3187  IF (xnpd .NE. xnpu) oa4(i,j,3) = 1. - xnpd / max(xnpu , 1.)
3188 
3189  xnpu = xnsum12 + (xnsum11+xnsum22)*0.5
3190  xnpd = xnsum21 + (xnsum11+xnsum22)*0.5
3191  IF (xnpd .NE. xnpu) oa4(i,j,4) = 1. - xnpd / max(xnpu , 1.)
3192 
3193 
3194  !--- calculate OL3 and OL4
3195  !--- lower left
3196  lon1 = lon-dlon*1.5
3197  lon2 = lon-dlon*0.5
3198  lat1 = lat-dlat*0.5
3199  lat2 = lat+dlat*0.5
3200  if(lat1<-90 .or. lat2>90) then
3201  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3202  endif
3203  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3204  & zavg,zslm,delxn, xnsum1_11, xnsum2_11, hc_11)
3205 
3206  !--- upper left
3207  lon1 = lon-dlon*1.5
3208  lon2 = lon-dlon*0.5
3209  lat1 = lat+dlat*0.5
3210  lat2 = lat+dlat*1.5
3211  if(lat1<-90 .or. lat2>90) then
3212  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3213  endif
3214  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3215  & zavg,zslm,delxn, xnsum1_12, xnsum2_12, hc_12)
3216 
3217  !--- lower right
3218  lon1 = lon-dlon*0.5
3219  lon2 = lon+dlon*0.5
3220  lat1 = lat-dlat*0.5
3221  lat2 = lat+dlat*0.5
3222  if(lat1<-90 .or. lat2>90) then
3223  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3224  endif
3225  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3226  & zavg,zslm,delxn, xnsum1_21, xnsum2_21, hc_21)
3227 
3228  !--- upper right
3229  lon1 = lon-dlon*0.5
3230  lon2 = lon+dlon*0.5
3231  lat1 = lat+dlat*0.5
3232  lat2 = lat+dlat*1.5
3233  if(lat1<-90 .or. lat2>90) then
3234  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3235  endif
3236  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3237  & zavg,zslm,delxn, xnsum1_22, xnsum2_22, hc_22)
3238 
3239  ol(i,j,3) = (xnsum1_22+xnsum1_11)/(xnsum2_22+xnsum2_11)
3240  ol(i,j,4) = (xnsum1_12+xnsum1_21)/(xnsum2_12+xnsum2_21)
3241 
3242  !--- calculate OL1 and OL2
3243  !--- lower left
3244  lon1 = lon-dlon*2.0
3245  lon2 = lon-dlon
3246  lat1 = lat
3247  lat2 = lat+dlat
3248  if(lat1<-90 .or. lat2>90) then
3249  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3250  endif
3251  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3252  & zavg,zslm,delxn, xnsum1_11, xnsum2_11, hc_11)
3253 
3254  !--- upper left
3255  lon1 = lon-dlon*2.0
3256  lon2 = lon-dlon
3257  lat1 = lat+dlat
3258  lat2 = lat+dlat*2.0
3259  if(lat1<-90 .or. lat2>90) then
3260  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3261  endif
3262 
3263  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3264  & zavg,zslm,delxn, xnsum1_12, xnsum2_12, hc_12)
3265 
3266  !--- lower right
3267  lon1 = lon-dlon
3268  lon2 = lon
3269  lat1 = lat
3270  lat2 = lat+dlat
3271  if(lat1<-90 .or. lat2>90) then
3272  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3273  endif
3274  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3275  & zavg,zslm,delxn, xnsum1_21, xnsum2_21, hc_21)
3276 
3277  !--- upper right
3278  lon1 = lon-dlon
3279  lon2 = lon
3280  lat1 = lat+dlat
3281  lat2 = lat+dlat*2.0
3282  if(lat1<-90 .or. lat2>90) then
3283  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3284  endif
3285 
3286  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3287  & zavg,zslm,delxn, xnsum1_22, xnsum2_22, hc_22)
3288 
3289  ol(i,j,1) = (xnsum1_11+xnsum1_21)/(xnsum2_11+xnsum2_21)
3290  ol(i,j,2) = (xnsum1_21+xnsum1_22)/(xnsum2_21+xnsum2_22)
3291  ENDIF
3292  ENDDO
3293  ENDDO
3294 !$omp end parallel do
3295  DO kwd=1,4
3296  DO j=1,jm
3297  DO i=1,im
3298  t = oa4(i,j,kwd)
3299  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3300  ENDDO
3301  ENDDO
3302  ENDDO
3303 C
3304  ns0 = 0
3305  ns1 = 0
3306  ns2 = 0
3307  ns3 = 0
3308  ns4 = 0
3309  ns5 = 0
3310  ns6 = 0
3311  DO kwd=1,4
3312  DO j=1,jm
3313  DO i=1,im
3314  t = abs( oa4(i,j,kwd) )
3315  IF(t .EQ. 0.) THEN
3316  ioa4(i,j,kwd) = 0
3317  ns0 = ns0 + 1
3318  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
3319  ioa4(i,j,kwd) = 1
3320  ns1 = ns1 + 1
3321  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
3322  ioa4(i,j,kwd) = 2
3323  ns2 = ns2 + 1
3324  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
3325  ioa4(i,j,kwd) = 3
3326  ns3 = ns3 + 1
3327  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
3328  ioa4(i,j,kwd) = 4
3329  ns4 = ns4 + 1
3330  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
3331  ioa4(i,j,kwd) = 5
3332  ns5 = ns5 + 1
3333  ELSE IF(t .GT. 10000.) THEN
3334  ioa4(i,j,kwd) = 6
3335  ns6 = ns6 + 1
3336  ENDIF
3337  ENDDO
3338  ENDDO
3339  ENDDO
3340 C
3341  WRITE(6,*) "! MAKEOA2 EXIT"
3342 C
3343  RETURN
3344  END SUBROUTINE makeoa2
3345 
3354  function spherical_distance(theta1,phi1,theta2,phi2)
3355 
3356  real, intent(in) :: theta1, phi1, theta2, phi2
3357  real :: spherical_distance, dot
3358 
3359  if(theta1 == theta2 .and. phi1 == phi2) then
3360  spherical_distance = 0.0
3361  return
3362  endif
3363 
3364  dot = cos(phi1)*cos(phi2)*cos(theta1-theta2) + sin(phi1)*sin(phi2)
3365  if(dot > 1. ) dot = 1.
3366  if(dot < -1.) dot = -1.
3367  spherical_distance = acos(dot)
3368 
3369  return
3370 
3371  end function spherical_distance
3372 
3389  subroutine get_mismatch_index(im_in, jm_in, geolon_in,geolat_in,
3390  & bitmap_in,num_out, lon_out,lat_out, iindx, jindx )
3391  integer, intent(in) :: im_in, jm_in, num_out
3392  real, intent(in) :: geolon_in(im_in,jm_in)
3393  real, intent(in) :: geolat_in(im_in,jm_in)
3394  logical*1, intent(in) :: bitmap_in(im_in,jm_in)
3395  real, intent(in) :: lon_out(num_out), lat_out(num_out)
3396  integer, intent(out):: iindx(num_out), jindx(num_out)
3397  real, parameter :: max_dist = 1.e+20
3398  integer, parameter :: numnbr = 20
3399  integer :: i_c,j_c,jstart,jend
3400  real :: lon,lat
3401 
3402  print*, "im_in,jm_in = ", im_in, jm_in
3403  print*, "num_out = ", num_out
3404  print*, "max and min of lon_in is", minval(geolon_in),
3405  & maxval(geolon_in)
3406  print*, "max and min of lat_in is", minval(geolat_in),
3407  & maxval(geolat_in)
3408  print*, "max and min of lon_out is", minval(lon_out),
3409  & maxval(lon_out)
3410  print*, "max and min of lat_out is", minval(lat_out),
3411  & maxval(lat_out)
3412  print*, "count(bitmap_in)= ", count(bitmap_in), max_dist
3413 
3414  do n = 1, num_out
3415  ! print*, "n = ", n
3416  lon = lon_out(n)
3417  lat = lat_out(n)
3418  !--- find the j-index for the nearest point
3419  i_c = 0; j_c = 0
3420  do j = 1, jm_in-1
3421  if(lat .LE. geolat_in(1,j) .and.
3422  & lat .GE. geolat_in(1,j+1)) then
3423  j_c = j
3424  endif
3425  enddo
3426  if(lat > geolat_in(1,1)) j_c = 1
3427  if(lat < geolat_in(1,jm_in)) j_c = jm_in
3428  ! print*, "lat =", lat, geolat_in(1,jm_in), geolat_in(1,jm_in-1)
3429  ! The input is Gaussian grid.
3430  jstart = max(j_c-numnbr,1)
3431  jend = min(j_c+numnbr,jm_in)
3432  dist = max_dist
3433  iindx(n) = 0
3434  jindx(n) = 0
3435  ! print*, "jstart, jend =", jstart, jend
3436  do j = jstart, jend; do i = 1,im_in
3437  if(bitmap_in(i,j) ) then
3438  ! print*, "bitmap_in is true"
3439  d = spherical_distance(lon_out(n),lat_out(n),
3440  & geolon_in(i,j), geolat_in(i,j))
3441  if( dist > d ) then
3442  iindx(n) = i; jindx(n) = j
3443  dist = d
3444  endif
3445  endif
3446  enddo; enddo
3447  if(iindx(n) ==0) then
3448  print*, "lon,lat=", lon,lat
3449  print*, "jstart, jend=", jstart, jend, dist
3450  print*, "ERROR in get mismatch_index: not find nearest points"
3451  call errexit(4)
3452  endif
3453  enddo
3454 
3455  end subroutine get_mismatch_index
3456 
3470  subroutine interpolate_mismatch(im_in, jm_in, data_in,
3471  & num_out, data_out, iindx, jindx)
3472  integer, intent(in) :: im_in, jm_in, num_out
3473  real, intent(in) :: data_in(im_in,jm_in)
3474  real, intent(out):: data_out(num_out)
3475  integer, intent(in) :: iindx(num_out), jindx(num_out)
3476 
3477  do n = 1, num_out
3478  data_out(n) = data_in(iindx(n),jindx(n))
3479  enddo
3480 
3481  end subroutine interpolate_mismatch
3482 
3527  SUBROUTINE makeoa3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
3528  1 oro,slm,oro1,xnsum,xnsum1,xnsum2,xnsum3,xnsum4,
3529  2 im,jm,imn,jmn,lon_c,lat_c,lon_t,lat_t,
3530  3 is_south_pole,is_north_pole,imi,jmi,oa_in,ol_in,
3531  4 slm_in,lon_in,lat_in)
3532  implicit none
3533  real, parameter :: missing_value = -9999.
3534  real, parameter :: d2r = 3.14159265358979/180.
3535  real, PARAMETER :: r2d=180./3.14159265358979
3536  integer im,jm,imn,jmn,imi,jmi
3537  real glat(jmn)
3538  INTEGER zavg(imn,jmn),zslm(imn,jmn)
3539  real slm(im,jm)
3540  real oro(im,jm),oro1(im,jm),elvmax(im,jm),zmax(im,jm)
3541  real oa4(im,jm,4)
3542  integer ioa4(im,jm,4)
3543  real oa_in(imi,jmi,4), ol_in(imi,jmi,4)
3544  real slm_in(imi,jmi)
3545  real lon_in(imi,jmi), lat_in(imi,jmi)
3546  real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
3547  real lon_t(im,jm), lat_t(im,jm)
3548  logical is_south_pole(im,jm), is_north_pole(im,jm)
3549  real xnsum(im,jm),xnsum1(im,jm),xnsum2(im,jm)
3550  real xnsum3(im,jm),xnsum4(im,jm)
3551  real var(im,jm),ol(im,jm,4)
3552  LOGICAL flag
3553  integer i,j,ilist(imn),numx,i1,j1,ii1
3554  integer kwd,ii,npts
3555  real lono(4),lato(4),loni,lati
3556  real delxn,hc,height,xnpu,xnpd,t
3557  integer ns0,ns1,ns2,ns3,ns4,ns5,ns6
3558  logical inside_a_polygon
3559  real lon,lat,dlon,dlat,dlat_old
3560  real lon1,lat1,lon2,lat2
3561  real xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
3562  real hc_11, hc_12, hc_21, hc_22
3563  real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
3564  real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
3565  real get_lon_angle, get_lat_angle, get_xnsum
3566  integer ist, ien, jst, jen
3567  real xland,xwatr,xl1,xs1,oroavg
3568  integer int_opt, ipopt(20), kgds_input(200), kgds_output(200)
3569  integer count_land_output
3570  integer ij, ijmdl_output, iret, num_mismatch_land, num
3571  integer ibo(1)
3572  logical*1, allocatable :: bitmap_input(:,:)
3573  logical*1, allocatable :: bitmap_output(:)
3574  integer, allocatable :: ijsav_land_output(:)
3575  real, allocatable :: lats_land_output(:)
3576  real, allocatable :: lons_land_output(:)
3577  real, allocatable :: output_data_land(:)
3578  real, allocatable :: lons_mismatch_output(:)
3579  real, allocatable :: lats_mismatch_output(:)
3580  real, allocatable :: data_mismatch_output(:)
3581  integer, allocatable :: iindx(:), jindx(:)
3582 C
3583 C---- GLOBAL XLAT AND XLON ( DEGREE )
3584 C
3585  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
3586 C
3587  ijmdl_output = im*jm
3588 
3589  DO j=1,jmn
3590  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3591  ENDDO
3592  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
3593 C
3594 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
3595 C
3596 C
3597  DO j=1,jm
3598  DO i=1,im
3599  xnsum(i,j) = 0.0
3600  elvmax(i,j) = oro(i,j)
3601  zmax(i,j) = 0.0
3602 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
3603 C IN A GRID BOX
3604  xnsum1(i,j) = 0.0
3605  xnsum2(i,j) = 0.0
3606  xnsum3(i,j) = 0.0
3607  xnsum4(i,j) = 0.0
3608  oro1(i,j) = oro(i,j)
3609  elvmax(i,j) = zmax(i,j)
3610  ENDDO
3611  ENDDO
3612 
3613 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3614 ! --- to JM or to JM1
3615  DO j=1,jm
3616 ! print*, "J=", J
3617  DO i=1,im
3618  hc = 1116.2 - 0.878 * var(i,j)
3619  lono(1) = lon_c(i,j)
3620  lono(2) = lon_c(i+1,j)
3621  lono(3) = lon_c(i+1,j+1)
3622  lono(4) = lon_c(i,j+1)
3623  lato(1) = lat_c(i,j)
3624  lato(2) = lat_c(i+1,j)
3625  lato(3) = lat_c(i+1,j+1)
3626  lato(4) = lat_c(i,j+1)
3627  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3628  do j1 = jst, jen; do ii1 = 1, numx
3629  i1 = ilist(ii1)
3630  loni = i1*delxn
3631  lati = -90 + j1*delxn
3632  if(inside_a_polygon(loni*d2r,lati*d2r,4,
3633  & lono*d2r,lato*d2r))then
3634 
3635  height = float(zavg(i1,j1))
3636  IF(height.LT.-990.) height = 0.0
3637  IF ( height .gt. oro(i,j) ) then
3638  if ( height .gt. zmax(i,j) )zmax(i,j) = height
3639  ENDIF
3640  endif
3641  ENDDO ; ENDDO
3642  ENDDO
3643  ENDDO
3644 
3645 C
3646 ! --- this will make work1 array take on oro's values on return
3647 ! --- this will make work1 array take on oro's values on return
3648  DO j=1,jm
3649  DO i=1,im
3650 
3651  oro1(i,j) = oro(i,j)
3652  elvmax(i,j) = zmax(i,j)
3653  ENDDO
3654  ENDDO
3655 
3656  DO kwd = 1, 4
3657  DO j=1,jm
3658  DO i=1,im
3659  oa4(i,j,kwd) = 0.0
3660  ol(i,j,kwd) = 0.0
3661  ENDDO
3662  ENDDO
3663  ENDDO
3664 
3665  !--- use the nearest point to do remapping.
3666  int_opt = 2
3667  ipopt=0
3668  kgds_input = 0
3669  kgds_input(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
3670  kgds_input(2) = imi ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
3671  kgds_input(3) = jmi ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
3672  kgds_input(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
3673  kgds_input(5) = 0 ! OCT 14-16 - LON OF ORIGIN
3674  kgds_input(6) = 128 ! OCT 17 - RESOLUTION FLAG
3675  kgds_input(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
3676  kgds_input(8) = nint(-360000./imi) ! OCT 21-23 - LON OF EXTREME POINT
3677  kgds_input(9) = nint((360.0 / float(imi))*1000.0)
3678  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
3679  kgds_input(10) = jmi /2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
3680  kgds_input(12) = 255 ! OCT 29 - RESERVED
3681  kgds_input(20) = 255 ! OCT 5 - NOT USED, SET TO 255
3682 
3683 
3684  kgds_output = -1
3685 ! KGDS_OUTPUT(1) = IDRT ! OCT 6 - TYPE OF GRID (GAUSSIAN)
3686  kgds_output(2) = im ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
3687  kgds_output(3) = jm ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
3688  kgds_output(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
3689  kgds_output(5) = 0 ! OCT 14-16 - LON OF ORIGIN
3690  kgds_output(6) = 128 ! OCT 17 - RESOLUTION FLAG
3691  kgds_output(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
3692  kgds_output(8) = nint(-360000./im) ! OCT 21-23 - LON OF EXTREME POINT
3693  kgds_output(9) = nint((360.0 / float(im))*1000.0)
3694  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
3695  kgds_output(10) = jm /2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
3696  kgds_output(12) = 255 ! OCT 29 - RESERVED
3697  kgds_output(20) = 255 ! OCT 5 - NOT USED, SET TO 255
3698 
3699  count_land_output=0
3700  print*, "sum(slm) = ", sum(slm)
3701  do ij = 1, ijmdl_output
3702  j = (ij-1)/im + 1
3703  i = mod(ij-1,im) + 1
3704  if (slm(i,j) > 0.0) then
3705  count_land_output=count_land_output+1
3706  endif
3707  enddo
3708  allocate(bitmap_input(imi,jmi))
3709  bitmap_input=.false.
3710  print*, "number of land input=", sum(slm_in)
3711  where(slm_in > 0.0) bitmap_input=.true.
3712  print*, "count(bitmap_input)", count(bitmap_input)
3713 
3714  allocate(bitmap_output(count_land_output))
3715  allocate(output_data_land(count_land_output))
3716  allocate(ijsav_land_output(count_land_output))
3717  allocate(lats_land_output(count_land_output))
3718  allocate(lons_land_output(count_land_output))
3719 
3720 
3721 
3722  count_land_output=0
3723  do ij = 1, ijmdl_output
3724  j = (ij-1)/im + 1
3725  i = mod(ij-1,im) + 1
3726  if (slm(i,j) > 0.0) then
3727  count_land_output=count_land_output+1
3728  ijsav_land_output(count_land_output)=ij
3729  lats_land_output(count_land_output)=lat_t(i,j)
3730  lons_land_output(count_land_output)=lon_t(i,j)
3731  endif
3732  enddo
3733 
3734  oa4 = 0.0
3735  ol = 0.0
3736 
3737  do kwd=1,4
3738  bitmap_output = .false.
3739  output_data_land = 0.0
3740  call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3741  & (imi*jmi), count_land_output,
3742  & 1, 1, bitmap_input, oa_in(:,:,kwd),
3743  & count_land_output, lats_land_output,
3744  & lons_land_output, ibo,
3745  & bitmap_output, output_data_land, iret)
3746  if (iret /= 0) then
3747  print*,'- ERROR IN IPOLATES ',iret
3748  call errexit(4)
3749  endif
3750 
3751  num_mismatch_land = 0
3752  do ij = 1, count_land_output
3753  if (bitmap_output(ij)) then
3754  j = (ijsav_land_output(ij)-1)/im + 1
3755  i = mod(ijsav_land_output(ij)-1,im) + 1
3756  oa4(i,j,kwd)=output_data_land(ij)
3757  else ! default value
3758  num_mismatch_land = num_mismatch_land + 1
3759  endif
3760  enddo
3761  print*, "num_mismatch_land = ", num_mismatch_land
3762 
3763  if(.not. allocated(data_mismatch_output)) then
3764  allocate(lons_mismatch_output(num_mismatch_land))
3765  allocate(lats_mismatch_output(num_mismatch_land))
3766  allocate(data_mismatch_output(num_mismatch_land))
3767  allocate(iindx(num_mismatch_land))
3768  allocate(jindx(num_mismatch_land))
3769 
3770  num = 0
3771  do ij = 1, count_land_output
3772  if (.not. bitmap_output(ij)) then
3773  num = num+1
3774  lons_mismatch_output(num) = lons_land_output(ij)
3775  lats_mismatch_output(num) = lats_land_output(ij)
3776  endif
3777  enddo
3778 
3779  ! For those land points that with bitmap_output=.false. use
3780  ! the nearest land points to interpolate.
3781  print*,"before get_mismatch_index", count(bitmap_input)
3782  call get_mismatch_index(imi,jmi,lon_in*d2r,
3783  & lat_in*d2r,bitmap_input,num_mismatch_land,
3784  & lons_mismatch_output*d2r,lats_mismatch_output*d2r,
3785  & iindx, jindx )
3786  endif
3787 
3788  data_mismatch_output = 0
3789  call interpolate_mismatch(imi,jmi,oa_in(:,:,kwd),
3790  & num_mismatch_land,data_mismatch_output,iindx,jindx)
3791 
3792  num = 0
3793  do ij = 1, count_land_output
3794  if (.not. bitmap_output(ij)) then
3795  num = num+1
3796  j = (ijsav_land_output(ij)-1)/im + 1
3797  i = mod(ijsav_land_output(ij)-1,im) + 1
3798  oa4(i,j,kwd) = data_mismatch_output(num)
3799  if(i==372 .and. j== 611) then
3800  print*, "ij=",ij, num, oa4(i,j,kwd),iindx(num),jindx(num)
3801  endif
3802  endif
3803  enddo
3804 
3805 
3806  enddo
3807 
3808  !OL
3809  do kwd=1,4
3810  bitmap_output = .false.
3811  output_data_land = 0.0
3812  call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3813  & (imi*jmi), count_land_output,
3814  & 1, 1, bitmap_input, ol_in(:,:,kwd),
3815  & count_land_output, lats_land_output,
3816  & lons_land_output, ibo,
3817  & bitmap_output, output_data_land, iret)
3818  if (iret /= 0) then
3819  print*,'- ERROR IN IPOLATES ',iret
3820  call errexit(4)
3821  endif
3822 
3823  num_mismatch_land = 0
3824  do ij = 1, count_land_output
3825  if (bitmap_output(ij)) then
3826  j = (ijsav_land_output(ij)-1)/im + 1
3827  i = mod(ijsav_land_output(ij)-1,im) + 1
3828  ol(i,j,kwd)=output_data_land(ij)
3829  else ! default value
3830  num_mismatch_land = num_mismatch_land + 1
3831  endif
3832  enddo
3833  print*, "num_mismatch_land = ", num_mismatch_land
3834 
3835  data_mismatch_output = 0
3836  call interpolate_mismatch(imi,jmi,ol_in(:,:,kwd),
3837  & num_mismatch_land,data_mismatch_output,iindx,jindx)
3838 
3839  num = 0
3840  do ij = 1, count_land_output
3841  if (.not. bitmap_output(ij)) then
3842  num = num+1
3843  j = (ijsav_land_output(ij)-1)/im + 1
3844  i = mod(ijsav_land_output(ij)-1,im) + 1
3845  ol(i,j,kwd) = data_mismatch_output(num)
3846  if(i==372 .and. j== 611) then
3847  print*, "ij=",ij, num, ol(i,j,kwd),iindx(num),jindx(num)
3848  endif
3849  endif
3850  enddo
3851 
3852 
3853  enddo
3854 
3855  deallocate(lons_mismatch_output,lats_mismatch_output)
3856  deallocate(data_mismatch_output,iindx,jindx)
3857  deallocate(bitmap_output,output_data_land)
3858  deallocate(ijsav_land_output,lats_land_output)
3859  deallocate(lons_land_output)
3860 
3861  DO kwd=1,4
3862  DO j=1,jm
3863  DO i=1,im
3864  t = oa4(i,j,kwd)
3865  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3866  ENDDO
3867  ENDDO
3868  ENDDO
3869 C
3870  ns0 = 0
3871  ns1 = 0
3872  ns2 = 0
3873  ns3 = 0
3874  ns4 = 0
3875  ns5 = 0
3876  ns6 = 0
3877  DO kwd=1,4
3878  DO j=1,jm
3879  DO i=1,im
3880  t = abs( oa4(i,j,kwd) )
3881  IF(t .EQ. 0.) THEN
3882  ioa4(i,j,kwd) = 0
3883  ns0 = ns0 + 1
3884  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
3885  ioa4(i,j,kwd) = 1
3886  ns1 = ns1 + 1
3887  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
3888  ioa4(i,j,kwd) = 2
3889  ns2 = ns2 + 1
3890  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
3891  ioa4(i,j,kwd) = 3
3892  ns3 = ns3 + 1
3893  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
3894  ioa4(i,j,kwd) = 4
3895  ns4 = ns4 + 1
3896  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
3897  ioa4(i,j,kwd) = 5
3898  ns5 = ns5 + 1
3899  ELSE IF(t .GT. 10000.) THEN
3900  ioa4(i,j,kwd) = 6
3901  ns6 = ns6 + 1
3902  ENDIF
3903  ENDDO
3904  ENDDO
3905  ENDDO
3906 C
3907  WRITE(6,*) "! MAKEOA3 EXIT"
3908 C
3909  RETURN
3910  END SUBROUTINE makeoa3
3911 
3922  SUBROUTINE revers(IM, JM, numi, F, WRK)
3923 !
3924  REAL f(im,jm), wrk(im,jm)
3925  integer numi(jm)
3926  imb2 = im / 2
3927  do i=1,im*jm
3928  wrk(i,1) = f(i,1)
3929  enddo
3930  do j=1,jm
3931  jr = jm - j + 1
3932  do i=1,im
3933  ir = i + imb2
3934  if (ir .gt. im) ir = ir - im
3935  f(ir,jr) = wrk(i,j)
3936  enddo
3937  enddo
3938 !
3939  tem = 0.0
3940  do i=1,im
3941  tem= tem + f(i,1)
3942  enddo
3943  tem = tem / im
3944  do i=1,im
3945  f(i,1) = tem
3946  enddo
3947 !
3948  RETURN
3949  END
3950 
3959  subroutine rg2gg(im,jm,numi,a)
3960  implicit none
3961  integer,intent(in):: im,jm,numi(jm)
3962  real,intent(inout):: a(im,jm)
3963  integer j,ir,ig
3964  real r,t(im)
3965  do j=1,jm
3966  r=real(numi(j))/real(im)
3967  do ig=1,im
3968  ir=mod(nint((ig-1)*r),numi(j))+1
3969  t(ig)=a(ir,j)
3970  enddo
3971  do ig=1,im
3972  a(ig,j)=t(ig)
3973  enddo
3974  enddo
3975  end subroutine
3976 
3985  subroutine gg2rg(im,jm,numi,a)
3986  implicit none
3987  integer,intent(in):: im,jm,numi(jm)
3988  real,intent(inout):: a(im,jm)
3989  integer j,ir,ig
3990  real r,t(im)
3991  do j=1,jm
3992  r=real(numi(j))/real(im)
3993  do ir=1,numi(j)
3994  ig=nint((ir-1)/r)+1
3995  t(ir)=a(ig,j)
3996  enddo
3997  do ir=1,numi(j)
3998  a(ir,j)=t(ir)
3999  enddo
4000  enddo
4001  end subroutine
4002 
4011  SUBROUTINE minmxj(IM,JM,A,title)
4012  implicit none
4013 
4014  real a(im,jm),rmin,rmax
4015  integer i,j,im,jm
4016  character*8 title
4017 
4018  rmin=1.e+10
4019  rmax=-rmin
4020 csela....................................................
4021 csela if(rmin.eq.1.e+10)return
4022 csela....................................................
4023  DO j=1,jm
4024  DO i=1,im
4025  if(a(i,j).ge.rmax)rmax=a(i,j)
4026  if(a(i,j).le.rmin)rmin=a(i,j)
4027  ENDDO
4028  ENDDO
4029  write(6,150)rmin,rmax,title
4030 150 format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
4031 C
4032  RETURN
4033  END
4034 
4046  SUBROUTINE mnmxja(IM,JM,A,imax,jmax,title)
4047  implicit none
4048 
4049  real a(im,jm),rmin,rmax
4050  integer i,j,im,jm,imax,jmax
4051  character*8 title
4052 
4053  rmin=1.e+10
4054  rmax=-rmin
4055 csela....................................................
4056 csela if(rmin.eq.1.e+10)return
4057 csela....................................................
4058  DO j=1,jm
4059  DO i=1,im
4060  if(a(i,j).ge.rmax)then
4061  rmax=a(i,j)
4062  imax=i
4063  jmax=j
4064  endif
4065  if(a(i,j).le.rmin)rmin=a(i,j)
4066  ENDDO
4067  ENDDO
4068  write(6,150)rmin,rmax,title
4069 150 format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
4070 C
4071  RETURN
4072  END
4073 
4108  SUBROUTINE spfft1(IMAX,INCW,INCG,KMAX,W,G,IDIR)
4109  IMPLICIT NONE
4110  INTEGER,INTENT(IN):: imax,incw,incg,kmax,idir
4111  COMPLEX,INTENT(INOUT):: w(incw,kmax)
4112  REAL,INTENT(INOUT):: g(incg,kmax)
4113  REAL:: aux1(25000+int(0.82*imax))
4114  REAL:: aux2(20000+int(0.57*imax))
4115  INTEGER:: naux1,naux2
4116 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4117  naux1=25000+int(0.82*imax)
4118  naux2=20000+int(0.57*imax)
4119 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4120 C FOURIER TO PHYSICAL TRANSFORM.
4121  SELECT CASE(idir)
4122  CASE(1:)
4123  SELECT CASE(digits(1.))
4124  CASE(digits(1._4))
4125  CALL scrft(1,w,incw,g,incg,imax,kmax,-1,1.,
4126  & aux1,naux1,aux2,naux2,0.,0)
4127  CALL scrft(0,w,incw,g,incg,imax,kmax,-1,1.,
4128  & aux1,naux1,aux2,naux2,0.,0)
4129  CASE(digits(1._8))
4130  CALL dcrft(1,w,incw,g,incg,imax,kmax,-1,1.,
4131  & aux1,naux1,aux2,naux2,0.,0)
4132  CALL dcrft(0,w,incw,g,incg,imax,kmax,-1,1.,
4133  & aux1,naux1,aux2,naux2,0.,0)
4134  END SELECT
4135 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4136 C PHYSICAL TO FOURIER TRANSFORM.
4137  CASE(:-1)
4138  SELECT CASE(digits(1.))
4139  CASE(digits(1._4))
4140  CALL srcft(1,g,incg,w,incw,imax,kmax,+1,1./imax,
4141  & aux1,naux1,aux2,naux2,0.,0)
4142  CALL srcft(0,g,incg,w,incw,imax,kmax,+1,1./imax,
4143  & aux1,naux1,aux2,naux2,0.,0)
4144  CASE(digits(1._8))
4145  CALL drcft(1,g,incg,w,incw,imax,kmax,+1,1./imax,
4146  & aux1,naux1,aux2,naux2,0.,0)
4147  CALL drcft(0,g,incg,w,incw,imax,kmax,+1,1./imax,
4148  & aux1,naux1,aux2,naux2,0.,0)
4149  END SELECT
4150  END SELECT
4151  END SUBROUTINE
4152 
4158  subroutine read_g(glob,ITOPO)
4159  implicit none
4160 cc
4161  integer*2 glob(360*120,180*120)
4162 cc
4163  integer ix,jx
4164  integer ia,ja
4165 cc
4166  parameter(ix=40*120,jx=50*120)
4167  parameter(ia=60*120,ja=30*120)
4168 cc
4169  integer*2 idat(ix,jx)
4170  integer itopo
4171 cc
4172  integer i,j,inttyp
4173 cc
4174  real(kind=8) dloin,dlain,rlon,rlat
4175 cc
4176  open(235, file="./fort.235", access='direct', recl=43200*21600*2)
4177  read(235,rec=1)glob
4178  close(235)
4179 cc
4180  print*,' '
4181  call maxmin(glob,360*120*180*120,'global0')
4182 cc
4183 cc
4184  dloin=1.d0/120.d0
4185  dlain=1.d0/120.d0
4186 cc
4187  rlon= -179.995833333333333333333333d0
4188  rlat= 89.995833333333333333333333d0
4189 cc
4190  inttyp=-1 ! average rectangular subset
4191 ccmr inttyp= 1 ! take closest grid point value
4192 ccmr inttyp= 0 ! interpolate from four closest grid point values
4193 cc
4194 ! call la2ga_gtopo30(glob,360*120,180*120,
4195 ! & dloin,dlain,rlon,rlat,inttyp,
4196 ! & .true.,glob,
4197 ! & 0,lonf,latg)
4198 cc
4199  return
4200  end
4201 
4209  subroutine maxmin(ia,len,tile)
4210 ccmr
4211  implicit none
4212 ccmr
4213  integer*2 ia(len)
4214  character*7 tile
4215  integer iaamax, iaamin, len, j, m, ja, kount
4216  integer(8) sum2,std,mean,isum
4217  integer i_count_notset,kount_9
4218 ! --- missing is -9999
4219 c
4220  isum = 0
4221  sum2 = 0
4222  kount = 0
4223  kount_9 = 0
4224  iaamax = -9999999
4225 ccmr iaamin = 1
4226  iaamin = 9999999
4227  i_count_notset=0
4228  do 10 m=1,len
4229  ja=ia(m)
4230 ccmr if ( ja .lt. 0 ) print *,' ja < 0:',ja
4231 ccmr if ( ja .eq. -9999 ) goto 10
4232  if ( ja .eq. -9999 ) then
4233  kount_9=kount_9+1
4234  goto 10
4235  endif
4236  if ( ja .eq. -12345 ) i_count_notset=i_count_notset+1
4237 ccmr if ( ja .eq. 0 ) goto 11
4238  iaamax = max0( iaamax, ja )
4239  iaamin = min0( iaamin, ja )
4240 ! iaamax = max0( iaamax, ia(m,j) )
4241 ! iaamin = min0( iaamin, ia(m,j) )
4242  11 continue
4243  kount = kount + 1
4244  isum = isum + ja
4245 ccmr sum2 = sum2 + ifix( float(ja) * float(ja) )
4246  sum2 = sum2 + ja*ja
4247  10 continue
4248 !
4249  mean = isum/kount
4250  std = ifix(sqrt(float((sum2/(kount))-mean**2)))
4251  print*,tile,' max=',iaamax,' min=',iaamin,' sum=',isum,
4252  & ' i_count_notset=',i_count_notset
4253  print*,tile,' mean=',mean,' std.dev=',std,
4254  & ' ko9s=',kount,kount_9,kount+kount_9
4255  return
4256  end
4257 
4267  SUBROUTINE minmaxj(IM,JM,A,title)
4268  implicit none
4269 
4270  real(kind=4) a(im,jm),rmin,rmax,undef
4271  integer i,j,im,jm,imax,jmax,imin,jmin,iundef
4272  character*8 title,chara
4273  data chara/' '/
4274  chara=title
4275  rmin=1.e+10
4276  rmax=-rmin
4277  imax=0
4278  imin=0
4279  jmax=0
4280  jmin=0
4281  iundef=0
4282  undef=-9999.
4283 csela....................................................
4284 csela if(rmin.eq.1.e+10)return
4285 csela....................................................
4286  DO j=1,jm
4287  DO i=1,im
4288  if(a(i,j).ge.rmax)then
4289  rmax=a(i,j)
4290  imax=i
4291  jmax=j
4292  endif
4293  if(a(i,j).le.rmin)then
4294  if ( a(i,j) .eq. undef ) then
4295  iundef = iundef + 1
4296  else
4297  rmin=a(i,j)
4298  imin=i
4299  jmin=j
4300  endif
4301  endif
4302  ENDDO
4303  ENDDO
4304  write(6,150)chara,rmin,imin,jmin,rmax,imax,jmax,iundef
4305 150 format(1x,a8,2x,'rmin=',e13.4,2i6,2x,'rmax=',e13.4,3i6)
4306 C
4307  RETURN
4308  END
4309 
4319  subroutine latlon2xyz(siz,lon, lat, x, y, z)
4320  implicit none
4321  integer, intent(in) :: siz
4322  real, intent(in) :: lon(siz), lat(siz)
4323  real, intent(out) :: x(siz), y(siz), z(siz)
4324 
4325  integer n
4326 
4327  do n = 1, siz
4328  x(n) = cos(lat(n))*cos(lon(n))
4329  y(n) = cos(lat(n))*sin(lon(n))
4330  z(n) = sin(lat(n))
4331  enddo
4332  end
4333 
4341  FUNCTION spherical_angle(v1, v2, v3)
4342  implicit none
4343  real, parameter :: epsln30 = 1.e-30
4344  real, parameter :: pi=3.1415926535897931
4345  real v1(3), v2(3), v3(3)
4346  real spherical_angle
4347 
4348  real px, py, pz, qx, qy, qz, ddd;
4349 
4350  ! vector product between v1 and v2
4351  px = v1(2)*v2(3) - v1(3)*v2(2)
4352  py = v1(3)*v2(1) - v1(1)*v2(3)
4353  pz = v1(1)*v2(2) - v1(2)*v2(1)
4354  ! vector product between v1 and v3
4355  qx = v1(2)*v3(3) - v1(3)*v3(2);
4356  qy = v1(3)*v3(1) - v1(1)*v3(3);
4357  qz = v1(1)*v3(2) - v1(2)*v3(1);
4358 
4359  ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
4360  if ( ddd <= 0.0 ) then
4361  spherical_angle = 0.
4362  else
4363  ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
4364  if( abs(ddd-1) < epsln30 ) ddd = 1;
4365  if( abs(ddd+1) < epsln30 ) ddd = -1;
4366  if ( ddd>1. .or. ddd<-1. ) then
4367  !FIX to correctly handle co-linear points (angle near pi or 0) */
4368  if (ddd < 0.) then
4369  spherical_angle = pi
4370  else
4371  spherical_angle = 0.
4372  endif
4373  else
4374  spherical_angle = acos( ddd )
4375  endif
4376  endif
4377 
4378  return
4379  END
4380 
4391  FUNCTION inside_a_polygon(lon1, lat1, npts, lon2, lat2)
4392  implicit none
4393  real, parameter :: epsln10 = 1.e-10
4394  real, parameter :: epsln8 = 1.e-8
4395  real, parameter :: pi=3.1415926535897931
4396  real, parameter :: range_check_criteria=0.05
4397  real :: anglesum, angle, spherical_angle
4398  integer i, ip1
4399  real lon1, lat1
4400  integer npts
4401  real lon2(npts), lat2(npts)
4402  real x2(npts), y2(npts), z2(npts)
4403  real lon1_1d(1), lat1_1d(1)
4404  real x1(1), y1(1), z1(1)
4405  real pnt0(3),pnt1(3),pnt2(3)
4406  logical inside_a_polygon
4407  real max_x2,min_x2,max_y2,min_y2,max_z2,min_z2
4408  !first convert to cartesian grid */
4409  call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
4410  lon1_1d(1) = lon1
4411  lat1_1d(1) = lat1
4412  call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
4413  inside_a_polygon = .false.
4414  max_x2 = maxval(x2)
4415  if( x1(1) > max_x2+range_check_criteria ) return
4416  min_x2 = minval(x2)
4417  if( x1(1)+range_check_criteria < min_x2 ) return
4418  max_y2 = maxval(y2)
4419  if( y1(1) > max_y2+range_check_criteria ) return
4420  min_y2 = minval(y2)
4421  if( y1(1)+range_check_criteria < min_y2 ) return
4422  max_z2 = maxval(z2)
4423  if( z1(1) > max_z2+range_check_criteria ) return
4424  min_z2 = minval(z2)
4425  if( z1(1)+range_check_criteria < min_z2 ) return
4426 
4427  pnt0(1) = x1(1)
4428  pnt0(2) = y1(1)
4429  pnt0(3) = z1(1)
4430 
4431  anglesum = 0;
4432  do i = 1, npts
4433  if(abs(x1(1)-x2(i)) < epsln10 .and.
4434  & abs(y1(1)-y2(i)) < epsln10 .and.
4435  & abs(z1(1)-z2(i)) < epsln10 ) then ! same as the corner point
4436  inside_a_polygon = .true.
4437  return
4438  endif
4439  ip1 = i+1
4440  if(ip1>npts) ip1 = 1
4441  pnt1(1) = x2(i)
4442  pnt1(2) = y2(i)
4443  pnt1(3) = z2(i)
4444  pnt2(1) = x2(ip1)
4445  pnt2(2) = y2(ip1)
4446  pnt2(3) = z2(ip1)
4447 
4448  angle = spherical_angle(pnt0, pnt2, pnt1);
4449 ! anglesum = anglesum + spherical_angle(pnt0, pnt2, pnt1);
4450  anglesum = anglesum + angle
4451  enddo
4452 
4453  if(abs(anglesum-2*pi) < epsln8) then
4454  inside_a_polygon = .true.
4455  else
4456  inside_a_polygon = .false.
4457  endif
4458 
4459  return
4460 
4461  end
4462 
4486  function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
4487  & glat,zavg,zslm,delxn)
4488  implicit none
4489 
4490  real get_xnsum
4491  logical verbose
4492  real lon1,lat1,lon2,lat2,oro,delxn
4493  integer imn,jmn
4494  real glat(jmn)
4495  integer zavg(imn,jmn),zslm(imn,jmn)
4496  integer i, j, ist, ien, jst, jen, i1
4497  real height
4498  real xland,xwatr,xl1,xs1,slm,xnsum
4499  !---figure out ist,ien,jst,jen
4500  do j = 1, jmn
4501  if( glat(j) .GT. lat1 ) then
4502  jst = j
4503  exit
4504  endif
4505  enddo
4506  do j = 1, jmn
4507  if( glat(j) .GT. lat2 ) then
4508  jen = j
4509  exit
4510  endif
4511  enddo
4512 
4513 
4514  ist = lon1/delxn + 1
4515  ien = lon2/delxn
4516  if(ist .le.0) ist = ist + imn
4517  if(ien < ist) ien = ien + imn
4518 ! if(verbose) print*, "ist,ien=",ist,ien,jst,jen
4519 
4520  !--- compute average oro
4521  oro = 0.0
4522  xnsum = 0
4523  xland = 0
4524  xwatr = 0
4525  xl1 = 0
4526  xs1 = 0
4527  do j = jst,jen
4528  do i1 = 1, ien - ist + 1
4529  i = ist + i1 -1
4530  if( i .LE. 0) i = i + imn
4531  if( i .GT. imn) i = i - imn
4532  xland = xland + float(zslm(i,j))
4533  xwatr = xwatr + float(1-zslm(i,j))
4534  xnsum = xnsum + 1.
4535  height = float(zavg(i,j))
4536  IF(height.LT.-990.) height = 0.0
4537  xl1 = xl1 + height * float(zslm(i,j))
4538  xs1 = xs1 + height * float(1-zslm(i,j))
4539  enddo
4540  enddo
4541  if( xnsum > 1.) THEN
4542  slm = float(nint(xland/xnsum))
4543  IF(slm.NE.0.) THEN
4544  oro= xl1 / xland
4545  ELSE
4546  oro = xs1 / xwatr
4547  ENDIF
4548  ENDIF
4549 
4550  get_xnsum = 0
4551  do j = jst, jen
4552  do i1= 1, ien-ist+1
4553  i = ist + i1 -1
4554  if( i .LE. 0) i = i + imn
4555  if( i .GT. imn) i = i - imn
4556  height = float(zavg(i,j))
4557  IF(height.LT.-990.) height = 0.0
4558  IF ( height .gt. oro ) get_xnsum = get_xnsum + 1
4559  enddo
4560  enddo
4561 ! if(verbose) print*, "get_xnsum=", get_xnsum, oro
4562 
4563  end function get_xnsum
4564 
4593  subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
4594  & glat,zavg,zslm,delxn,xnsum1,xnsum2,hc)
4595  implicit none
4596 
4597  real, intent(out) :: xnsum1,xnsum2,hc
4598  logical verbose
4599  real lon1,lat1,lon2,lat2,oro,delxn
4600  integer imn,jmn
4601  real glat(jmn)
4602  integer zavg(imn,jmn),zslm(imn,jmn)
4603  integer i, j, ist, ien, jst, jen, i1
4604  real height, var
4605  real xw1,xw2,slm,xnsum
4606  !---figure out ist,ien,jst,jen
4607  do j = 1, jmn
4608  if( glat(j) .GT. lat1 ) then
4609  jst = j
4610  exit
4611  endif
4612  enddo
4613  do j = 1, jmn
4614  if( glat(j) .GT. lat2 ) then
4615  jen = j
4616  exit
4617  endif
4618  enddo
4619 
4620 
4621  ist = lon1/delxn + 1
4622  ien = lon2/delxn
4623  if(ist .le.0) ist = ist + imn
4624  if(ien < ist) ien = ien + imn
4625 ! if(verbose) print*, "ist,ien=",ist,ien,jst,jen
4626 
4627  !--- compute average oro
4628  xnsum = 0
4629  xw1 = 0
4630  xw2 = 0
4631  do j = jst,jen
4632  do i1 = 1, ien - ist + 1
4633  i = ist + i1 -1
4634  if( i .LE. 0) i = i + imn
4635  if( i .GT. imn) i = i - imn
4636  xnsum = xnsum + 1.
4637  height = float(zavg(i,j))
4638  IF(height.LT.-990.) height = 0.0
4639  xw1 = xw1 + height
4640  xw2 = xw2 + height ** 2
4641  enddo
4642  enddo
4643  var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
4644  hc = 1116.2 - 0.878 * var
4645  xnsum1 = 0
4646  xnsum2 = 0
4647  do j = jst, jen
4648  do i1= 1, ien-ist+1
4649  i = ist + i1 -1
4650  if( i .LE. 0) i = i + imn
4651  if( i .GT. imn) i = i - imn
4652  height = float(zavg(i,j))
4653  IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4654  xnsum2 = xnsum2 + 1
4655  enddo
4656  enddo
4657 
4658  end subroutine get_xnsum2
4659 
4689  subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
4690  & glat,zavg,zslm,delxn,xnsum1,xnsum2,hc)
4691  implicit none
4692 
4693  real, intent(out) :: xnsum1,xnsum2
4694  real lon1,lat1,lon2,lat2,oro,delxn
4695  integer imn,jmn
4696  real glat(jmn)
4697  integer zavg(imn,jmn),zslm(imn,jmn)
4698  integer i, j, ist, ien, jst, jen, i1
4699  real height, hc
4700  real xw1,xw2,slm,xnsum
4701  !---figure out ist,ien,jst,jen
4702  ! if lat1 or lat 2 is 90 degree. set jst = JMN
4703  jst = jmn
4704  jen = jmn
4705  do j = 1, jmn
4706  if( glat(j) .GT. lat1 ) then
4707  jst = j
4708  exit
4709  endif
4710  enddo
4711  do j = 1, jmn
4712  if( glat(j) .GT. lat2 ) then
4713  jen = j
4714  exit
4715  endif
4716  enddo
4717 
4718 
4719  ist = lon1/delxn + 1
4720  ien = lon2/delxn
4721  if(ist .le.0) ist = ist + imn
4722  if(ien < ist) ien = ien + imn
4723 ! if(verbose) print*, "ist,ien=",ist,ien,jst,jen
4724 
4725  xnsum1 = 0
4726  xnsum2 = 0
4727  do j = jst, jen
4728  do i1= 1, ien-ist+1
4729  i = ist + i1 -1
4730  if( i .LE. 0) i = i + imn
4731  if( i .GT. imn) i = i - imn
4732  height = float(zavg(i,j))
4733  IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4734  xnsum2 = xnsum2 + 1
4735  enddo
4736  enddo
4737 
4738  end subroutine get_xnsum3
4739 
4753  subroutine nanc(a,l,c)
4754  integer inan1,inan2,inan3,inan4,inaq1,inaq2,inaq3,inaq4
4755  real word
4756  integer itest
4757  equivalence(itest,word)
4758 c
4759 c signaling NaN
4760  data inan1/x'7F800001'/
4761  data inan2/x'7FBFFFFF'/
4762  data inan3/x'FF800001'/
4763  data inan4/x'FFBFFFFF'/
4764 c
4765 c quiet NaN
4766 c
4767  data inaq1/x'7FC00000'/
4768  data inaq2/x'7FFFFFFF'/
4769  data inaq3/x'FFC00000'/
4770  data inaq4/x'FFFFFFFF'/
4771 c
4772  real(kind=8)a(l),rtc,t1,t2
4773  character*24 cn
4774  character*(*) c
4775 c t1=rtc()
4776 cgwv print *, ' nanc call ',c
4777  do k=1,l
4778  word=a(k)
4779  if( (itest .GE. inan1 .AND. itest .LE. inan2) .OR.
4780  * (itest .GE. inan3 .AND. itest .LE. inan4) ) then
4781  print *,' NaNs detected at word',k,' ',c
4782  return
4783  endif
4784  if( (itest .GE. inaq1 .AND. itest .LE. inaq2) .OR.
4785  * (itest .GE. inaq3 .AND. itest .LE. inaq4) ) then
4786  print *,' NaNq detected at word',k,' ',c
4787  return
4788  endif
4789 
4790  101 format(e20.10)
4791  end do
4792 c t2=rtc()
4793 cgwv print 102,l,t2-t1,c
4794  102 format(' time to check ',i9,' words is ',f10.4,' ',a24)
4795  return
4796  end
4797 
4802  real function timef()
4803  character(8) :: date
4804  character(10) :: time
4805  character(5) :: zone
4806  integer,dimension(8) :: values
4807  integer :: total
4808  real :: elapsed
4809  call date_and_time(date,time,zone,values)
4810  total=(3600*values(5))+(60*values(6))
4811  * +values(7)
4812  elapsed=float(total) + (1.0e-3*float(values(8)))
4813  timef=elapsed
4814  return
4815  end
get_xnsum3
subroutine get_xnsum3(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, zslm, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
Definition: mtnlm7_oclsm.f:4689
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:4391
spherical_distance
real function spherical_distance(theta1, phi1, theta2, phi2)
Compute a great circle distance between two points.
Definition: mtnlm7_oclsm.f:3354
tersub
subroutine tersub(IMN, JMN, IM, JM, NM, NR, NF0, NF1, NW, EFAC, BLAT, OUTGRID, INPUTOROG)
Driver routine to compute terrain.
Definition: mtnlm7_oclsm.f:185
timef
real function timef()
Get the date/time for the system clock.
Definition: mtnlm7_oclsm.f:4802
netcdf_err
subroutine netcdf_err(err, string)
Check NetCDF error code and output the error message.
Definition: netcdf_io.F90:218
latlon2xyz
subroutine latlon2xyz(siz, lon, lat, x, y, z)
Convert from latitude and longitude to x,y,z coordinates.
Definition: mtnlm7_oclsm.f:4319
spfft1
subroutine spfft1(IMAX, INCW, INCG, KMAX, W, G, IDIR)
Perform multiple fast fourier transforms.
Definition: mtnlm7_oclsm.f:4108
gg2rg
subroutine gg2rg(im, jm, numi, a)
Convert from a full grid to a reduced grid.
Definition: mtnlm7_oclsm.f:3985
rg2gg
subroutine rg2gg(im, jm, numi, a)
Convert from a reduced grid to a full grid.
Definition: mtnlm7_oclsm.f:3959
spherical_angle
real function spherical_angle(v1, v2, v3)
Compute spherical angle.
Definition: mtnlm7_oclsm.f:4341
makeoa3
subroutine makeoa3(ZAVG, zslm, 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, is_south_pole, is_north_pole, 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:3527
nanc
subroutine nanc(a, l, c)
Report NaNS and NaNQ within an address range for 8-byte real words.
Definition: mtnlm7_oclsm.f:4753
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:4267
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:2018
makemt2
subroutine makemt2(ZAVG, ZSLM, ORO, SLM, land_frac, VAR, VAR4, GLAT, IM, JM, IMN, JMN, lon_c, lat_c)
Create the orography, land-mask, land fraction, standard deviation of orography and the convexity on ...
Definition: mtnlm7_oclsm.f:1858
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:21
minmxj
subroutine minmxj(IM, JM, A, title)
Print out the maximum and minimum values of an array.
Definition: mtnlm7_oclsm.f:4011
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:2913
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:1764
get_xnsum2
subroutine get_xnsum2(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, zslm, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
Definition: mtnlm7_oclsm.f:4593
mnmxja
subroutine mnmxja(IM, JM, A, imax, jmax, title)
Print out the maximum and minimum values of an array.
Definition: mtnlm7_oclsm.f:4046
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:4486
maxmin
subroutine maxmin(ia, len, tile)
Print the maximum, mininum, mean and standard deviation of an array.
Definition: mtnlm7_oclsm.f:4209
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:2566
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:3470
read_g
subroutine read_g(glob, ITOPO)
Read input global 30-arc second orography data.
Definition: mtnlm7_oclsm.f:4158
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:3389
revers
subroutine revers(IM, JM, numi, F, WRK)
Reverse the east-west and north-south axes in a two-dimensional array.
Definition: mtnlm7_oclsm.f:3922
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:1598
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:2894
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:2960
makepc2
subroutine makepc2(ZAVG, ZSLM, THETA, GAMMA, SIGMA, GLAT, IM, JM, IMN, JMN, lon_c, lat_c)
Make the principle coordinates - slope of orography, anisotropy, angle of mountain range with respect...
Definition: mtnlm7_oclsm.f:2292