orog_mask_tools  1.8.0
mtnlm7_oclsm.f
Go to the documentation of this file.
1 C> @file
2 C> Terrain maker for global spectral model.
3 C> @author Mark Iredell @date 92-04-16
4 
5 C> This program creates 7 terrain-related files computed from the 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
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,ii,jj,i1,numx,i2
1900  integer ilist(IMN)
1901  real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1,XW2,XW4
1902 !jaa
1903  real :: xnsum_j,xland_j,xwatr_j
1904  logical inside_a_polygon
1905 C==== DATA DEBUG/.TRUE./
1906  DATA debug/.false./
1907 C
1908 ! ---- OCLSM holds the ocean (im,jm) grid
1909 ! --- mskocn=1 Use ocean model sea land mask, OK and present,
1910 ! --- mskocn=0 dont use Ocean model sea land mask, not OK, not present
1911  print *,' _____ SUBROUTINE MAKEMT2 '
1912  allocate(hgt_1d(maxsum))
1913 C---- GLOBAL XLAT AND XLON ( DEGREE )
1914 C
1915  jm1 = jm - 1
1916  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1917 C
1918  DO j=1,jmn
1919  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1920  ENDDO
1921  DO i=1,imn
1922  glon(i) = 0. + (i-1) * delxn + delxn * 0.5
1923  ENDDO
1924 
1925  land_frac(:,:) = 0.0
1926 C
1927 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1928 C
1929 C (*j*) for hard wired zero offset (lambda s =0) for terr05
1930 !$omp parallel do
1931 !$omp* private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,xw2,xw4,lono,
1932 !$omp* lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,height,
1933 !$omp* hgt_1d)
1934  DO j=1,jm
1935 ! print*, "J=", J
1936  DO i=1,im
1937  oro(i,j) = 0.0
1938  var(i,j) = 0.0
1939  var4(i,j) = 0.0
1940  xnsum = 0.0
1941  xland = 0.0
1942  xwatr = 0.0
1943  nsum = 0
1944  xl1 = 0.0
1945  xs1 = 0.0
1946  xw1 = 0.0
1947  xw2 = 0.0
1948  xw4 = 0.0
1949 
1950  lono(1) = lon_c(i,j)
1951  lono(2) = lon_c(i+1,j)
1952  lono(3) = lon_c(i+1,j+1)
1953  lono(4) = lon_c(i,j+1)
1954  lato(1) = lat_c(i,j)
1955  lato(2) = lat_c(i+1,j)
1956  lato(3) = lat_c(i+1,j+1)
1957  lato(4) = lat_c(i,j+1)
1958  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
1959  do jj = jst, jen; do i2 = 1, numx
1960  ii = ilist(i2)
1961  loni = ii*delxn
1962  lati = -90 + jj*delxn
1963  if(inside_a_polygon(loni*d2r,lati*d2r,4,
1964  & lono*d2r,lato*d2r))then
1965 
1966  xland = xland + float(zslm(ii,jj))
1967  xwatr = xwatr + float(1-zslm(ii,jj))
1968  xnsum = xnsum + 1.
1969  height = float(zavg(ii,jj))
1970  nsum = nsum+1
1971  if(nsum > maxsum) then
1972  print*, "nsum is greater than MAXSUM, increase MAXSUM"
1973  call abort()
1974  endif
1975  hgt_1d(nsum) = height
1976  IF(height.LT.-990.) height = 0.0
1977  xl1 = xl1 + height * float(zslm(ii,jj))
1978  xs1 = xs1 + height * float(1-zslm(ii,jj))
1979  xw1 = xw1 + height
1980  xw2 = xw2 + height ** 2
1981  endif
1982  enddo ; enddo
1983 
1984 
1985  IF(xnsum.GT.1.) THEN
1986 ! --- SLM initialized with OCLSM calc from all land points except ....
1987 ! --- 0 is ocean and 1 is land for slm
1988 ! --- Step 1 is to only change SLM after GFS SLM is applied
1989  land_frac(i,j) = xland/xnsum
1990  slm(i,j) = float(nint(xland/xnsum))
1991  IF(slm(i,j).NE.0.) THEN
1992  oro(i,j)= xl1 / xland
1993  ELSE
1994  oro(i,j)= xs1 / xwatr
1995  ENDIF
1996  var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
1997  do i1 = 1, nsum
1998  xw4 = xw4 + (hgt_1d(i1) - oro(i,j)) ** 4
1999  enddo
2000 
2001  IF(var(i,j).GT.1.) THEN
2002  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
2003  ENDIF
2004  ENDIF
2005  ENDDO
2006  ENDDO
2007 !$omp end parallel do
2008  WRITE(6,*) "! MAKEMT2 ORO SLM VAR VAR4 DONE"
2009 C
2010  deallocate(hgt_1d)
2011  RETURN
2012  END
2013 
2042  SUBROUTINE makepc(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2043  1 GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
2044 C
2045 C=== PC: principal coordinates of each Z avg orog box for L&M
2046 C
2047  parameter(rearth=6.3712e+6)
2048  dimension glat(jmn),xlat(jm),deltax(jmn)
2049  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2050  DIMENSION ORO(IM,JM),SLM(IM,JM),HL(IM,JM),HK(IM,JM)
2051  DIMENSION HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
2052  dimension theta(im,jm),gamma(im,jm),sigma2(im,jm),sigma(im,jm)
2053  dimension ist(im,jm),ien(im,jm),jst(jm),jen(jm),numi(jm)
2054  LOGICAL FLAG, DEBUG
2055 C=== DATA DEBUG/.TRUE./
2056  DATA debug/.false./
2057 C
2058  pi = 4.0 * atan(1.0)
2059  certh = pi * rearth
2060 C---- GLOBAL XLAT AND XLON ( DEGREE )
2061 C
2062  jm1 = jm - 1
2063  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2064  deltay = certh / float(jmn)
2065  print *, 'MAKEPC: DELTAY=',deltay
2066 C
2067  DO j=1,jmn
2068  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2069  deltax(j) = deltay * cos(glat(j)*pi/180.0)
2070  ENDDO
2071 C
2072 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2073 C
2074  DO j=1,jm
2075  DO i=1,numi(j)
2076 C IM1 = numi(j) - 1
2077  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
2078  faclon = delx / delxn
2079  ist(i,j) = faclon * float(i-1) - faclon * 0.5
2080  ist(i,j) = ist(i,j) + 1
2081  ien(i,j) = faclon * float(i) - faclon * 0.5
2082 C if (debug) then
2083 C if ( I .lt. 10 .and. J .lt. 10 )
2084 C 1 PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
2085 C endif
2086 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
2087 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
2088  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2089  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2090  if (debug) then
2091  if ( i .lt. 10 .and. j .lt. 10 )
2092  1 print*, ' I j IST IEN ',i,j,ist(i,j),ien(i,j)
2093  endif
2094  IF (ien(i,j) .LT. ist(i,j))
2095  1 print *,' MAKEPC: IEN < IST: I,J,IST(I,J),IEN(I,J)',
2096  2 i,j,ist(i,j),ien(i,j)
2097  ENDDO
2098  ENDDO
2099  DO j=1,jm-1
2100  flag=.true.
2101  DO j1=1,jmn
2102  xxlat = (xlat(j)+xlat(j+1))/2.
2103  IF(flag.AND.glat(j1).GT.xxlat) THEN
2104  jst(j) = j1
2105  jen(j+1) = j1 - 1
2106  flag = .false.
2107  ENDIF
2108  ENDDO
2109  ENDDO
2110  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2111  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2112  if (debug) then
2113  print*, ' IST,IEN(1,1-numi(1,JM))',ist(1,1),ien(1,1),
2114  1 ist(numi(jm),jm),ien(numi(jm),jm), numi(jm)
2115  print*, ' JST,JEN(1,JM) ',jst(1),jen(1),jst(jm),jen(jm)
2116  endif
2117 C
2118 C... DERIVITIVE TENSOR OF HEIGHT
2119 C
2120  DO j=1,jm
2121  DO i=1,numi(j)
2122  oro(i,j) = 0.0
2123  hx2(i,j) = 0.0
2124  hy2(i,j) = 0.0
2125  hxy(i,j) = 0.0
2126  xnsum = 0.0
2127  xland = 0.0
2128  xwatr = 0.0
2129  xl1 = 0.0
2130  xs1 = 0.0
2131  xfp = 0.0
2132  yfp = 0.0
2133  xfpyfp = 0.0
2134  xfp2 = 0.0
2135  yfp2 = 0.0
2136  hl(i,j) = 0.0
2137  hk(i,j) = 0.0
2138  hlprim(i,j) = 0.0
2139  theta(i,j) = 0.0
2140  gamma(i,j) = 0.
2141  sigma2(i,j) = 0.
2142  sigma(i,j) = 0.
2143 C
2144  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2145  i1 = ist(i,j) + ii1 - 1
2146  IF(i1.LE.0.) i1 = i1 + imn
2147  IF(i1.GT.imn) i1 = i1 - imn
2148 C
2149 C=== set the rest of the indexs for ave: 2pt staggered derivitive
2150 C
2151  i0 = i1 - 1
2152  if (i1 - 1 .le. 0 ) i0 = i0 + imn
2153  if (i1 - 1 .gt. imn) i0 = i0 - imn
2154 C
2155  ip1 = i1 + 1
2156  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2157  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2158 C
2159  DO j1=jst(j),jen(j)
2160  if (debug) then
2161  if ( i1 .eq. ist(i,j) .and. j1 .eq. jst(j) )
2162  1 print*, ' J, J1,IST,JST,DELTAX,GLAT ',
2163  2 j,j1,ist(i,j),jst(j),deltax(j1),glat(j1)
2164  if ( i1 .eq. ien(i,j) .and. j1 .eq. jen(j) )
2165  1 print*, ' J, J1,IEN,JEN,DELTAX,GLAT ',
2166  2 j,j1,ien(i,j),jen(j),deltax(j1),glat(j1)
2167  endif
2168  xland = xland + float(zslm(i1,j1))
2169  xwatr = xwatr + float(1-zslm(i1,j1))
2170  xnsum = xnsum + 1.
2171 C
2172  height = float(zavg(i1,j1))
2173  hi0 = float(zavg(i0,j1))
2174  hip1 = float(zavg(ip1,j1))
2175 C
2176  IF(height.LT.-990.) height = 0.0
2177  if(hi0 .lt. -990.) hi0 = 0.0
2178  if(hip1 .lt. -990.) hip1 = 0.0
2179 C........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
2180  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2181  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2182 C
2183 ! --- not at boundaries
2184 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
2185  if ( j1 .ne. jst(jm) .and. j1 .ne. jen(1) ) then
2186  hj0 = float(zavg(i1,j1-1))
2187  hjp1 = float(zavg(i1,j1+1))
2188  if(hj0 .lt. -990.) hj0 = 0.0
2189  if(hjp1 .lt. -990.) hjp1 = 0.0
2190 C....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
2191  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2192  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2193 C
2194 C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
2195 C === the NH pole: NB J1 goes from High at NP to Low toward SP
2196 C
2197 !RAB elseif ( J1 .eq. JST(1) ) then
2198  elseif ( j1 .eq. jst(jm) ) then
2199  ijax = i1 + imn/2
2200  if (ijax .le. 0 ) ijax = ijax + imn
2201  if (ijax .gt. imn) ijax = ijax - imn
2202 C..... at N pole we stay at the same latitude j1 but cross to opp side
2203  hijax = float(zavg(ijax,j1))
2204  hi1j1 = float(zavg(i1,j1))
2205  if(hijax .lt. -990.) hijax = 0.0
2206  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2207 C....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
2208  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2209  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2210  1 / deltay )**2
2211 C
2212 C === the SH pole: NB J1 goes from High at NP to Low toward SP
2213 C
2214 !RAB elseif ( J1 .eq. JEN(JM) ) then
2215  elseif ( j1 .eq. jen(1) ) then
2216  ijax = i1 + imn/2
2217  if (ijax .le. 0 ) ijax = ijax + imn
2218  if (ijax .gt. imn) ijax = ijax - imn
2219  hijax = float(zavg(ijax,j1))
2220  hi1j1 = float(zavg(i1,j1))
2221  if(hijax .lt. -990.) hijax = 0.0
2222  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2223  if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,hijax,hi1j1
2224 C..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
2225  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2226  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2227  1 / deltay )**2
2228  endif
2229 C
2230 C === The above does an average across the pole for the bndry in j.
2231 C23456789012345678901234567890123456789012345678901234567890123456789012......
2232 C
2233  xfpyfp = xfpyfp + xfp * yfp
2234  xl1 = xl1 + height * float(zslm(i1,j1))
2235  xs1 = xs1 + height * float(1-zslm(i1,j1))
2236 C
2237 C === average the HX2, HY2 and HXY
2238 C === This will be done over all land
2239 C
2240  ENDDO
2241  ENDDO
2242 C
2243 C === HTENSR
2244 C
2245  IF(xnsum.GT.1.) THEN
2246  slm(i,j) = float(nint(xland/xnsum))
2247  IF(slm(i,j).NE.0.) THEN
2248  oro(i,j)= xl1 / xland
2249  hx2(i,j) = xfp2 / xland
2250  hy2(i,j) = yfp2 / xland
2251  hxy(i,j) = xfpyfp / xland
2252  ELSE
2253  oro(i,j)= xs1 / xwatr
2254  ENDIF
2255 C=== degub testing
2256  if (debug) then
2257  print *," I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2258  1 xland,slm(i,j)
2259  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2260  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2261  ENDIF
2262 C
2263 C === make the principal axes, theta, and the degree of anisotropy,
2264 C === and sigma2, the slope parameter
2265 C
2266  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2267  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2268  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2269  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
2270 C
2271  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) * 180.0 / pi
2272 C === for testing print out in degrees
2273 C THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
2274  ENDIF
2275  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2276  if ( sigma2(i,j) .GE. 0. ) then
2277  sigma(i,j) = sqrt(sigma2(i,j) )
2278  if (sigma2(i,j) .ne. 0. .and.
2279  & hk(i,j) .GE. hlprim(i,j) )
2280  1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2281  else
2282  sigma(i,j)=0.
2283  endif
2284  ENDIF
2285  if (debug) then
2286  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2287  1 sigma(i,j),gamma(i,j)
2288  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2289  endif
2290  ENDDO
2291  ENDDO
2292  WRITE(6,*) "! MAKE Principal Coord DONE"
2293 C
2294  RETURN
2295  END
2296 
2316  SUBROUTINE makepc2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
2317  1 GLAT,IM,JM,IMN,JMN,lon_c,lat_c)
2318 C
2319 C=== PC: principal coordinates of each Z avg orog box for L&M
2320 C
2321  implicit none
2322  real, parameter :: REARTH=6.3712e+6
2323  real, parameter :: D2R = 3.14159265358979/180.
2324  integer :: im,jm,imn,jmn
2325  real :: GLAT(JMN),DELTAX(JMN)
2326  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2327  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
2328  real ORO(IM,JM),SLM(IM,JM),HL(IM,JM),HK(IM,JM)
2329  real HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
2330  real THETA(IM,JM),GAMMA(IM,JM),SIGMA2(IM,JM),SIGMA(IM,JM)
2331  real PI,CERTH,DELXN,DELTAY,XNSUM,XLAND,XWATR,XL1,XS1
2332  real xfp,yfp,xfpyfp,xfp2,yfp2,HEIGHT
2333  real hi0,hip1,hj0,hjp1,hijax,hi1j1
2334  real LONO(4),LATO(4),LONI,LATI
2335  integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
2336  integer ilist(IMN)
2337  logical inside_a_polygon
2338  LOGICAL FLAG, DEBUG
2339 C=== DATA DEBUG/.TRUE./
2340  DATA debug/.false./
2341 C
2342  pi = 4.0 * atan(1.0)
2343  certh = pi * rearth
2344 C---- GLOBAL XLAT AND XLON ( DEGREE )
2345 C
2346  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2347  deltay = certh / float(jmn)
2348  print *, 'MAKEPC2: DELTAY=',deltay
2349 C
2350  DO j=1,jmn
2351  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2352  deltax(j) = deltay * cos(glat(j)*d2r)
2353  ENDDO
2354 C
2355 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2356 C
2357 
2358 C... DERIVITIVE TENSOR OF HEIGHT
2359 C
2360 !$omp parallel do
2361 !$omp* private (j,i,xnsum,xland,xwatr,xl1,xs1,xfp,yfp,xfpyfp,
2362 !$omp* xfp2,yfp2,lono,lato,jst,jen,ilist,numx,j1,i2,i1,
2363 !$omp* loni,lati,i0,ip1,height,hi0,hip1,hj0,hjp1,ijax,
2364 !$omp* hijax,hi1j1)
2365  DO j=1,jm
2366 ! print*, "J=", J
2367  DO i=1,im
2368  oro(i,j) = 0.0
2369  hx2(i,j) = 0.0
2370  hy2(i,j) = 0.0
2371  hxy(i,j) = 0.0
2372  xnsum = 0.0
2373  xland = 0.0
2374  xwatr = 0.0
2375  xl1 = 0.0
2376  xs1 = 0.0
2377  xfp = 0.0
2378  yfp = 0.0
2379  xfpyfp = 0.0
2380  xfp2 = 0.0
2381  yfp2 = 0.0
2382  hl(i,j) = 0.0
2383  hk(i,j) = 0.0
2384  hlprim(i,j) = 0.0
2385  theta(i,j) = 0.0
2386  gamma(i,j) = 0.
2387  sigma2(i,j) = 0.
2388  sigma(i,j) = 0.
2389 
2390  lono(1) = lon_c(i,j)
2391  lono(2) = lon_c(i+1,j)
2392  lono(3) = lon_c(i+1,j+1)
2393  lono(4) = lon_c(i,j+1)
2394  lato(1) = lat_c(i,j)
2395  lato(2) = lat_c(i+1,j)
2396  lato(3) = lat_c(i+1,j+1)
2397  lato(4) = lat_c(i,j+1)
2398  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
2399 
2400  do j1 = jst, jen; do i2 = 1, numx
2401  i1 = ilist(i2)
2402  loni = i1*delxn
2403  lati = -90 + j1*delxn
2404  if(inside_a_polygon(loni*d2r,lati*d2r,4,
2405  & lono*d2r,lato*d2r))then
2406 
2407 C=== set the rest of the indexs for ave: 2pt staggered derivitive
2408 C
2409  i0 = i1 - 1
2410  if (i1 - 1 .le. 0 ) i0 = i0 + imn
2411  if (i1 - 1 .gt. imn) i0 = i0 - imn
2412 C
2413  ip1 = i1 + 1
2414  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
2415  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
2416 
2417  xland = xland + float(zslm(i1,j1))
2418  xwatr = xwatr + float(1-zslm(i1,j1))
2419  xnsum = xnsum + 1.
2420 C
2421  height = float(zavg(i1,j1))
2422  hi0 = float(zavg(i0,j1))
2423  hip1 = float(zavg(ip1,j1))
2424 C
2425  IF(height.LT.-990.) height = 0.0
2426  if(hi0 .lt. -990.) hi0 = 0.0
2427  if(hip1 .lt. -990.) hip1 = 0.0
2428 C........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
2429  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
2430  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
2431 C
2432 ! --- not at boundaries
2433 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
2434  if ( j1 .ne. 1 .and. j1 .ne. jmn ) then
2435  hj0 = float(zavg(i1,j1-1))
2436  hjp1 = float(zavg(i1,j1+1))
2437  if(hj0 .lt. -990.) hj0 = 0.0
2438  if(hjp1 .lt. -990.) hjp1 = 0.0
2439 C....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
2440  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
2441  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
2442 C
2443 C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
2444 C === the NH pole: NB J1 goes from High at NP to Low toward SP
2445 C
2446 !RAB elseif ( J1 .eq. JST(1) ) then
2447  elseif ( j1 .eq. 1 ) then
2448  ijax = i1 + imn/2
2449  if (ijax .le. 0 ) ijax = ijax + imn
2450  if (ijax .gt. imn) ijax = ijax - imn
2451 C..... at N pole we stay at the same latitude j1 but cross to opp side
2452  hijax = float(zavg(ijax,j1))
2453  hi1j1 = float(zavg(i1,j1))
2454  if(hijax .lt. -990.) hijax = 0.0
2455  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2456 C....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
2457  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
2458  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) )
2459  1 / deltay )**2
2460 C
2461 C === the SH pole: NB J1 goes from High at NP to Low toward SP
2462 C
2463 !RAB elseif ( J1 .eq. JEN(JM) ) then
2464  elseif ( j1 .eq. jmn ) then
2465  ijax = i1 + imn/2
2466  if (ijax .le. 0 ) ijax = ijax + imn
2467  if (ijax .gt. imn) ijax = ijax - imn
2468  hijax = float(zavg(ijax,j1))
2469  hi1j1 = float(zavg(i1,j1))
2470  if(hijax .lt. -990.) hijax = 0.0
2471  if(hi1j1 .lt. -990.) hi1j1 = 0.0
2472  if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,
2473  & hijax,hi1j1
2474 C..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
2475  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
2476  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) )
2477  1 / deltay )**2
2478  endif
2479 C
2480 C === The above does an average across the pole for the bndry in j.
2481 C23456789012345678901234567890123456789012345678901234567890123456789012......
2482 C
2483  xfpyfp = xfpyfp + xfp * yfp
2484  xl1 = xl1 + height * float(zslm(i1,j1))
2485  xs1 = xs1 + height * float(1-zslm(i1,j1))
2486  ENDIF
2487 C
2488 C === average the HX2, HY2 and HXY
2489 C === This will be done over all land
2490 C
2491  ENDDO
2492  ENDDO
2493 C
2494 C === HTENSR
2495 C
2496  IF(xnsum.GT.1.) THEN
2497  slm(i,j) = float(nint(xland/xnsum))
2498  IF(slm(i,j).NE.0.) THEN
2499  oro(i,j)= xl1 / xland
2500  hx2(i,j) = xfp2 / xland
2501  hy2(i,j) = yfp2 / xland
2502  hxy(i,j) = xfpyfp / xland
2503  ELSE
2504  oro(i,j)= xs1 / xwatr
2505  ENDIF
2506 C=== degub testing
2507  if (debug) then
2508  print *," I,J,i1,j1,HEIGHT:", i,j,i1,j1,height,
2509  1 xland,slm(i,j)
2510  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
2511  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
2512  ENDIF
2513 C
2514 C === make the principal axes, theta, and the degree of anisotropy,
2515 C === and sigma2, the slope parameter
2516 C
2517  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
2518  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
2519  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
2520  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
2521 C
2522  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) / d2r
2523 C === for testing print out in degrees
2524 C THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
2525  ENDIF
2526  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
2527  if ( sigma2(i,j) .GE. 0. ) then
2528  sigma(i,j) = sqrt(sigma2(i,j) )
2529  if (sigma2(i,j) .ne. 0. .and.
2530  & hk(i,j) .GE. hlprim(i,j) )
2531  1 gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
2532  else
2533  sigma(i,j)=0.
2534  endif
2535  ENDIF
2536  if (debug) then
2537  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),
2538  1 sigma(i,j),gamma(i,j)
2539  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
2540  endif
2541  ENDDO
2542  ENDDO
2543 !$omp end parallel do
2544  WRITE(6,*) "! MAKE Principal Coord DONE"
2545 C
2546  RETURN
2547  END SUBROUTINE makepc2
2548 
2590  SUBROUTINE makeoa(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2591  1 ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
2592  2 IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
2593  dimension glat(jmn),xlat(jm)
2594  INTEGER ZAVG(IMN,JMN)
2595  DIMENSION ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
2596  DIMENSION OA4(IM,JM,4),IOA4(IM,JM,4)
2597  DIMENSION IST(IM,jm),IEN(IM,jm),JST(JM),JEN(JM)
2598  dimension xnsum(im,jm),xnsum1(im,jm),xnsum2(im,jm)
2599  dimension xnsum3(im,jm),xnsum4(im,jm)
2600  dimension var(im,jm),ol(im,jm,4),numi(jm)
2601  LOGICAL FLAG
2602 C
2603 C---- GLOBAL XLAT AND XLON ( DEGREE )
2604 C
2605 ! --- IM1 = IM - 1 removed (not used in this sub)
2606  jm1 = jm - 1
2607  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
2608 C
2609  DO j=1,jmn
2610  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
2611  ENDDO
2612  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
2613 C
2614 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
2615 C
2616  DO j=1,jm
2617  DO i=1,numi(j)
2618  delx = 360./numi(j) ! GAUSSIAN GRID RESOLUTION
2619  faclon = delx / delxn
2620 C --- minus sign here in IST and IEN as in MAKEMT!
2621  ist(i,j) = faclon * float(i-1) - faclon * 0.5
2622  ist(i,j) = ist(i,j) + 1
2623  ien(i,j) = faclon * float(i) - faclon * 0.5
2624 ! IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
2625 ! IEN(I,j) = FACLON * FLOAT(I) + 0.0001
2626  IF (ist(i,j) .LE. 0) ist(i,j) = ist(i,j) + imn
2627  IF (ien(i,j) .LT. ist(i,j)) ien(i,j) = ien(i,j) + imn
2628 cx PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
2629  if ( i .lt. 3 .and. j .lt. 3 )
2630  1print*,' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2631  if ( i .lt. 3 .and. j .ge. jm-1 )
2632  1print*,' MAKEOA: I j IST IEN ',i,j,ist(i,j),ien(i,j)
2633  ENDDO
2634  ENDDO
2635  print *,'MAKEOA: DELXN,DELX,FACLON',delxn,delx,faclon
2636  print *, ' ***** ready to start JST JEN section '
2637  DO j=1,jm-1
2638  flag=.true.
2639  DO j1=1,jmn
2640 ! --- XXLAT added as in MAKEMT and in next line as well
2641  xxlat = (xlat(j)+xlat(j+1))/2.
2642  IF(flag.AND.glat(j1).GT.xxlat) THEN
2643  jst(j) = j1
2644 ! --- JEN(J+1) = J1 - 1
2645  flag = .false.
2646  if ( j .eq. 1 )
2647  1print*,' MAKEOA: XX j JST JEN ',j,jst(j),jen(j)
2648  ENDIF
2649  ENDDO
2650  if ( j .lt. 3 )
2651  1print*,' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2652  if ( j .ge. jm-2 )
2653  1print*,' MAKEOA: j JST JEN ',j,jst(j),jen(j)
2654 C FLAG=.TRUE.
2655 C DO J1=JST(J),JMN
2656 C IF(FLAG.AND.GLAT(J1).GT.XLAT(J)) THEN
2657 C JEN(J) = J1 - 1
2658 C FLAG = .FALSE.
2659 C ENDIF
2660 C ENDDO
2661  ENDDO
2662  jst(jm) = max(jst(jm-1) - (jen(jm-1)-jst(jm-1)),1)
2663  jen(1) = min(jen(2) + (jen(2)-jst(2)),jmn)
2664  print *,' ***** JST(1) JEN(1) ',jst(1),jen(1)
2665  print *,' ***** JST(JM) JEN(JM) ',jst(jm),jen(jm)
2666 C
2667  DO j=1,jm
2668  DO i=1,numi(j)
2669  xnsum(i,j) = 0.0
2670  elvmax(i,j) = oro(i,j)
2671  zmax(i,j) = 0.0
2672  ENDDO
2673  ENDDO
2674 !
2675 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
2676 ! --- to JM or to JM1
2677  DO j=1,jm
2678  DO i=1,numi(j)
2679  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2680  i1 = ist(i,j) + ii1 - 1
2681 ! --- next line as in makemt (I1 not II1) (*j*) 20070701
2682  IF(i1.LE.0.) i1 = i1 + imn
2683  IF (i1 .GT. imn) i1 = i1 - imn
2684  DO j1=jst(j),jen(j)
2685  height = float(zavg(i1,j1))
2686  IF(height.LT.-990.) height = 0.0
2687  IF ( height .gt. oro(i,j) ) then
2688  if ( height .gt. zmax(i,j) )zmax(i,j) = height
2689  xnsum(i,j) = xnsum(i,j) + 1
2690  ENDIF
2691  ENDDO
2692  ENDDO
2693  if ( i .lt. 5 .and. j .ge. jm-5 ) then
2694  print *,' I,J,ORO(I,J),XNSUM(I,J),ZMAX(I,J):',
2695  1 i,j,oro(i,j),xnsum(i,j),zmax(i,j)
2696  endif
2697  ENDDO
2698  ENDDO
2699 !
2700 C.... make ELVMAX ORO from MAKEMT sub
2701 C
2702 ! --- this will make work1 array take on oro's values on return
2703  DO j=1,jm
2704  DO i=1,numi(j)
2705 
2706  oro1(i,j) = oro(i,j)
2707  elvmax(i,j) = zmax(i,j)
2708  ENDDO
2709  ENDDO
2710 C........
2711 C The MAX elev peak (no averaging)
2712 C........
2713 ! DO J=1,JM
2714 ! DO I=1,numi(j)
2715 ! DO II1 = 1, IEN(I,J) - IST(I,J) + 1
2716 ! I1 = IST(I,J) + II1 - 1
2717 ! IF(I1.LE.0.) I1 = I1 + IMN
2718 ! IF(I1.GT.IMN) I1 = I1 - IMN
2719 ! DO J1=JST(J),JEN(J)
2720 ! if ( ELVMAX(I,J) .lt. ZMAX(I1,J1))
2721 ! 1 ELVMAX(I,J) = ZMAX(I1,J1)
2722 ! ENDDO
2723 ! ENDDO
2724 ! ENDDO
2725 ! ENDDO
2726 C
2727 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
2728 C IN A GRID BOX
2729  DO j=1,jm
2730  DO i=1,numi(j)
2731  xnsum1(i,j) = 0.0
2732  xnsum2(i,j) = 0.0
2733  xnsum3(i,j) = 0.0
2734  xnsum4(i,j) = 0.0
2735  ENDDO
2736  ENDDO
2737 ! --- loop
2738  DO j=1,jm1
2739  DO i=1,numi(j)
2740  hc = 1116.2 - 0.878 * var(i,j)
2741 ! print *,' I,J,HC,VAR:',I,J,HC,VAR(I,J)
2742  DO ii1 = 1, ien(i,j) - ist(i,j) + 1
2743  i1 = ist(i,j) + ii1 - 1
2744 ! IF (I1.LE.0.) print *,' I1 less than 0',I1,II1,IST(I,J),IEN(I,J)
2745 ! if ( J .lt. 3 .or. J .gt. JM-2 ) then
2746 ! IF(I1 .GT. IMN)print *,' I1 > IMN',J,I1,II1,IMN,IST(I,J),IEN(I,J)
2747 ! endif
2748  IF(i1.GT.imn) i1 = i1 - imn
2749  DO j1=jst(j),jen(j)
2750  IF(float(zavg(i1,j1)) .GT. hc)
2751  1 xnsum1(i,j) = xnsum1(i,j) + 1
2752  xnsum2(i,j) = xnsum2(i,j) + 1
2753  ENDDO
2754  ENDDO
2755 C
2756  inci = nint((ien(i,j)-ist(i,j)) * 0.5)
2757  isttt = min(max(ist(i,j)-inci,1),imn)
2758  ieddd = min(max(ien(i,j)-inci,1),imn)
2759 C
2760  incj = nint((jen(j)-jst(j)) * 0.5)
2761  jsttt = min(max(jst(j)-incj,1),jmn)
2762  jeddd = min(max(jen(j)-incj,1),jmn)
2763 ! if ( J .lt. 3 .or. J .gt. JM-3 ) then
2764 ! if(I .lt. 3 .or. I .gt. IM-3) then
2765 ! print *,' INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD:',
2766 ! 1 I,J,INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD
2767 ! endif
2768 ! endif
2769 C
2770  DO i1=isttt,ieddd
2771  DO j1=jsttt,jeddd
2772  IF(float(zavg(i1,j1)) .GT. hc)
2773  1 xnsum3(i,j) = xnsum3(i,j) + 1
2774  xnsum4(i,j) = xnsum4(i,j) + 1
2775  ENDDO
2776  ENDDO
2777 cx print*,' i j hc var ',i,j,hc,var(i,j)
2778 cx print*,'xnsum12 ',xnsum1(i,j),xnsum2(i,j)
2779 cx print*,'xnsum34 ',xnsum3(i,j),xnsum4(i,j)
2780  ENDDO
2781  ENDDO
2782 C
2783 C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
2784 C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
2785 C (KWD = 1 2 3 4)
2786 C ( WD = W S SW NW)
2787 C
2788 C
2789  DO kwd = 1, 4
2790  DO j=1,jm
2791  DO i=1,numi(j)
2792  oa4(i,j,kwd) = 0.0
2793  ENDDO
2794  ENDDO
2795  ENDDO
2796 C
2797  DO j=1,jm-2
2798  DO i=1,numi(j)
2799  ii = i + 1
2800  IF (ii .GT. numi(j)) ii = ii - numi(j)
2801  xnpu = xnsum(i,j) + xnsum(i,j+1)
2802  xnpd = xnsum(ii,j) + xnsum(ii,j+1)
2803  IF (xnpd .NE. xnpu) oa4(ii,j+1,1) = 1. - xnpd / max(xnpu , 1.)
2804  ol(ii,j+1,1) = (xnsum3(i,j+1)+xnsum3(ii,j+1))/
2805  1 (xnsum4(i,j+1)+xnsum4(ii,j+1))
2806 ! if ( I .lt. 20 .and. J .ge. JM-19 ) then
2807 ! PRINT*,' MAKEOA: I J IST IEN ',I,j,IST(I,J),IEN(I,J)
2808 ! PRINT*,' HC VAR ',HC,VAR(i,j)
2809 ! PRINT*,' MAKEOA: XNSUM(I,J)=',XNSUM(I,J),XNPU, XNPD
2810 ! PRINT*,' MAKEOA: XNSUM3(I,J+1),XNSUM3(II,J+1)',
2811 ! 1 XNSUM3(I,J+1),XNSUM3(II,J+1)
2812 ! PRINT*,' MAKEOA: II, OA4(II,J+1,1), OL(II,J+1,1):',
2813 ! 1 II, OA4(II,J+1,1), OL(II,J+1,1)
2814 ! endif
2815  ENDDO
2816  ENDDO
2817  DO j=1,jm-2
2818  DO i=1,numi(j)
2819  ii = i + 1
2820  IF (ii .GT. numi(j)) ii = ii - numi(j)
2821  xnpu = xnsum(i,j+1) + xnsum(ii,j+1)
2822  xnpd = xnsum(i,j) + xnsum(ii,j)
2823  IF (xnpd .NE. xnpu) oa4(ii,j+1,2) = 1. - xnpd / max(xnpu , 1.)
2824  ol(ii,j+1,2) = (xnsum3(ii,j)+xnsum3(ii,j+1))/
2825  1 (xnsum4(ii,j)+xnsum4(ii,j+1))
2826  ENDDO
2827  ENDDO
2828  DO j=1,jm-2
2829  DO i=1,numi(j)
2830  ii = i + 1
2831  IF (ii .GT. numi(j)) ii = ii - numi(j)
2832  xnpu = xnsum(i,j+1) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2833  xnpd = xnsum(ii,j) + ( xnsum(i,j) + xnsum(ii,j+1) )*0.5
2834  IF (xnpd .NE. xnpu) oa4(ii,j+1,3) = 1. - xnpd / max(xnpu , 1.)
2835  ol(ii,j+1,3) = (xnsum1(ii,j)+xnsum1(i,j+1))/
2836  1 (xnsum2(ii,j)+xnsum2(i,j+1))
2837  ENDDO
2838  ENDDO
2839  DO j=1,jm-2
2840  DO i=1,numi(j)
2841  ii = i + 1
2842  IF (ii .GT. numi(j)) ii = ii - numi(j)
2843  xnpu = xnsum(i,j) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2844  xnpd = xnsum(ii,j+1) + ( xnsum(ii,j) + xnsum(i,j+1) )*0.5
2845  IF (xnpd .NE. xnpu) oa4(ii,j+1,4) = 1. - xnpd / max(xnpu , 1.)
2846  ol(ii,j+1,4) = (xnsum1(i,j)+xnsum1(ii,j+1))/
2847  1 (xnsum2(i,j)+xnsum2(ii,j+1))
2848  ENDDO
2849  ENDDO
2850 C
2851  DO kwd = 1, 4
2852  DO i=1,numi(j)
2853  ol(i,1,kwd) = ol(i,2,kwd)
2854  ol(i,jm,kwd) = ol(i,jm-1,kwd)
2855  ENDDO
2856  ENDDO
2857 C
2858  DO kwd=1,4
2859  DO j=1,jm
2860  DO i=1,numi(j)
2861  t = oa4(i,j,kwd)
2862  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
2863  ENDDO
2864  ENDDO
2865  ENDDO
2866 C
2867  ns0 = 0
2868  ns1 = 0
2869  ns2 = 0
2870  ns3 = 0
2871  ns4 = 0
2872  ns5 = 0
2873  ns6 = 0
2874  DO kwd=1,4
2875  DO j=1,jm
2876  DO i=1,numi(j)
2877  t = abs( oa4(i,j,kwd) )
2878  IF(t .EQ. 0.) THEN
2879  ioa4(i,j,kwd) = 0
2880  ns0 = ns0 + 1
2881  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
2882  ioa4(i,j,kwd) = 1
2883  ns1 = ns1 + 1
2884  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
2885  ioa4(i,j,kwd) = 2
2886  ns2 = ns2 + 1
2887  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
2888  ioa4(i,j,kwd) = 3
2889  ns3 = ns3 + 1
2890  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
2891  ioa4(i,j,kwd) = 4
2892  ns4 = ns4 + 1
2893  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
2894  ioa4(i,j,kwd) = 5
2895  ns5 = ns5 + 1
2896  ELSE IF(t .GT. 10000.) THEN
2897  ioa4(i,j,kwd) = 6
2898  ns6 = ns6 + 1
2899  ENDIF
2900  ENDDO
2901  ENDDO
2902  ENDDO
2903 C
2904  WRITE(6,*) "! MAKEOA EXIT"
2905 C
2906  RETURN
2907  END SUBROUTINE makeoa
2908 
2918  function get_lon_angle(dx,lat, DEGRAD)
2919  implicit none
2920  real dx, lat, degrad
2921 
2922  real get_lon_angle
2923  real, parameter :: radius = 6371200
2924 
2925  get_lon_angle = 2*asin( sin(dx/radius*0.5)/cos(lat) )*degrad
2926 
2927  end function get_lon_angle
2928 
2937  function get_lat_angle(dy, DEGRAD)
2938  implicit none
2939  real dy, degrad
2940 
2941  real get_lat_angle
2942  real, parameter :: radius = 6371200
2943 
2944  get_lat_angle = dy/radius*degrad
2945 
2946  end function get_lat_angle
2947 
2984  SUBROUTINE makeoa2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
2985  1 ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
2986  2 IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,dx,dy,
2987  3 is_south_pole,is_north_pole )
2988  implicit none
2989  real, parameter :: MISSING_VALUE = -9999.
2990  real, parameter :: d2r = 3.14159265358979/180.
2991  real, PARAMETER :: R2D=180./3.14159265358979
2992  integer im,jm,imn,jmn
2993  real GLAT(JMN)
2994  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
2995  real ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
2996  real OA4(IM,JM,4)
2997  integer IOA4(IM,JM,4)
2998  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
2999  real lon_t(IM,JM), lat_t(IM,JM)
3000  real dx(IM,JM), dy(IM,JM)
3001  logical is_south_pole(IM,JM), is_north_pole(IM,JM)
3002  real XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
3003  real XNSUM3(IM,JM),XNSUM4(IM,JM)
3004  real VAR(IM,JM),OL(IM,JM,4)
3005  LOGICAL FLAG
3006  integer i,j,ilist(IMN),numx,i1,j1,ii1
3007  integer KWD,II,npts
3008  real LONO(4),LATO(4),LONI,LATI
3009  real DELXN,HC,HEIGHT,XNPU,XNPD,T
3010  integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
3011  logical inside_a_polygon
3012  real lon,lat,dlon,dlat,dlat_old
3013  real lon1,lat1,lon2,lat2
3014  real xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
3015  real HC_11, HC_12, HC_21, HC_22
3016  real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
3017  real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
3018  real get_lon_angle, get_lat_angle, get_xnsum
3019  integer ist, ien, jst, jen
3020  real xland,xwatr,xl1,xs1,oroavg,slm
3021 C
3022 C---- GLOBAL XLAT AND XLON ( DEGREE )
3023 C
3024  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
3025 C
3026  DO j=1,jmn
3027  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3028  ENDDO
3029  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
3030 C
3031 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
3032 C
3033 C
3034  DO j=1,jm
3035  DO i=1,im
3036  xnsum(i,j) = 0.0
3037  elvmax(i,j) = oro(i,j)
3038  zmax(i,j) = 0.0
3039 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
3040 C IN A GRID BOX
3041  xnsum1(i,j) = 0.0
3042  xnsum2(i,j) = 0.0
3043  xnsum3(i,j) = 0.0
3044  xnsum4(i,j) = 0.0
3045  oro1(i,j) = oro(i,j)
3046  elvmax(i,j) = zmax(i,j)
3047  ENDDO
3048  ENDDO
3049 
3050 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3051 ! --- to JM or to JM1
3052 !$omp parallel do
3053 !$omp* private (j,i,hc,lono,lato,jst,jen,ilist,numx,j1,ii1,i1,loni,
3054 !$omp* lati,height)
3055  DO j=1,jm
3056 ! print*, "J=", J
3057  DO i=1,im
3058  hc = 1116.2 - 0.878 * var(i,j)
3059  lono(1) = lon_c(i,j)
3060  lono(2) = lon_c(i+1,j)
3061  lono(3) = lon_c(i+1,j+1)
3062  lono(4) = lon_c(i,j+1)
3063  lato(1) = lat_c(i,j)
3064  lato(2) = lat_c(i+1,j)
3065  lato(3) = lat_c(i+1,j+1)
3066  lato(4) = lat_c(i,j+1)
3067  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3068  do j1 = jst, jen; do ii1 = 1, numx
3069  i1 = ilist(ii1)
3070  loni = i1*delxn
3071  lati = -90 + j1*delxn
3072  if(inside_a_polygon(loni*d2r,lati*d2r,4,
3073  & lono*d2r,lato*d2r))then
3074 
3075  height = float(zavg(i1,j1))
3076  IF(height.LT.-990.) height = 0.0
3077  IF ( height .gt. oro(i,j) ) then
3078  if ( height .gt. zmax(i,j) )zmax(i,j) = height
3079  ENDIF
3080  endif
3081  ENDDO ; ENDDO
3082  ENDDO
3083  ENDDO
3084 !$omp end parallel do
3085 C
3086 ! --- this will make work1 array take on oro's values on return
3087 ! --- this will make work1 array take on oro's values on return
3088  DO j=1,jm
3089  DO i=1,im
3090 
3091  oro1(i,j) = oro(i,j)
3092  elvmax(i,j) = zmax(i,j)
3093  ENDDO
3094  ENDDO
3095 
3096  DO kwd = 1, 4
3097  DO j=1,jm
3098  DO i=1,im
3099  oa4(i,j,kwd) = 0.0
3100  ol(i,j,kwd) = 0.0
3101  ENDDO
3102  ENDDO
3103  ENDDO
3104  !
3105 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3106 C
3107 C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
3108 C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
3109 C (KWD = 1 2 3 4)
3110 C ( WD = W S SW NW)
3111 C
3112 C
3113 !$omp parallel do
3114 !$omp* private (j,i,lon,lat,kwd,dlon,dlat,lon1,lon2,lat1,lat2,
3115 !$omp* xnsum11,xnsum12,xnsum21,xnsum22,xnpu,xnpd,
3116 !$omp* xnsum1_11,xnsum2_11,hc_11, xnsum1_12,xnsum2_12,
3117 !$omp* hc_12,xnsum1_21,xnsum2_21,hc_21, xnsum1_22,
3118 !$omp* xnsum2_22,hc_22)
3119  DO j=1,jm
3120 ! print*, "j = ", j
3121  DO i=1,im
3122  lon = lon_t(i,j)
3123  lat = lat_t(i,j)
3124  !--- for around north pole, oa and ol are all 0
3125 
3126  if(is_north_pole(i,j)) then
3127  print*, "set oa1 = 0 and ol=0 at i,j=", i,j
3128  do kwd = 1, 4
3129  oa4(i,j,kwd) = 0.
3130  ol(i,j,kwd) = 0.
3131  enddo
3132  else if(is_south_pole(i,j)) then
3133  print*, "set oa1 = 0 and ol=1 at i,j=", i,j
3134  do kwd = 1, 4
3135  oa4(i,j,kwd) = 0.
3136  ol(i,j,kwd) = 1.
3137  enddo
3138  else
3139 
3140  !--- for each point, find a lat-lon grid box with same dx and dy as the cubic grid box
3141  dlon = get_lon_angle(dx(i,j), lat*d2r, r2d )
3142  dlat = get_lat_angle(dy(i,j), r2d)
3143  !--- adjust dlat if the points are close to pole.
3144  if( lat-dlat*0.5<-90.) then
3145  print*, "at i,j =", i,j, lat, dlat, lat-dlat*0.5
3146  print*, "ERROR: lat-dlat*0.5<-90."
3147  call errexit(4)
3148  endif
3149  if( lat+dlat*2 > 90.) then
3150  dlat_old = dlat
3151  dlat = (90-lat)*0.5
3152  print*, "at i,j=",i,j," adjust dlat from ",
3153  & dlat_old, " to ", dlat
3154  endif
3155  !--- lower left
3156  lon1 = lon-dlon*1.5
3157  lon2 = lon-dlon*0.5
3158  lat1 = lat-dlat*0.5
3159  lat2 = lat+dlat*0.5
3160 
3161  if(lat1<-90 .or. lat2>90) then
3162  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3163  endif
3164  xnsum11 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3165  & zavg,zslm,delxn)
3166 
3167  !--- upper left
3168  lon1 = lon-dlon*1.5
3169  lon2 = lon-dlon*0.5
3170  lat1 = lat+dlat*0.5
3171  lat2 = lat+dlat*1.5
3172  if(lat1<-90 .or. lat2>90) then
3173  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3174  endif
3175  xnsum12 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3176  & zavg,zslm,delxn)
3177 
3178  !--- lower right
3179  lon1 = lon-dlon*0.5
3180  lon2 = lon+dlon*0.5
3181  lat1 = lat-dlat*0.5
3182  lat2 = lat+dlat*0.5
3183  if(lat1<-90 .or. lat2>90) then
3184  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3185  endif
3186  xnsum21 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3187  & zavg,zslm,delxn)
3188 
3189  !--- upper right
3190  lon1 = lon-dlon*0.5
3191  lon2 = lon+dlon*0.5
3192  lat1 = lat+dlat*0.5
3193  lat2 = lat+dlat*1.5
3194  if(lat1<-90 .or. lat2>90) then
3195  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3196  endif
3197 
3198  xnsum22 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
3199  & zavg,zslm,delxn)
3200 
3201  xnpu = xnsum11 + xnsum12
3202  xnpd = xnsum21 + xnsum22
3203  IF (xnpd .NE. xnpu) oa4(i,j,1) = 1. - xnpd / max(xnpu , 1.)
3204 
3205  xnpu = xnsum11 + xnsum21
3206  xnpd = xnsum12 + xnsum22
3207  IF (xnpd .NE. xnpu) oa4(i,j,2) = 1. - xnpd / max(xnpu , 1.)
3208 
3209  xnpu = xnsum11 + (xnsum12+xnsum21)*0.5
3210  xnpd = xnsum22 + (xnsum12+xnsum21)*0.5
3211  IF (xnpd .NE. xnpu) oa4(i,j,3) = 1. - xnpd / max(xnpu , 1.)
3212 
3213  xnpu = xnsum12 + (xnsum11+xnsum22)*0.5
3214  xnpd = xnsum21 + (xnsum11+xnsum22)*0.5
3215  IF (xnpd .NE. xnpu) oa4(i,j,4) = 1. - xnpd / max(xnpu , 1.)
3216 
3217 
3218  !--- calculate OL3 and OL4
3219  !--- lower left
3220  lon1 = lon-dlon*1.5
3221  lon2 = lon-dlon*0.5
3222  lat1 = lat-dlat*0.5
3223  lat2 = lat+dlat*0.5
3224  if(lat1<-90 .or. lat2>90) then
3225  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3226  endif
3227  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3228  & zavg,zslm,delxn, xnsum1_11, xnsum2_11, hc_11)
3229 
3230  !--- upper left
3231  lon1 = lon-dlon*1.5
3232  lon2 = lon-dlon*0.5
3233  lat1 = lat+dlat*0.5
3234  lat2 = lat+dlat*1.5
3235  if(lat1<-90 .or. lat2>90) then
3236  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3237  endif
3238  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3239  & zavg,zslm,delxn, xnsum1_12, xnsum2_12, hc_12)
3240 
3241  !--- lower right
3242  lon1 = lon-dlon*0.5
3243  lon2 = lon+dlon*0.5
3244  lat1 = lat-dlat*0.5
3245  lat2 = lat+dlat*0.5
3246  if(lat1<-90 .or. lat2>90) then
3247  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3248  endif
3249  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3250  & zavg,zslm,delxn, xnsum1_21, xnsum2_21, hc_21)
3251 
3252  !--- upper right
3253  lon1 = lon-dlon*0.5
3254  lon2 = lon+dlon*0.5
3255  lat1 = lat+dlat*0.5
3256  lat2 = lat+dlat*1.5
3257  if(lat1<-90 .or. lat2>90) then
3258  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3259  endif
3260  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
3261  & zavg,zslm,delxn, xnsum1_22, xnsum2_22, hc_22)
3262 
3263  ol(i,j,3) = (xnsum1_22+xnsum1_11)/(xnsum2_22+xnsum2_11)
3264  ol(i,j,4) = (xnsum1_12+xnsum1_21)/(xnsum2_12+xnsum2_21)
3265 
3266  !--- calculate OL1 and OL2
3267  !--- lower left
3268  lon1 = lon-dlon*2.0
3269  lon2 = lon-dlon
3270  lat1 = lat
3271  lat2 = lat+dlat
3272  if(lat1<-90 .or. lat2>90) then
3273  print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
3274  endif
3275  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3276  & zavg,zslm,delxn, xnsum1_11, xnsum2_11, hc_11)
3277 
3278  !--- upper left
3279  lon1 = lon-dlon*2.0
3280  lon2 = lon-dlon
3281  lat1 = lat+dlat
3282  lat2 = lat+dlat*2.0
3283  if(lat1<-90 .or. lat2>90) then
3284  print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
3285  endif
3286 
3287  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3288  & zavg,zslm,delxn, xnsum1_12, xnsum2_12, hc_12)
3289 
3290  !--- lower right
3291  lon1 = lon-dlon
3292  lon2 = lon
3293  lat1 = lat
3294  lat2 = lat+dlat
3295  if(lat1<-90 .or. lat2>90) then
3296  print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
3297  endif
3298  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3299  & zavg,zslm,delxn, xnsum1_21, xnsum2_21, hc_21)
3300 
3301  !--- upper right
3302  lon1 = lon-dlon
3303  lon2 = lon
3304  lat1 = lat+dlat
3305  lat2 = lat+dlat*2.0
3306  if(lat1<-90 .or. lat2>90) then
3307  print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
3308  endif
3309 
3310  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
3311  & zavg,zslm,delxn, xnsum1_22, xnsum2_22, hc_22)
3312 
3313  ol(i,j,1) = (xnsum1_11+xnsum1_21)/(xnsum2_11+xnsum2_21)
3314  ol(i,j,2) = (xnsum1_21+xnsum1_22)/(xnsum2_21+xnsum2_22)
3315  ENDIF
3316  ENDDO
3317  ENDDO
3318 !$omp end parallel do
3319  DO kwd=1,4
3320  DO j=1,jm
3321  DO i=1,im
3322  t = oa4(i,j,kwd)
3323  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3324  ENDDO
3325  ENDDO
3326  ENDDO
3327 C
3328  ns0 = 0
3329  ns1 = 0
3330  ns2 = 0
3331  ns3 = 0
3332  ns4 = 0
3333  ns5 = 0
3334  ns6 = 0
3335  DO kwd=1,4
3336  DO j=1,jm
3337  DO i=1,im
3338  t = abs( oa4(i,j,kwd) )
3339  IF(t .EQ. 0.) THEN
3340  ioa4(i,j,kwd) = 0
3341  ns0 = ns0 + 1
3342  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
3343  ioa4(i,j,kwd) = 1
3344  ns1 = ns1 + 1
3345  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
3346  ioa4(i,j,kwd) = 2
3347  ns2 = ns2 + 1
3348  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
3349  ioa4(i,j,kwd) = 3
3350  ns3 = ns3 + 1
3351  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
3352  ioa4(i,j,kwd) = 4
3353  ns4 = ns4 + 1
3354  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
3355  ioa4(i,j,kwd) = 5
3356  ns5 = ns5 + 1
3357  ELSE IF(t .GT. 10000.) THEN
3358  ioa4(i,j,kwd) = 6
3359  ns6 = ns6 + 1
3360  ENDIF
3361  ENDDO
3362  ENDDO
3363  ENDDO
3364 C
3365  WRITE(6,*) "! MAKEOA2 EXIT"
3366 C
3367  RETURN
3368  END SUBROUTINE makeoa2
3369 
3378  function spherical_distance(theta1,phi1,theta2,phi2)
3379 
3380  real, intent(in) :: theta1, phi1, theta2, phi2
3381  real :: spherical_distance, dot
3382 
3383  if(theta1 == theta2 .and. phi1 == phi2) then
3384  spherical_distance = 0.0
3385  return
3386  endif
3387 
3388  dot = cos(phi1)*cos(phi2)*cos(theta1-theta2) + sin(phi1)*sin(phi2)
3389  if(dot > 1. ) dot = 1.
3390  if(dot < -1.) dot = -1.
3391  spherical_distance = acos(dot)
3392 
3393  return
3394 
3395  end function spherical_distance
3396 
3413  subroutine get_mismatch_index(im_in, jm_in, geolon_in,geolat_in,
3414  & bitmap_in,num_out, lon_out,lat_out, iindx, jindx )
3415  integer, intent(in) :: im_in, jm_in, num_out
3416  real, intent(in) :: geolon_in(im_in,jm_in)
3417  real, intent(in) :: geolat_in(im_in,jm_in)
3418  logical*1, intent(in) :: bitmap_in(im_in,jm_in)
3419  real, intent(in) :: lon_out(num_out), lat_out(num_out)
3420  integer, intent(out):: iindx(num_out), jindx(num_out)
3421  real, parameter :: MAX_DIST = 1.e+20
3422  integer, parameter :: NUMNBR = 20
3423  integer :: i_c,j_c,jstart,jend
3424  real :: lon,lat
3425 
3426  print*, "im_in,jm_in = ", im_in, jm_in
3427  print*, "num_out = ", num_out
3428  print*, "max and min of lon_in is", minval(geolon_in),
3429  & maxval(geolon_in)
3430  print*, "max and min of lat_in is", minval(geolat_in),
3431  & maxval(geolat_in)
3432  print*, "max and min of lon_out is", minval(lon_out),
3433  & maxval(lon_out)
3434  print*, "max and min of lat_out is", minval(lat_out),
3435  & maxval(lat_out)
3436  print*, "count(bitmap_in)= ", count(bitmap_in), max_dist
3437 
3438  do n = 1, num_out
3439  ! print*, "n = ", n
3440  lon = lon_out(n)
3441  lat = lat_out(n)
3442  !--- find the j-index for the nearest point
3443  i_c = 0; j_c = 0
3444  do j = 1, jm_in-1
3445  if(lat .LE. geolat_in(1,j) .and.
3446  & lat .GE. geolat_in(1,j+1)) then
3447  j_c = j
3448  endif
3449  enddo
3450  if(lat > geolat_in(1,1)) j_c = 1
3451  if(lat < geolat_in(1,jm_in)) j_c = jm_in
3452  ! print*, "lat =", lat, geolat_in(1,jm_in), geolat_in(1,jm_in-1)
3453  ! The input is Gaussian grid.
3454  jstart = max(j_c-numnbr,1)
3455  jend = min(j_c+numnbr,jm_in)
3456  dist = max_dist
3457  iindx(n) = 0
3458  jindx(n) = 0
3459  ! print*, "jstart, jend =", jstart, jend
3460  do j = jstart, jend; do i = 1,im_in
3461  if(bitmap_in(i,j) ) then
3462  ! print*, "bitmap_in is true"
3463  d = spherical_distance(lon_out(n),lat_out(n),
3464  & geolon_in(i,j), geolat_in(i,j))
3465  if( dist > d ) then
3466  iindx(n) = i; jindx(n) = j
3467  dist = d
3468  endif
3469  endif
3470  enddo; enddo
3471  if(iindx(n) ==0) then
3472  print*, "lon,lat=", lon,lat
3473  print*, "jstart, jend=", jstart, jend, dist
3474  print*, "ERROR in get mismatch_index: not find nearest points"
3475  call errexit(4)
3476  endif
3477  enddo
3478 
3479  end subroutine get_mismatch_index
3480 
3494  subroutine interpolate_mismatch(im_in, jm_in, data_in,
3495  & num_out, data_out, iindx, jindx)
3496  integer, intent(in) :: im_in, jm_in, num_out
3497  real, intent(in) :: data_in(im_in,jm_in)
3498  real, intent(out):: data_out(num_out)
3499  integer, intent(in) :: iindx(num_out), jindx(num_out)
3500 
3501  do n = 1, num_out
3502  data_out(n) = data_in(iindx(n),jindx(n))
3503  enddo
3504 
3505  end subroutine interpolate_mismatch
3506 
3551  SUBROUTINE makeoa3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
3552  1 ORO,SLM,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
3553  2 IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,
3554  3 is_south_pole,is_north_pole,IMI,JMI,OA_IN,OL_IN,
3555  4 slm_in,lon_in,lat_in)
3556  implicit none
3557  real, parameter :: MISSING_VALUE = -9999.
3558  real, parameter :: d2r = 3.14159265358979/180.
3559  real, PARAMETER :: R2D=180./3.14159265358979
3560  integer im,jm,imn,jmn,imi,jmi
3561  real GLAT(JMN)
3562  INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
3563  real SLM(IM,JM)
3564  real ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
3565  real OA4(IM,JM,4)
3566  integer IOA4(IM,JM,4)
3567  real OA_IN(IMI,JMI,4), OL_IN(IMI,JMI,4)
3568  real slm_in(IMI,JMI)
3569  real lon_in(IMI,JMI), lat_in(IMI,JMI)
3570  real lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
3571  real lon_t(IM,JM), lat_t(IM,JM)
3572  logical is_south_pole(IM,JM), is_north_pole(IM,JM)
3573  real XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
3574  real XNSUM3(IM,JM),XNSUM4(IM,JM)
3575  real VAR(IM,JM),OL(IM,JM,4)
3576  LOGICAL FLAG
3577  integer i,j,ilist(IMN),numx,i1,j1,ii1
3578  integer KWD,II,npts
3579  real LONO(4),LATO(4),LONI,LATI
3580  real DELXN,HC,HEIGHT,XNPU,XNPD,T
3581  integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
3582  logical inside_a_polygon
3583  real lon,lat,dlon,dlat,dlat_old
3584  real lon1,lat1,lon2,lat2
3585  real xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
3586  real HC_11, HC_12, HC_21, HC_22
3587  real xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
3588  real xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
3589  real get_lon_angle, get_lat_angle, get_xnsum
3590  integer ist, ien, jst, jen
3591  real xland,xwatr,xl1,xs1,oroavg
3592  integer int_opt, ipopt(20), kgds_input(200), kgds_output(200)
3593  integer count_land_output
3594  integer ij, ijmdl_output, iret, num_mismatch_land, num
3595  integer ibo(1)
3596  logical*1, allocatable :: bitmap_input(:,:)
3597  logical*1, allocatable :: bitmap_output(:)
3598  integer, allocatable :: ijsav_land_output(:)
3599  real, allocatable :: lats_land_output(:)
3600  real, allocatable :: lons_land_output(:)
3601  real, allocatable :: output_data_land(:)
3602  real, allocatable :: lons_mismatch_output(:)
3603  real, allocatable :: lats_mismatch_output(:)
3604  real, allocatable :: data_mismatch_output(:)
3605  integer, allocatable :: iindx(:), jindx(:)
3606 C
3607 C---- GLOBAL XLAT AND XLON ( DEGREE )
3608 C
3609  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
3610 C
3611  ijmdl_output = im*jm
3612 
3613  DO j=1,jmn
3614  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
3615  ENDDO
3616  print *,' IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
3617 C
3618 C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
3619 C
3620 C
3621  DO j=1,jm
3622  DO i=1,im
3623  xnsum(i,j) = 0.0
3624  elvmax(i,j) = oro(i,j)
3625  zmax(i,j) = 0.0
3626 C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
3627 C IN A GRID BOX
3628  xnsum1(i,j) = 0.0
3629  xnsum2(i,j) = 0.0
3630  xnsum3(i,j) = 0.0
3631  xnsum4(i,j) = 0.0
3632  oro1(i,j) = oro(i,j)
3633  elvmax(i,j) = zmax(i,j)
3634  ENDDO
3635  ENDDO
3636 
3637 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
3638 ! --- to JM or to JM1
3639  DO j=1,jm
3640 ! print*, "J=", J
3641  DO i=1,im
3642  hc = 1116.2 - 0.878 * var(i,j)
3643  lono(1) = lon_c(i,j)
3644  lono(2) = lon_c(i+1,j)
3645  lono(3) = lon_c(i+1,j+1)
3646  lono(4) = lon_c(i,j+1)
3647  lato(1) = lat_c(i,j)
3648  lato(2) = lat_c(i+1,j)
3649  lato(3) = lat_c(i+1,j+1)
3650  lato(4) = lat_c(i,j+1)
3651  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
3652  do j1 = jst, jen; do ii1 = 1, numx
3653  i1 = ilist(ii1)
3654  loni = i1*delxn
3655  lati = -90 + j1*delxn
3656  if(inside_a_polygon(loni*d2r,lati*d2r,4,
3657  & lono*d2r,lato*d2r))then
3658 
3659  height = float(zavg(i1,j1))
3660  IF(height.LT.-990.) height = 0.0
3661  IF ( height .gt. oro(i,j) ) then
3662  if ( height .gt. zmax(i,j) )zmax(i,j) = height
3663  ENDIF
3664  endif
3665  ENDDO ; ENDDO
3666  ENDDO
3667  ENDDO
3668 
3669 C
3670 ! --- this will make work1 array take on oro's values on return
3671 ! --- this will make work1 array take on oro's values on return
3672  DO j=1,jm
3673  DO i=1,im
3674 
3675  oro1(i,j) = oro(i,j)
3676  elvmax(i,j) = zmax(i,j)
3677  ENDDO
3678  ENDDO
3679 
3680  DO kwd = 1, 4
3681  DO j=1,jm
3682  DO i=1,im
3683  oa4(i,j,kwd) = 0.0
3684  ol(i,j,kwd) = 0.0
3685  ENDDO
3686  ENDDO
3687  ENDDO
3688 
3689  !--- use the nearest point to do remapping.
3690  int_opt = 2
3691  ipopt=0
3692  kgds_input = 0
3693  kgds_input(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
3694  kgds_input(2) = imi ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
3695  kgds_input(3) = jmi ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
3696  kgds_input(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
3697  kgds_input(5) = 0 ! OCT 14-16 - LON OF ORIGIN
3698  kgds_input(6) = 128 ! OCT 17 - RESOLUTION FLAG
3699  kgds_input(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
3700  kgds_input(8) = nint(-360000./imi) ! OCT 21-23 - LON OF EXTREME POINT
3701  kgds_input(9) = nint((360.0 / float(imi))*1000.0)
3702  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
3703  kgds_input(10) = jmi /2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
3704  kgds_input(12) = 255 ! OCT 29 - RESERVED
3705  kgds_input(20) = 255 ! OCT 5 - NOT USED, SET TO 255
3706 
3707 
3708  kgds_output = -1
3709 ! KGDS_OUTPUT(1) = IDRT ! OCT 6 - TYPE OF GRID (GAUSSIAN)
3710  kgds_output(2) = im ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
3711  kgds_output(3) = jm ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
3712  kgds_output(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
3713  kgds_output(5) = 0 ! OCT 14-16 - LON OF ORIGIN
3714  kgds_output(6) = 128 ! OCT 17 - RESOLUTION FLAG
3715  kgds_output(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
3716  kgds_output(8) = nint(-360000./im) ! OCT 21-23 - LON OF EXTREME POINT
3717  kgds_output(9) = nint((360.0 / float(im))*1000.0)
3718  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
3719  kgds_output(10) = jm /2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
3720  kgds_output(12) = 255 ! OCT 29 - RESERVED
3721  kgds_output(20) = 255 ! OCT 5 - NOT USED, SET TO 255
3722 
3723  count_land_output=0
3724  print*, "sum(slm) = ", sum(slm)
3725  do ij = 1, ijmdl_output
3726  j = (ij-1)/im + 1
3727  i = mod(ij-1,im) + 1
3728  if (slm(i,j) > 0.0) then
3729  count_land_output=count_land_output+1
3730  endif
3731  enddo
3732  allocate(bitmap_input(imi,jmi))
3733  bitmap_input=.false.
3734  print*, "number of land input=", sum(slm_in)
3735  where(slm_in > 0.0) bitmap_input=.true.
3736  print*, "count(bitmap_input)", count(bitmap_input)
3737 
3738  allocate(bitmap_output(count_land_output))
3739  allocate(output_data_land(count_land_output))
3740  allocate(ijsav_land_output(count_land_output))
3741  allocate(lats_land_output(count_land_output))
3742  allocate(lons_land_output(count_land_output))
3743 
3744 
3745 
3746  count_land_output=0
3747  do ij = 1, ijmdl_output
3748  j = (ij-1)/im + 1
3749  i = mod(ij-1,im) + 1
3750  if (slm(i,j) > 0.0) then
3751  count_land_output=count_land_output+1
3752  ijsav_land_output(count_land_output)=ij
3753  lats_land_output(count_land_output)=lat_t(i,j)
3754  lons_land_output(count_land_output)=lon_t(i,j)
3755  endif
3756  enddo
3757 
3758  oa4 = 0.0
3759  ol = 0.0
3760 
3761  do kwd=1,4
3762  bitmap_output = .false.
3763  output_data_land = 0.0
3764  call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3765  & (imi*jmi), count_land_output,
3766  & 1, 1, bitmap_input, oa_in(:,:,kwd),
3767  & count_land_output, lats_land_output,
3768  & lons_land_output, ibo,
3769  & bitmap_output, output_data_land, iret)
3770  if (iret /= 0) then
3771  print*,'- ERROR IN IPOLATES ',iret
3772  call errexit(4)
3773  endif
3774 
3775  num_mismatch_land = 0
3776  do ij = 1, count_land_output
3777  if (bitmap_output(ij)) then
3778  j = (ijsav_land_output(ij)-1)/im + 1
3779  i = mod(ijsav_land_output(ij)-1,im) + 1
3780  oa4(i,j,kwd)=output_data_land(ij)
3781  else ! default value
3782  num_mismatch_land = num_mismatch_land + 1
3783  endif
3784  enddo
3785  print*, "num_mismatch_land = ", num_mismatch_land
3786 
3787  if(.not. allocated(data_mismatch_output)) then
3788  allocate(lons_mismatch_output(num_mismatch_land))
3789  allocate(lats_mismatch_output(num_mismatch_land))
3790  allocate(data_mismatch_output(num_mismatch_land))
3791  allocate(iindx(num_mismatch_land))
3792  allocate(jindx(num_mismatch_land))
3793 
3794  num = 0
3795  do ij = 1, count_land_output
3796  if (.not. bitmap_output(ij)) then
3797  num = num+1
3798  lons_mismatch_output(num) = lons_land_output(ij)
3799  lats_mismatch_output(num) = lats_land_output(ij)
3800  endif
3801  enddo
3802 
3803  ! For those land points that with bitmap_output=.false. use
3804  ! the nearest land points to interpolate.
3805  print*,"before get_mismatch_index", count(bitmap_input)
3806  call get_mismatch_index(imi,jmi,lon_in*d2r,
3807  & lat_in*d2r,bitmap_input,num_mismatch_land,
3808  & lons_mismatch_output*d2r,lats_mismatch_output*d2r,
3809  & iindx, jindx )
3810  endif
3811 
3812  data_mismatch_output = 0
3813  call interpolate_mismatch(imi,jmi,oa_in(:,:,kwd),
3814  & num_mismatch_land,data_mismatch_output,iindx,jindx)
3815 
3816  num = 0
3817  do ij = 1, count_land_output
3818  if (.not. bitmap_output(ij)) then
3819  num = num+1
3820  j = (ijsav_land_output(ij)-1)/im + 1
3821  i = mod(ijsav_land_output(ij)-1,im) + 1
3822  oa4(i,j,kwd) = data_mismatch_output(num)
3823  if(i==372 .and. j== 611) then
3824  print*, "ij=",ij, num, oa4(i,j,kwd),iindx(num),jindx(num)
3825  endif
3826  endif
3827  enddo
3828 
3829 
3830  enddo
3831 
3832  !OL
3833  do kwd=1,4
3834  bitmap_output = .false.
3835  output_data_land = 0.0
3836  call ipolates(int_opt, ipopt, kgds_input, kgds_output,
3837  & (imi*jmi), count_land_output,
3838  & 1, 1, bitmap_input, ol_in(:,:,kwd),
3839  & count_land_output, lats_land_output,
3840  & lons_land_output, ibo,
3841  & bitmap_output, output_data_land, iret)
3842  if (iret /= 0) then
3843  print*,'- ERROR IN IPOLATES ',iret
3844  call errexit(4)
3845  endif
3846 
3847  num_mismatch_land = 0
3848  do ij = 1, count_land_output
3849  if (bitmap_output(ij)) then
3850  j = (ijsav_land_output(ij)-1)/im + 1
3851  i = mod(ijsav_land_output(ij)-1,im) + 1
3852  ol(i,j,kwd)=output_data_land(ij)
3853  else ! default value
3854  num_mismatch_land = num_mismatch_land + 1
3855  endif
3856  enddo
3857  print*, "num_mismatch_land = ", num_mismatch_land
3858 
3859  data_mismatch_output = 0
3860  call interpolate_mismatch(imi,jmi,ol_in(:,:,kwd),
3861  & num_mismatch_land,data_mismatch_output,iindx,jindx)
3862 
3863  num = 0
3864  do ij = 1, count_land_output
3865  if (.not. bitmap_output(ij)) then
3866  num = num+1
3867  j = (ijsav_land_output(ij)-1)/im + 1
3868  i = mod(ijsav_land_output(ij)-1,im) + 1
3869  ol(i,j,kwd) = data_mismatch_output(num)
3870  if(i==372 .and. j== 611) then
3871  print*, "ij=",ij, num, ol(i,j,kwd),iindx(num),jindx(num)
3872  endif
3873  endif
3874  enddo
3875 
3876 
3877  enddo
3878 
3879  deallocate(lons_mismatch_output,lats_mismatch_output)
3880  deallocate(data_mismatch_output,iindx,jindx)
3881  deallocate(bitmap_output,output_data_land)
3882  deallocate(ijsav_land_output,lats_land_output)
3883  deallocate(lons_land_output)
3884 
3885  DO kwd=1,4
3886  DO j=1,jm
3887  DO i=1,im
3888  t = oa4(i,j,kwd)
3889  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
3890  ENDDO
3891  ENDDO
3892  ENDDO
3893 C
3894  ns0 = 0
3895  ns1 = 0
3896  ns2 = 0
3897  ns3 = 0
3898  ns4 = 0
3899  ns5 = 0
3900  ns6 = 0
3901  DO kwd=1,4
3902  DO j=1,jm
3903  DO i=1,im
3904  t = abs( oa4(i,j,kwd) )
3905  IF(t .EQ. 0.) THEN
3906  ioa4(i,j,kwd) = 0
3907  ns0 = ns0 + 1
3908  ELSE IF(t .GT. 0. .AND. t .LE. 1.) THEN
3909  ioa4(i,j,kwd) = 1
3910  ns1 = ns1 + 1
3911  ELSE IF(t .GT. 1. .AND. t .LE. 10.) THEN
3912  ioa4(i,j,kwd) = 2
3913  ns2 = ns2 + 1
3914  ELSE IF(t .GT. 10. .AND. t .LE. 100.) THEN
3915  ioa4(i,j,kwd) = 3
3916  ns3 = ns3 + 1
3917  ELSE IF(t .GT. 100. .AND. t .LE. 1000.) THEN
3918  ioa4(i,j,kwd) = 4
3919  ns4 = ns4 + 1
3920  ELSE IF(t .GT. 1000. .AND. t .LE. 10000.) THEN
3921  ioa4(i,j,kwd) = 5
3922  ns5 = ns5 + 1
3923  ELSE IF(t .GT. 10000.) THEN
3924  ioa4(i,j,kwd) = 6
3925  ns6 = ns6 + 1
3926  ENDIF
3927  ENDDO
3928  ENDDO
3929  ENDDO
3930 C
3931  WRITE(6,*) "! MAKEOA3 EXIT"
3932 C
3933  RETURN
3934  END SUBROUTINE makeoa3
3935 
3946  SUBROUTINE revers(IM, JM, numi, F, WRK)
3947 !
3948  REAL F(IM,JM), WRK(IM,JM)
3949  integer numi(jm)
3950  imb2 = im / 2
3951  do i=1,im*jm
3952  wrk(i,1) = f(i,1)
3953  enddo
3954  do j=1,jm
3955  jr = jm - j + 1
3956  do i=1,im
3957  ir = i + imb2
3958  if (ir .gt. im) ir = ir - im
3959  f(ir,jr) = wrk(i,j)
3960  enddo
3961  enddo
3962 !
3963  tem = 0.0
3964  do i=1,im
3965  tem= tem + f(i,1)
3966  enddo
3967  tem = tem / im
3968  do i=1,im
3969  f(i,1) = tem
3970  enddo
3971 !
3972  RETURN
3973  END
3974 
3983  subroutine rg2gg(im,jm,numi,a)
3984  implicit none
3985  integer,intent(in):: im,jm,numi(jm)
3986  real,intent(inout):: a(im,jm)
3987  integer j,ir,ig
3988  real r,t(im)
3989  do j=1,jm
3990  r=real(numi(j))/real(im)
3991  do ig=1,im
3992  ir=mod(nint((ig-1)*r),numi(j))+1
3993  t(ig)=a(ir,j)
3994  enddo
3995  do ig=1,im
3996  a(ig,j)=t(ig)
3997  enddo
3998  enddo
3999  end subroutine
4000 
4009  subroutine gg2rg(im,jm,numi,a)
4010  implicit none
4011  integer,intent(in):: im,jm,numi(jm)
4012  real,intent(inout):: a(im,jm)
4013  integer j,ir,ig
4014  real r,t(im)
4015  do j=1,jm
4016  r=real(numi(j))/real(im)
4017  do ir=1,numi(j)
4018  ig=nint((ir-1)/r)+1
4019  t(ir)=a(ig,j)
4020  enddo
4021  do ir=1,numi(j)
4022  a(ir,j)=t(ir)
4023  enddo
4024  enddo
4025  end subroutine
4026 
4035  SUBROUTINE minmxj(IM,JM,A,title)
4036  implicit none
4037 
4038  real A(IM,JM),rmin,rmax
4039  integer i,j,IM,JM
4040  character*8 title
4041 
4042  rmin=1.e+10
4043  rmax=-rmin
4044 csela....................................................
4045 csela if(rmin.eq.1.e+10)return
4046 csela....................................................
4047  DO j=1,jm
4048  DO i=1,im
4049  if(a(i,j).ge.rmax)rmax=a(i,j)
4050  if(a(i,j).le.rmin)rmin=a(i,j)
4051  ENDDO
4052  ENDDO
4053  write(6,150)rmin,rmax,title
4054 150 format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
4055 C
4056  RETURN
4057  END
4058 
4070  SUBROUTINE mnmxja(IM,JM,A,imax,jmax,title)
4071  implicit none
4072 
4073  real A(IM,JM),rmin,rmax
4074  integer i,j,IM,JM,imax,jmax
4075  character*8 title
4076 
4077  rmin=1.e+10
4078  rmax=-rmin
4079 csela....................................................
4080 csela if(rmin.eq.1.e+10)return
4081 csela....................................................
4082  DO j=1,jm
4083  DO i=1,im
4084  if(a(i,j).ge.rmax)then
4085  rmax=a(i,j)
4086  imax=i
4087  jmax=j
4088  endif
4089  if(a(i,j).le.rmin)rmin=a(i,j)
4090  ENDDO
4091  ENDDO
4092  write(6,150)rmin,rmax,title
4093 150 format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
4094 C
4095  RETURN
4096  END
4097 
4132  SUBROUTINE spfft1(IMAX,INCW,INCG,KMAX,W,G,IDIR)
4133  IMPLICIT NONE
4134  INTEGER,INTENT(IN):: IMAX,INCW,INCG,KMAX,IDIR
4135  COMPLEX,INTENT(INOUT):: W(INCW,KMAX)
4136  REAL,INTENT(INOUT):: G(INCG,KMAX)
4137  REAL:: AUX1(25000+INT(0.82*IMAX))
4138  REAL:: AUX2(20000+INT(0.57*IMAX))
4139  INTEGER:: NAUX1,NAUX2
4140 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4141  NAUX1=25000+int(0.82*imax)
4142  naux2=20000+int(0.57*imax)
4143 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4144 C FOURIER TO PHYSICAL TRANSFORM.
4145  SELECT CASE(idir)
4146  CASE(1:)
4147  SELECT CASE(digits(1.))
4148  CASE(digits(1._4))
4149  CALL scrft(1,w,incw,g,incg,imax,kmax,-1,1.,
4150  & aux1,naux1,aux2,naux2,0.,0)
4151  CALL scrft(0,w,incw,g,incg,imax,kmax,-1,1.,
4152  & aux1,naux1,aux2,naux2,0.,0)
4153  CASE(digits(1._8))
4154  CALL dcrft(1,w,incw,g,incg,imax,kmax,-1,1.,
4155  & aux1,naux1,aux2,naux2,0.,0)
4156  CALL dcrft(0,w,incw,g,incg,imax,kmax,-1,1.,
4157  & aux1,naux1,aux2,naux2,0.,0)
4158  END SELECT
4159 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4160 C PHYSICAL TO FOURIER TRANSFORM.
4161  CASE(:-1)
4162  SELECT CASE(digits(1.))
4163  CASE(digits(1._4))
4164  CALL srcft(1,g,incg,w,incw,imax,kmax,+1,1./imax,
4165  & aux1,naux1,aux2,naux2,0.,0)
4166  CALL srcft(0,g,incg,w,incw,imax,kmax,+1,1./imax,
4167  & aux1,naux1,aux2,naux2,0.,0)
4168  CASE(digits(1._8))
4169  CALL drcft(1,g,incg,w,incw,imax,kmax,+1,1./imax,
4170  & aux1,naux1,aux2,naux2,0.,0)
4171  CALL drcft(0,g,incg,w,incw,imax,kmax,+1,1./imax,
4172  & aux1,naux1,aux2,naux2,0.,0)
4173  END SELECT
4174  END SELECT
4175  END SUBROUTINE
4176 
4182  subroutine read_g(glob,ITOPO)
4183  implicit none
4184 cc
4185  integer*2 glob(360*120,180*120)
4186 cc
4187  integer ix,jx
4188  integer ia,ja
4189 cc
4190  parameter(ix=40*120,jx=50*120)
4191  parameter(ia=60*120,ja=30*120)
4192 cc
4193  integer*2 idat(ix,jx)
4194  integer itopo
4195 cc
4196  integer i,j,inttyp
4197 cc
4198  real(kind=8) dloin,dlain,rlon,rlat
4199 cc
4200  open(235, file="./fort.235", access='direct', recl=43200*21600*2)
4201  read(235,rec=1)glob
4202  close(235)
4203 cc
4204  print*,' '
4205  call maxmin (glob,360*120*180*120,'global0')
4206 cc
4207 cc
4208  dloin=1.d0/120.d0
4209  dlain=1.d0/120.d0
4210 cc
4211  rlon= -179.995833333333333333333333d0
4212  rlat= 89.995833333333333333333333d0
4213 cc
4214  inttyp=-1 ! average rectangular subset
4215 ccmr inttyp= 1 ! take closest grid point value
4216 ccmr inttyp= 0 ! interpolate from four closest grid point values
4217 cc
4218 ! call la2ga_gtopo30(glob,360*120,180*120,
4219 ! & dloin,dlain,rlon,rlat,inttyp,
4220 ! & .true.,glob,
4221 ! & 0,lonf,latg)
4222 cc
4223  return
4224  end
4225 
4233  subroutine maxmin(ia,len,tile)
4234 ccmr
4235  implicit none
4236 ccmr
4237  integer*2 ia(len)
4238  character*7 tile
4239  integer iaamax, iaamin, len, j, m, ja, kount
4240  integer(8) sum2,std,mean,isum
4241  integer i_count_notset,kount_9
4242 ! --- missing is -9999
4243 c
4244  isum = 0
4245  sum2 = 0
4246  kount = 0
4247  kount_9 = 0
4248  iaamax = -9999999
4249 ccmr iaamin = 1
4250  iaamin = 9999999
4251  i_count_notset=0
4252  do 10 m=1,len
4253  ja=ia(m)
4254 ccmr if ( ja .lt. 0 ) print *,' ja < 0:',ja
4255 ccmr if ( ja .eq. -9999 ) goto 10
4256  if ( ja .eq. -9999 ) then
4257  kount_9=kount_9+1
4258  goto 10
4259  endif
4260  if ( ja .eq. -12345 ) i_count_notset=i_count_notset+1
4261 ccmr if ( ja .eq. 0 ) goto 11
4262  iaamax = max0( iaamax, ja )
4263  iaamin = min0( iaamin, ja )
4264 ! iaamax = max0( iaamax, ia(m,j) )
4265 ! iaamin = min0( iaamin, ia(m,j) )
4266  11 continue
4267  kount = kount + 1
4268  isum = isum + ja
4269 ccmr sum2 = sum2 + ifix( float(ja) * float(ja) )
4270  sum2 = sum2 + ja*ja
4271  10 continue
4272 !
4273  mean = isum/kount
4274  std = ifix(sqrt(float((sum2/(kount))-mean**2)))
4275  print*,tile,' max=',iaamax,' min=',iaamin,' sum=',isum,
4276  & ' i_count_notset=',i_count_notset
4277  print*,tile,' mean=',mean,' std.dev=',std,
4278  & ' ko9s=',kount,kount_9,kount+kount_9
4279  return
4280  end
4281 
4291  SUBROUTINE minmaxj(IM,JM,A,title)
4292  implicit none
4293 
4294  real(kind=4) A(IM,JM),rmin,rmax,undef
4295  integer i,j,IM,JM,imax,jmax,imin,jmin,iundef
4296  character*8 title,chara
4297  data chara/' '/
4298  chara=title
4299  rmin=1.e+10
4300  rmax=-rmin
4301  imax=0
4302  imin=0
4303  jmax=0
4304  jmin=0
4305  iundef=0
4306  undef=-9999.
4307 csela....................................................
4308 csela if(rmin.eq.1.e+10)return
4309 csela....................................................
4310  DO j=1,jm
4311  DO i=1,im
4312  if(a(i,j).ge.rmax)then
4313  rmax=a(i,j)
4314  imax=i
4315  jmax=j
4316  endif
4317  if(a(i,j).le.rmin)then
4318  if ( a(i,j) .eq. undef ) then
4319  iundef = iundef + 1
4320  else
4321  rmin=a(i,j)
4322  imin=i
4323  jmin=j
4324  endif
4325  endif
4326  ENDDO
4327  ENDDO
4328  write(6,150)chara,rmin,imin,jmin,rmax,imax,jmax,iundef
4329 150 format(1x,a8,2x,'rmin=',e13.4,2i6,2x,'rmax=',e13.4,3i6)
4330 C
4331  RETURN
4332  END
4333 
4343  subroutine latlon2xyz(siz,lon, lat, x, y, z)
4344  implicit none
4345  integer, intent(in) :: siz
4346  real, intent(in) :: lon(siz), lat(siz)
4347  real, intent(out) :: x(siz), y(siz), z(siz)
4348 
4349  integer n
4350 
4351  do n = 1, siz
4352  x(n) = cos(lat(n))*cos(lon(n))
4353  y(n) = cos(lat(n))*sin(lon(n))
4354  z(n) = sin(lat(n))
4355  enddo
4356  end
4357 
4365  FUNCTION spherical_angle(v1, v2, v3)
4366  implicit none
4367  real, parameter :: epsln30 = 1.e-30
4368  real, parameter :: pi=3.1415926535897931
4369  real v1(3), v2(3), v3(3)
4370  real spherical_angle
4371 
4372  real px, py, pz, qx, qy, qz, ddd;
4373 
4374  ! vector product between v1 and v2
4375  px = v1(2)*v2(3) - v1(3)*v2(2)
4376  py = v1(3)*v2(1) - v1(1)*v2(3)
4377  pz = v1(1)*v2(2) - v1(2)*v2(1)
4378  ! vector product between v1 and v3
4379  qx = v1(2)*v3(3) - v1(3)*v3(2);
4380  qy = v1(3)*v3(1) - v1(1)*v3(3);
4381  qz = v1(1)*v3(2) - v1(2)*v3(1);
4382 
4383  ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
4384  if ( ddd <= 0.0 ) then
4385  spherical_angle = 0.
4386  else
4387  ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
4388  if( abs(ddd-1) < epsln30 ) ddd = 1;
4389  if( abs(ddd+1) < epsln30 ) ddd = -1;
4390  if ( ddd>1. .or. ddd<-1. ) then
4391  !FIX to correctly handle co-linear points (angle near pi or 0) */
4392  if (ddd < 0.) then
4393  spherical_angle = pi
4394  else
4395  spherical_angle = 0.
4396  endif
4397  else
4398  spherical_angle = acos( ddd )
4399  endif
4400  endif
4401 
4402  return
4403  END
4404 
4415  FUNCTION inside_a_polygon(lon1, lat1, npts, lon2, lat2)
4416  implicit none
4417  real, parameter :: epsln10 = 1.e-10
4418  real, parameter :: epsln8 = 1.e-8
4419  real, parameter :: pi=3.1415926535897931
4420  real, parameter :: range_check_criteria=0.05
4421  real :: anglesum, angle, spherical_angle
4422  integer i, ip1
4423  real lon1, lat1
4424  integer npts
4425  real lon2(npts), lat2(npts)
4426  real x2(npts), y2(npts), z2(npts)
4427  real lon1_1d(1), lat1_1d(1)
4428  real x1(1), y1(1), z1(1)
4429  real pnt0(3),pnt1(3),pnt2(3)
4430  logical inside_a_polygon
4431  real max_x2,min_x2,max_y2,min_y2,max_z2,min_z2
4432  !first convert to cartesian grid */
4433  call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
4434  lon1_1d(1) = lon1
4435  lat1_1d(1) = lat1
4436  call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
4437  inside_a_polygon = .false.
4438  max_x2 = maxval(x2)
4439  if( x1(1) > max_x2+range_check_criteria ) return
4440  min_x2 = minval(x2)
4441  if( x1(1)+range_check_criteria < min_x2 ) return
4442  max_y2 = maxval(y2)
4443  if( y1(1) > max_y2+range_check_criteria ) return
4444  min_y2 = minval(y2)
4445  if( y1(1)+range_check_criteria < min_y2 ) return
4446  max_z2 = maxval(z2)
4447  if( z1(1) > max_z2+range_check_criteria ) return
4448  min_z2 = minval(z2)
4449  if( z1(1)+range_check_criteria < min_z2 ) return
4450 
4451  pnt0(1) = x1(1)
4452  pnt0(2) = y1(1)
4453  pnt0(3) = z1(1)
4454 
4455  anglesum = 0;
4456  do i = 1, npts
4457  if(abs(x1(1)-x2(i)) < epsln10 .and.
4458  & abs(y1(1)-y2(i)) < epsln10 .and.
4459  & abs(z1(1)-z2(i)) < epsln10 ) then ! same as the corner point
4460  inside_a_polygon = .true.
4461  return
4462  endif
4463  ip1 = i+1
4464  if(ip1>npts) ip1 = 1
4465  pnt1(1) = x2(i)
4466  pnt1(2) = y2(i)
4467  pnt1(3) = z2(i)
4468  pnt2(1) = x2(ip1)
4469  pnt2(2) = y2(ip1)
4470  pnt2(3) = z2(ip1)
4471 
4472  angle = spherical_angle(pnt0, pnt2, pnt1);
4473 ! anglesum = anglesum + spherical_angle(pnt0, pnt2, pnt1);
4474  anglesum = anglesum + angle
4475  enddo
4476 
4477  if(abs(anglesum-2*pi) < epsln8) then
4478  inside_a_polygon = .true.
4479  else
4480  inside_a_polygon = .false.
4481  endif
4482 
4483  return
4484 
4485  end
4486 
4510  function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
4511  & glat,zavg,zslm,delxn)
4512  implicit none
4513 
4514  real get_xnsum
4515  logical verbose
4516  real lon1,lat1,lon2,lat2,oro,delxn
4517  integer imn,jmn
4518  real glat(jmn)
4519  integer zavg(imn,jmn),zslm(imn,jmn)
4520  integer i, j, ist, ien, jst, jen, i1
4521  real height
4522  real xland,xwatr,xl1,xs1,slm,xnsum
4523  !---figure out ist,ien,jst,jen
4524  do j = 1, jmn
4525  if( glat(j) .GT. lat1 ) then
4526  jst = j
4527  exit
4528  endif
4529  enddo
4530  do j = 1, jmn
4531  if( glat(j) .GT. lat2 ) then
4532  jen = j
4533  exit
4534  endif
4535  enddo
4536 
4537 
4538  ist = lon1/delxn + 1
4539  ien = lon2/delxn
4540  if(ist .le.0) ist = ist + imn
4541  if(ien < ist) ien = ien + imn
4542 ! if(verbose) print*, "ist,ien=",ist,ien,jst,jen
4543 
4544  !--- compute average oro
4545  oro = 0.0
4546  xnsum = 0
4547  xland = 0
4548  xwatr = 0
4549  xl1 = 0
4550  xs1 = 0
4551  do j = jst,jen
4552  do i1 = 1, ien - ist + 1
4553  i = ist + i1 -1
4554  if( i .LE. 0) i = i + imn
4555  if( i .GT. imn) i = i - imn
4556  xland = xland + float(zslm(i,j))
4557  xwatr = xwatr + float(1-zslm(i,j))
4558  xnsum = xnsum + 1.
4559  height = float(zavg(i,j))
4560  IF(height.LT.-990.) height = 0.0
4561  xl1 = xl1 + height * float(zslm(i,j))
4562  xs1 = xs1 + height * float(1-zslm(i,j))
4563  enddo
4564  enddo
4565  if( xnsum > 1.) THEN
4566  slm = float(nint(xland/xnsum))
4567  IF(slm.NE.0.) THEN
4568  oro= xl1 / xland
4569  ELSE
4570  oro = xs1 / xwatr
4571  ENDIF
4572  ENDIF
4573 
4574  get_xnsum = 0
4575  do j = jst, jen
4576  do i1= 1, ien-ist+1
4577  i = ist + i1 -1
4578  if( i .LE. 0) i = i + imn
4579  if( i .GT. imn) i = i - imn
4580  height = float(zavg(i,j))
4581  IF(height.LT.-990.) height = 0.0
4582  IF ( height .gt. oro ) get_xnsum = get_xnsum + 1
4583  enddo
4584  enddo
4585 ! if(verbose) print*, "get_xnsum=", get_xnsum, oro
4586 
4587  end function get_xnsum
4588 
4617  subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
4618  & glat,zavg,zslm,delxn,xnsum1,xnsum2,HC)
4619  implicit none
4620 
4621  real, intent(out) :: xnsum1,xnsum2,HC
4622  logical verbose
4623  real lon1,lat1,lon2,lat2,oro,delxn
4624  integer IMN,JMN
4625  real glat(JMN)
4626  integer zavg(IMN,JMN),zslm(IMN,JMN)
4627  integer i, j, ist, ien, jst, jen, i1
4628  real HEIGHT, var
4629  real XW1,XW2,slm,xnsum
4630  !---figure out ist,ien,jst,jen
4631  do j = 1, jmn
4632  if( glat(j) .GT. lat1 ) then
4633  jst = j
4634  exit
4635  endif
4636  enddo
4637  do j = 1, jmn
4638  if( glat(j) .GT. lat2 ) then
4639  jen = j
4640  exit
4641  endif
4642  enddo
4643 
4644 
4645  ist = lon1/delxn + 1
4646  ien = lon2/delxn
4647  if(ist .le.0) ist = ist + imn
4648  if(ien < ist) ien = ien + imn
4649 ! if(verbose) print*, "ist,ien=",ist,ien,jst,jen
4650 
4651  !--- compute average oro
4652  xnsum = 0
4653  xw1 = 0
4654  xw2 = 0
4655  do j = jst,jen
4656  do i1 = 1, ien - ist + 1
4657  i = ist + i1 -1
4658  if( i .LE. 0) i = i + imn
4659  if( i .GT. imn) i = i - imn
4660  xnsum = xnsum + 1.
4661  height = float(zavg(i,j))
4662  IF(height.LT.-990.) height = 0.0
4663  xw1 = xw1 + height
4664  xw2 = xw2 + height ** 2
4665  enddo
4666  enddo
4667  var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
4668  hc = 1116.2 - 0.878 * var
4669  xnsum1 = 0
4670  xnsum2 = 0
4671  do j = jst, jen
4672  do i1= 1, ien-ist+1
4673  i = ist + i1 -1
4674  if( i .LE. 0) i = i + imn
4675  if( i .GT. imn) i = i - imn
4676  height = float(zavg(i,j))
4677  IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4678  xnsum2 = xnsum2 + 1
4679  enddo
4680  enddo
4681 
4682  end subroutine get_xnsum2
4683 
4713  subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
4714  & glat,zavg,zslm,delxn,xnsum1,xnsum2,HC)
4715  implicit none
4716 
4717  real, intent(out) :: xnsum1,xnsum2
4718  real lon1,lat1,lon2,lat2,oro,delxn
4719  integer IMN,JMN
4720  real glat(JMN)
4721  integer zavg(IMN,JMN),zslm(IMN,JMN)
4722  integer i, j, ist, ien, jst, jen, i1
4723  real HEIGHT, HC
4724  real XW1,XW2,slm,xnsum
4725  !---figure out ist,ien,jst,jen
4726  ! if lat1 or lat 2 is 90 degree. set jst = JMN
4727  jst = jmn
4728  jen = jmn
4729  do j = 1, jmn
4730  if( glat(j) .GT. lat1 ) then
4731  jst = j
4732  exit
4733  endif
4734  enddo
4735  do j = 1, jmn
4736  if( glat(j) .GT. lat2 ) then
4737  jen = j
4738  exit
4739  endif
4740  enddo
4741 
4742 
4743  ist = lon1/delxn + 1
4744  ien = lon2/delxn
4745  if(ist .le.0) ist = ist + imn
4746  if(ien < ist) ien = ien + imn
4747 ! if(verbose) print*, "ist,ien=",ist,ien,jst,jen
4748 
4749  xnsum1 = 0
4750  xnsum2 = 0
4751  do j = jst, jen
4752  do i1= 1, ien-ist+1
4753  i = ist + i1 -1
4754  if( i .LE. 0) i = i + imn
4755  if( i .GT. imn) i = i - imn
4756  height = float(zavg(i,j))
4757  IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
4758  xnsum2 = xnsum2 + 1
4759  enddo
4760  enddo
4761 
4762  end subroutine get_xnsum3
4763 
4777  subroutine nanc(a,l,c)
4778  integer inan1,inan2,inan3,inan4,inaq1,inaq2,inaq3,inaq4
4779  real word
4780  integer itest
4781  equivalence (itest,word)
4782 c
4783 c signaling NaN
4784  data inan1/x'7F800001'/
4785  data inan2/x'7FBFFFFF'/
4786  data inan3/x'FF800001'/
4787  data inan4/x'FFBFFFFF'/
4788 c
4789 c quiet NaN
4790 c
4791  data inaq1/x'7FC00000'/
4792  data inaq2/x'7FFFFFFF'/
4793  data inaq3/x'FFC00000'/
4794  data inaq4/x'FFFFFFFF'/
4795 c
4796  real(kind=8)a(l),rtc,t1,t2
4797  character*24 cn
4798  character*(*) c
4799 c t1=rtc()
4800 cgwv print *, ' nanc call ',c
4801  do k=1,l
4802  word=a(k)
4803  if( (itest .GE. inan1 .AND. itest .LE. inan2) .OR.
4804  * (itest .GE. inan3 .AND. itest .LE. inan4) ) then
4805  print *,' NaNs detected at word',k,' ',c
4806  return
4807  endif
4808  if( (itest .GE. inaq1 .AND. itest .LE. inaq2) .OR.
4809  * (itest .GE. inaq3 .AND. itest .LE. inaq4) ) then
4810  print *,' NaNq detected at word',k,' ',c
4811  return
4812  endif
4813 
4814  101 format(e20.10)
4815  end do
4816 c t2=rtc()
4817 cgwv print 102,l,t2-t1,c
4818  102 format(' time to check ',i9,' words is ',f10.4,' ',a24)
4819  return
4820  end
4821 
4826  real function timef()
4827  character(8) :: date
4828  character(10) :: time
4829  character(5) :: zone
4830  integer,dimension(8) :: values
4831  integer :: total
4832  real :: elapsed
4833  call date_and_time(date,time,zone,values)
4834  total=(3600*values(5))+(60*values(6))
4835  * +values(7)
4836  elapsed=float(total) + (1.0e-3*float(values(8)))
4837  timef=elapsed
4838  return
4839  end
get_xnsum3
subroutine get_xnsum3(lon1, lat1, lon2, lat2, IMN, JMN, glat, zavg, zslm, delxn, xnsum1, xnsum2, HC)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
Definition: mtnlm7_oclsm.f:4715
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:4416
spherical_distance
real function spherical_distance(theta1, phi1, theta2, phi2)
Compute a great circle distance between two points.
Definition: mtnlm7_oclsm.f:3379
timef
real function timef()
Get the date/time for the system clock.
Definition: mtnlm7_oclsm.f:4827
netcdf_err
subroutine netcdf_err(err, string)
Check NetCDF error code and output the error message.
Definition: netcdf_io.F90:219
latlon2xyz
subroutine latlon2xyz(siz, lon, lat, x, y, z)
Convert from latitude and longitude to x,y,z coordinates.
Definition: mtnlm7_oclsm.f:4344
gg2rg
subroutine gg2rg(im, jm, numi, a)
Convert from a full grid to a reduced grid.
Definition: mtnlm7_oclsm.f:4010
rg2gg
subroutine rg2gg(im, jm, numi, a)
Convert from a reduced grid to a full grid.
Definition: mtnlm7_oclsm.f:3984
spherical_angle
real function spherical_angle(v1, v2, v3)
Compute spherical angle.
Definition: mtnlm7_oclsm.f:4366
nanc
subroutine nanc(a, l, c)
Report NaNS and NaNQ within an address range for 8-byte real words.
Definition: mtnlm7_oclsm.f:4778
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:4292
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
minmxj
subroutine minmxj(IM, JM, A, title)
Print out the maximum and minimum values of an array.
Definition: mtnlm7_oclsm.f:4036
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:2938
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
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:4619
mnmxja
subroutine mnmxja(IM, JM, A, imax, jmax, title)
Print out the maximum and minimum values of an array.
Definition: mtnlm7_oclsm.f:4071
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:4512
maxmin
subroutine maxmin(ia, len, tile)
Print the maximum, mininum, mean and standard deviation of an array.
Definition: mtnlm7_oclsm.f:4234
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:3496
read_g
subroutine read_g(glob, ITOPO)
Read input global 30-arc second orography data.
Definition: mtnlm7_oclsm.f:4183
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:3415
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:2919