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