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