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