orog_mask_tools  1.13.0
•All Data Structures Namespaces Files Functions Variables Pages
mtnlm7_oclsm.F90
Go to the documentation of this file.
1 
4 
71 
72  use io_utils, only : read_mdl_dims
73  implicit none
74 
75  character(len=256) :: mdl_grid_file = "none"
76  character(len=256) :: external_mask_file = "none"
77  integer :: im, jm, efac
78  logical :: mask_only = .false.
79 
80  print*,"- BEGIN OROGRAPHY PROGRAM."
81 
82  read(5,*) mdl_grid_file
83  read(5,*) mask_only
84  read(5,*) external_mask_file
85 
86  efac = 0
87 
88  if (mask_only) then
89  print*,"- WILL COMPUTE LANDMASK ONLY."
90  endif
91 
92  if (trim(external_mask_file) /= "none") then
93  print*,"- WILL USE EXTERNAL LANDMASK FROM FILE: ", trim(external_mask_file)
94  endif
95 
96  call read_mdl_dims(mdl_grid_file, im, jm)
97 
98  call tersub(im,jm,efac,mdl_grid_file,mask_only,external_mask_file)
99 
100  print*,"- NORMAL TERMINATION."
101 
102  stop
103  end
104 
118  SUBROUTINE tersub(IM,JM,EFAC,OUTGRID,MASK_ONLY,EXTERNAL_MASK_FILE)
119 
120  use io_utils, only : qc_orog_by_ramp, write_mask_netcdf, &
121  read_global_mask, read_global_orog, &
122  read_mask, write_netcdf, &
123  read_mdl_grid_file
124  use orog_utils, only : minmax, timef, remove_isolated_pts
125 
126  implicit none
127 
128  integer, parameter :: imn = 360*120
129  integer, parameter :: jmn = 180*120
130 
131  integer, intent(in) :: IM,JM,efac
132  character(len=*), intent(in) :: OUTGRID
133  character(len=*), intent(in) :: EXTERNAL_MASK_FILE
134 
135  logical, intent(in) :: mask_only
136 
137  integer :: i,j
138  integer :: itest,jtest
139 
140  integer, allocatable :: ZAVG(:,:),ZSLM(:,:)
141  integer(1), allocatable :: UMD(:,:)
142  integer(2), allocatable :: glob(:,:)
143 
144  real :: tbeg,tend,tbeg1
145 
146  real, allocatable :: XLAT(:),XLON(:)
147  real, allocatable :: GEOLON(:,:),GEOLON_C(:,:),DX(:,:)
148  real, allocatable :: GEOLAT(:,:),GEOLAT_C(:,:),DY(:,:)
149  real, allocatable :: SLM(:,:),ORO(:,:),VAR(:,:)
150  real, allocatable :: land_frac(:,:),lake_frac(:,:)
151  real, allocatable :: THETA(:,:),GAMMA(:,:),SIGMA(:,:),ELVMAX(:,:)
152  real, allocatable :: VAR4(:,:)
153  real, allocatable :: OA(:,:,:),OL(:,:,:),HPRIME(:,:,:)
154 
155  logical :: is_south_pole(IM,JM), is_north_pole(IM,JM)
156 
157  tbeg1=timef()
158  tbeg=timef()
159 
160  allocate (glob(imn,jmn))
161  allocate (zavg(imn,jmn))
162  allocate (zslm(imn,jmn))
163  allocate (umd(imn,jmn))
164 
165 ! Read global mask data.
166 
167  call read_global_mask(imn,jmn,umd)
168 
169 ! Read global orography data.
170 
171  call read_global_orog(imn,jmn,glob)
172 
173 ! ZSLM initialize with all land (1). Ocean is '0'.
174 
175  zslm=1
176 
177 ! ZAVG initialize from glob
178 
179  zavg=glob
180 
181  do j=1,jmn
182  do i=1,imn
183  if ( umd(i,j) .eq. 0 ) zslm(i,j) = 0
184  enddo
185  enddo
186 
187  deallocate (umd,glob)
188 
189 ! Fixing an error in the topo 30" data set at pole (-9999).
190 
191  do i=1,imn
192  zslm(i,1)=0
193  zslm(i,jmn)=1
194  enddo
195 
196 ! Quality control the global topography data over Antarctica
197 ! using RAMP data.
198 
199  call qc_orog_by_ramp(imn, jmn, zavg, zslm)
200 
201  allocate (geolon(im,jm),geolon_c(im+1,jm+1),dx(im,jm))
202  allocate (geolat(im,jm),geolat_c(im+1,jm+1),dy(im,jm))
203  allocate (slm(im,jm))
204  allocate (land_frac(im,jm),lake_frac(im,jm))
205 
206 ! Reading grid file.
207 
208  call read_mdl_grid_file(outgrid,im,jm,geolon,geolon_c, &
209  geolat,geolat_c,dx,dy,is_north_pole,is_south_pole)
210 
211  tend=timef()
212  print*,"- TIMING: READING INPUT DATA ",tend-tbeg
213  !
214  tbeg=timef()
215 
216  IF (external_mask_file == 'none') then
217  CALL make_mask(zslm,slm,land_frac, &
218  im,jm,imn,jmn,geolon_c,geolat_c)
219  lake_frac=9999.9
220  ELSE
221  CALL read_mask(external_mask_file,slm,land_frac, &
222  lake_frac,im,jm)
223  ENDIF
224 
225  IF (mask_only) THEN
226  print*,'- WILL COMPUTE LANDMASK ONLY.'
227  CALL write_mask_netcdf(im,jm,slm,land_frac, &
228  1,1,geolon,geolat)
229 
230  DEALLOCATE(zavg, zslm, slm, land_frac, lake_frac)
231  DEALLOCATE(geolon, geolon_c, geolat, geolat_c)
232  print*,'- NORMAL TERMINATION.'
233  stop
234  END IF
235 
236  allocate (var(im,jm),var4(im,jm),oro(im,jm))
237 
238  CALL makemt2(zavg,zslm,oro,slm,var,var4, &
239  im,jm,imn,jmn,geolon_c,geolat_c,lake_frac,land_frac)
240 
241  tend=timef()
242  print*,"- TIMING: MASK AND OROG CREATION ", tend-tbeg
243 
244  call minmax(im,jm,oro,'ORO ')
245  call minmax(im,jm,slm,'SLM ')
246  call minmax(im,jm,var,'VAR ')
247  call minmax(im,jm,var4,'VAR4 ')
248 
249 ! Compute mtn principal coord HTENSR: THETA,GAMMA,SIGMA
250 
251  allocate (theta(im,jm),gamma(im,jm),sigma(im,jm),elvmax(im,jm))
252 
253  tbeg=timef()
254  CALL makepc2(zavg,zslm,theta,gamma,sigma, &
255  im,jm,imn,jmn,geolon_c,geolat_c,slm)
256  tend=timef()
257 
258  print*,"- TIMING: CREATE PRINCIPLE COORDINATE ",tend-tbeg
259 
260  call minmax(im,jm,theta,'THETA ')
261  call minmax(im,jm,gamma,'GAMMA ')
262  call minmax(im,jm,sigma,'SIGMA ')
263 
264 ! COMPUTE MOUNTAIN DATA : OA OL
265 
266  allocate (oa(im,jm,4),ol(im,jm,4))
267 
268  tbeg=timef()
269  CALL makeoa2(zavg,zslm,var,oa,ol,elvmax,oro, &
270  im,jm,imn,jmn,geolon_c,geolat_c, &
271  geolon,geolat,dx,dy,is_south_pole,is_north_pole)
272 
273  tend=timef()
274 
275  print*,"- TIMING: CREATE ASYMETRY AND LENGTH SCALE ",tend-tbeg
276 
277  deallocate (zslm,zavg)
278  deallocate (dx,dy)
279 
280  tbeg=timef()
281  call minmax(im,jm,oa,'OA ')
282  call minmax(im,jm,ol,'OL ')
283  call minmax(im,jm,elvmax,'ELVMAX ')
284  call minmax(im,jm,oro,'ORO ')
285 
286 ! Replace maximum elevation with max elevation minus orography.
287 ! If maximum elevation is less than the orography, replace with
288 ! a proxy.
289 
290  print*,"- QC MAXIMUM ELEVATION."
291  DO j = 1,jm
292  DO i = 1,im
293  if (elvmax(i,j) .lt. oro(i,j) ) then
294  elvmax(i,j) = max( 3. * var(i,j),0.)
295  else
296  elvmax(i,j) = max( elvmax(i,j) - oro(i,j),0.)
297  endif
298  ENDDO
299  ENDDO
300 
301  call minmax(im,jm,elvmax,'ELVMAX ',itest,jtest)
302 
303  print*,"- ZERO FIELDS OVER OCEAN."
304 
305  DO j = 1,jm
306  DO i = 1,im
307  IF(slm(i,j).EQ.0.) THEN
308 ! VAR(I,J) = 0.
309  var4(i,j) = 0.
310  oa(i,j,1) = 0.
311  oa(i,j,2) = 0.
312  oa(i,j,3) = 0.
313  oa(i,j,4) = 0.
314  ol(i,j,1) = 0.
315  ol(i,j,2) = 0.
316  ol(i,j,3) = 0.
317  ol(i,j,4) = 0.
318 ! THETA(I,J) =0.
319 ! GAMMA(I,J) =0.
320 ! SIGMA(I,J) =0.
321 ! ELVMAX(I,J)=0.
322 ! --- the sub-grid scale parameters for mtn blocking and gwd retain
323 ! --- properties even if over ocean but there is elevation within the
324 ! --- gaussian grid box.
325  ENDIF
326  ENDDO
327  ENDDO
328 
329  IF (external_mask_file == 'none') then
330 
331  call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol)
332 
333  endif
334 
335  allocate(hprime(im,jm,14))
336 
337  DO j=1,jm
338  DO i=1,im
339  oro(i,j) = oro(i,j) + efac*var(i,j)
340  hprime(i,j,1) = var(i,j)
341  hprime(i,j,2) = var4(i,j)
342  hprime(i,j,3) = oa(i,j,1)
343  hprime(i,j,4) = oa(i,j,2)
344  hprime(i,j,5) = oa(i,j,3)
345  hprime(i,j,6) = oa(i,j,4)
346  hprime(i,j,7) = ol(i,j,1)
347  hprime(i,j,8) = ol(i,j,2)
348  hprime(i,j,9) = ol(i,j,3)
349  hprime(i,j,10)= ol(i,j,4)
350  hprime(i,j,11)= theta(i,j)
351  hprime(i,j,12)= gamma(i,j)
352  hprime(i,j,13)= sigma(i,j)
353  hprime(i,j,14)= elvmax(i,j)
354  ENDDO
355  ENDDO
356 
357  deallocate(var4)
358 
359  call minmax(im,jm,elvmax,'ELVMAX ',itest,jtest)
360  call minmax(im,jm,oro,'ORO ')
361 
362  print *,'- ORO(itest,jtest),itest,jtest:', &
363  oro(itest,jtest),itest,jtest
364  print *,'- ELVMAX(',itest,jtest,')=',elvmax(itest,jtest)
365 
366  tend=timef()
367  print*,"- TIMING: FINAL QUALITY CONTROL ", tend-tbeg
368 
369  allocate(xlat(jm), xlon(im))
370  do j = 1, jm
371  xlat(j) = geolat(1,j)
372  enddo
373  do i = 1, im
374  xlon(i) = geolon(i,1)
375  enddo
376 
377  tbeg=timef()
378  CALL write_netcdf(im,jm,slm,land_frac,oro,hprime,1,1, &
379  geolon(1:im,1:jm),geolat(1:im,1:jm), xlon,xlat)
380  tend=timef()
381  print*,"- TIMING: WRITE OUTPUT FILE ", tend-tbeg
382 
383  deallocate(xlat,xlon)
384  deallocate (geolon,geolon_c,geolat,geolat_c)
385  deallocate (slm,oro,var,land_frac)
386  deallocate (theta,gamma,sigma,elvmax,hprime)
387 
388  tend=timef()
389  print*,"- TIMING: TOTAL RUNTIME ", tend-tbeg1
390 
391  return
392  END SUBROUTINE tersub
393 
407  SUBROUTINE make_mask(zslm,slm,land_frac, &
408  im,jm,imn,jmn,lon_c,lat_c)
409 
410  use orog_utils, only : inside_a_polygon, get_index
411 
412  implicit none
413 
414  integer, intent(in) :: zslm(imn,jmn)
415  integer, intent(in) :: im, jm, imn, jmn
416 
417  real, intent(in) :: lon_c(im+1,jm+1), lat_c(im+1,jm+1)
418 
419  real, intent(out) :: slm(im,jm)
420  real, intent(out) :: land_frac(im,jm)
421 
422  real, parameter :: D2R = 3.14159265358979/180.
423 
424  integer jst, jen
425  real GLAT(JMN), GLON(IMN)
426  real LONO(4),LATO(4),LONI,LATI
427  real LONO_RAD(4), LATO_RAD(4)
428  integer JM1,i,j,ii,jj,numx,i2
429  integer ilist(IMN)
430  real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1
431  real XNSUM_ALL,XLAND_ALL,XWATR_ALL
432 
433  print *,'- CREATE LANDMASK AND LAND FRACTION.'
434 !---- GLOBAL XLAT AND XLON ( DEGREE )
435 
436  jm1 = jm - 1
437  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
438 
439  DO j=1,jmn
440  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
441  ENDDO
442  DO i=1,imn
443  glon(i) = 0. + (i-1) * delxn + delxn * 0.5
444  ENDDO
445 
446  land_frac(:,:) = 0.0
447 !
448 !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
449 !
450 ! (*j*) for hard wired zero offset (lambda s =0) for terr05
451 !$omp parallel do &
452 !$omp private (j,i,xnsum,xland,xwatr,xl1,xs1,xw1,lono, &
453 !$omp lato,lono_rad,lato_rad,jst,jen,ilist,numx,jj,i2,ii,loni,lati, &
454 !$omp xnsum_all,xland_all,xwatr_all)
455 !
456  DO j=1,jm
457  DO i=1,im
458  xnsum = 0.0
459  xland = 0.0
460  xwatr = 0.0
461  xnsum_all = 0.0
462  xland_all = 0.0
463  xwatr_all = 0.0
464 
465  lono(1) = lon_c(i,j)
466  lono(2) = lon_c(i+1,j)
467  lono(3) = lon_c(i+1,j+1)
468  lono(4) = lon_c(i,j+1)
469  lato(1) = lat_c(i,j)
470  lato(2) = lat_c(i+1,j)
471  lato(3) = lat_c(i+1,j+1)
472  lato(4) = lat_c(i,j+1)
473  lono_rad=lono*d2r
474  lato_rad=lato*d2r
475  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
476  do jj = jst, jen; do i2 = 1, numx
477  ii = ilist(i2)
478  loni = ii*delxn
479  lati = -90 + jj*delxn
480 
481  xland_all = xland_all + float(zslm(ii,jj))
482  xwatr_all = xwatr_all + float(1-zslm(ii,jj))
483  xnsum_all = xnsum_all + 1.
484 
485  if(inside_a_polygon(loni*d2r,lati*d2r,4, &
486  lono_rad,lato_rad))then
487 
488  xland = xland + float(zslm(ii,jj))
489  xwatr = xwatr + float(1-zslm(ii,jj))
490  xnsum = xnsum + 1.
491  endif
492  enddo ; enddo
493 
494 
495  IF(xnsum.GT.1.) THEN
496  land_frac(i,j) = xland/xnsum
497  slm(i,j) = float(nint(xland/xnsum))
498  ELSEIF(xnsum_all.GT.1.) THEN
499  land_frac(i,j) = xland_all/xnsum_all
500  slm(i,j) = float(nint(xland_all/xnsum_all))
501  ELSE
502  print*, "FATAL ERROR: no source points in MAKE_MASK."
503  call abort()
504  ENDIF
505  ENDDO
506  ENDDO
507 !$omp end parallel do
508 
509  RETURN
510  END SUBROUTINE make_mask
529  SUBROUTINE makemt2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, &
530  IM,JM,IMN,JMN,lon_c,lat_c,lake_frac,land_frac)
531 
532  use orog_utils, only : inside_a_polygon, get_index
533 
534  implicit none
535 
536  integer, intent(in) :: zavg(imn,jmn),zslm(imn,jmn)
537  integer, intent(in) :: im, jm, imn, jmn
538 
539  real, intent(in) :: slm(im,jm)
540  real, intent(in) :: lake_frac(im,jm),land_frac(im,jm)
541  real, intent(in) :: lon_c(im+1,jm+1), lat_c(im+1,jm+1)
542 
543  real, intent(out) :: oro(im,jm)
544  real, intent(out) :: var(im,jm),var4(im,jm)
545 
546  real, parameter :: D2R = 3.14159265358979/180.
547 
548  real, dimension(:), allocatable :: hgt_1d, hgt_1d_all
549 
550  real GLAT(JMN), GLON(IMN)
551  integer JST, JEN, maxsum
552  real LONO(4),LATO(4),LONI,LATI
553  real LONO_RAD(4), LATO_RAD(4)
554  real HEIGHT
555  integer JM1,i,j,nsum,nsum_all,ii,jj,i1,numx,i2
556  integer ilist(IMN)
557  real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1,XW2,XW4
558  real XNSUM_ALL,XLAND_ALL,XWATR_ALL,HEIGHT_ALL
559  real XL1_ALL,XS1_ALL,XW1_ALL,XW2_ALL,XW4_ALL
560 
561  print*,'- CREATE OROGRAPHY AND CONVEXITY.'
562 !---- GLOBAL XLAT AND XLON ( DEGREE )
563 !
564  jm1 = jm - 1
565  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
566 
567  DO j=1,jmn
568  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
569  ENDDO
570  DO i=1,imn
571  glon(i) = 0. + (i-1) * delxn + delxn * 0.5
572  ENDDO
573 
574  maxsum=0
575  DO j=1,jm
576  DO i=1,im
577  lono(1) = lon_c(i,j)
578  lono(2) = lon_c(i+1,j)
579  lono(3) = lon_c(i+1,j+1)
580  lono(4) = lon_c(i,j+1)
581  lato(1) = lat_c(i,j)
582  lato(2) = lat_c(i+1,j)
583  lato(3) = lat_c(i+1,j+1)
584  lato(4) = lat_c(i,j+1)
585  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
586  maxsum=max(maxsum,(jen-jst+1)*imn)
587  ENDDO
588  ENDDO
589  print*,"- MAXSUM IS ", maxsum
590  allocate(hgt_1d(maxsum))
591  allocate(hgt_1d_all(maxsum))
592 !
593 !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
594 !
595 ! (*j*) for hard wired zero offset (lambda s =0) for terr05
596 !$omp parallel do &
597 !$omp private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,xw2,xw4,lono, &
598 !$omp lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,height, &
599 !$omp lato_rad,lono_rad,hgt_1d, &
600 !$omp xnsum_all,xland_all,xwatr_all,nsum_all, &
601 !$omp xl1_all,xs1_all,xw1_all,xw2_all,xw4_all, &
602 !$omp height_all,hgt_1d_all)
603  DO j=1,jm
604  DO i=1,im
605  oro(i,j) = 0.0
606  var(i,j) = 0.0
607  var4(i,j) = 0.0
608  xnsum = 0.0
609  xland = 0.0
610  xwatr = 0.0
611  nsum = 0
612  xl1 = 0.0
613  xs1 = 0.0
614  xw1 = 0.0
615  xw2 = 0.0
616  xw4 = 0.0
617  xnsum_all = 0.0
618  xland_all = 0.0
619  xwatr_all = 0.0
620  nsum_all = 0
621  xl1_all = 0.0
622  xs1_all = 0.0
623  xw1_all = 0.0
624  xw2_all = 0.0
625  xw4_all = 0.0
626 
627  lono(1) = lon_c(i,j)
628  lono(2) = lon_c(i+1,j)
629  lono(3) = lon_c(i+1,j+1)
630  lono(4) = lon_c(i,j+1)
631  lato(1) = lat_c(i,j)
632  lato(2) = lat_c(i+1,j)
633  lato(3) = lat_c(i+1,j+1)
634  lato(4) = lat_c(i,j+1)
635  lono_rad = lono*d2r
636  lato_rad = lato*d2r
637  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
638  do jj = jst, jen; do i2 = 1, numx
639  ii = ilist(i2)
640  loni = ii*delxn
641  lati = -90 + jj*delxn
642 
643  xland_all = xland_all + float(zslm(ii,jj))
644  xwatr_all = xwatr_all + float(1-zslm(ii,jj))
645  xnsum_all = xnsum_all + 1.
646  height_all = float(zavg(ii,jj))
647  nsum_all = nsum_all+1
648  if(nsum_all > maxsum) then
649  print*, "FATAL ERROR: nsum_all is greater than MAXSUM,"
650  print*, "increase MAXSUM.", jst,jen
651  call abort()
652  endif
653  hgt_1d_all(nsum_all) = height_all
654  IF(height_all.LT.-990.) height_all = 0.0
655  xl1_all = xl1_all + height_all * float(zslm(ii,jj))
656  xs1_all = xs1_all + height_all * float(1-zslm(ii,jj))
657  xw1_all = xw1_all + height_all
658  xw2_all = xw2_all + height_all ** 2
659 
660  if(inside_a_polygon(loni*d2r,lati*d2r,4,lono_rad,lato_rad))then
661 
662  xland = xland + float(zslm(ii,jj))
663  xwatr = xwatr + float(1-zslm(ii,jj))
664  xnsum = xnsum + 1.
665  height = float(zavg(ii,jj))
666  nsum = nsum+1
667  if(nsum > maxsum) then
668  print*, "FATAL ERROR: nsum is greater than MAXSUM,"
669  print*, "increase MAXSUM.", jst,jen
670  call abort()
671  endif
672  hgt_1d(nsum) = height
673  IF(height.LT.-990.) height = 0.0
674  xl1 = xl1 + height * float(zslm(ii,jj))
675  xs1 = xs1 + height * float(1-zslm(ii,jj))
676  xw1 = xw1 + height
677  xw2 = xw2 + height ** 2
678  endif
679  enddo ; enddo
680 
681  IF(xnsum.GT.1.) THEN
682  IF(slm(i,j) .NE. 0. .OR. land_frac(i,j) > 0.) THEN
683  IF (xland > 0) THEN
684  oro(i,j)= xl1 / xland
685  ELSE
686  oro(i,j)= xs1 / xwatr
687  ENDIF
688  ELSE
689  IF (xwatr > 0) THEN
690  oro(i,j)= xs1 / xwatr
691  ELSE
692  oro(i,j)= xl1 / xland
693  ENDIF
694  ENDIF
695 
696  var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
697  do i1 = 1, nsum
698  xw4 = xw4 + (hgt_1d(i1) - oro(i,j)) ** 4
699  enddo
700 
701  IF(var(i,j).GT.1.) THEN
702  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
703  ENDIF
704 
705  ELSEIF(xnsum_all.GT.1.) THEN
706 
707  !IF(SLM(I,J).NE.0.) THEN
708  IF(slm(i,j) .NE. 0. .OR. land_frac(i,j) > 0.) THEN
709  IF (xland_all > 0) THEN
710  oro(i,j)= xl1_all / xland_all
711  ELSE
712  oro(i,j)= xs1_all / xwatr_all
713  ENDIF
714  ELSE
715  IF (xwatr_all > 0) THEN
716  oro(i,j)= xs1_all / xwatr_all
717  ELSE
718  oro(i,j)= xl1_all / xland_all
719  ENDIF
720  ENDIF
721 
722  var(i,j)=sqrt(max(xw2_all/xnsum_all-(xw1_all/xnsum_all)**2,0.))
723  do i1 = 1, nsum_all
724  xw4_all = xw4_all + (hgt_1d_all(i1) - oro(i,j)) ** 4
725  enddo
726 
727  IF(var(i,j).GT.1.) THEN
728  var4(i,j) = min(xw4_all/xnsum_all/var(i,j) **4,10.)
729  ENDIF
730  ELSE
731  print*, "FATAL ERROR: no source points in MAKEMT2."
732  call abort()
733  ENDIF
734 
735 ! set orog to 0 meters at ocean.
736 ! IF (LAKE_FRAC(I,J) .EQ. 0. .AND. SLM(I,J) .EQ. 0.)THEN
737  IF (lake_frac(i,j) .EQ. 0. .AND. land_frac(i,j) .EQ. 0.)THEN
738  oro(i,j) = 0.0
739  ENDIF
740 
741  ENDDO
742  ENDDO
743 !$omp end parallel do
744 
745  deallocate(hgt_1d)
746  deallocate(hgt_1d_all)
747  RETURN
748  END SUBROUTINE makemt2
749 
768  SUBROUTINE makepc2(ZAVG,ZSLM,THETA,GAMMA,SIGMA, &
769  IM,JM,IMN,JMN,lon_c,lat_c,SLM)
770 !
771 !=== PC: principal coordinates of each Z avg orog box for L&M
772 !
773  use orog_utils, only : get_index, inside_a_polygon
774 
775  implicit none
776 
777  integer, intent(in) :: zavg(imn,jmn),zslm(imn,jmn)
778  integer, intent(in) :: im,jm,imn,jmn
779 
780  real, intent(in) :: lon_c(im+1,jm+1), lat_c(im+1,jm+1)
781  real, intent(in) :: slm(im,jm)
782 
783  real, intent(out) :: theta(im,jm), gamma(im,jm), sigma(im,jm)
784 
785  real, parameter :: REARTH=6.3712e+6
786  real, parameter :: D2R = 3.14159265358979/180.
787 
788  real GLAT(JMN),DELTAX(JMN)
789  real HL(IM,JM),HK(IM,JM)
790  real HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
791  real SIGMA2(IM,JM)
792  real PI,CERTH,DELXN,DELTAY,XNSUM,XLAND
793  real xfp,yfp,xfpyfp,xfp2,yfp2
794  real hi0,hip1,hj0,hjp1,hijax,hi1j1
795  real LONO(4),LATO(4),LONI,LATI
796  real LONO_RAD(4), LATO_RAD(4)
797  integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
798  integer ilist(IMN)
799  LOGICAL DEBUG
800 !=== DATA DEBUG/.TRUE./
801  DATA debug/.false./
802 
803  print*,"- CREATE PRINCIPLE COORDINATES."
804  pi = 4.0 * atan(1.0)
805  certh = pi * rearth
806 !---- GLOBAL XLAT AND XLON ( DEGREE )
807 !
808  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
809  deltay = certh / float(jmn)
810 
811  DO j=1,jmn
812  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
813  deltax(j) = deltay * cos(glat(j)*d2r)
814  ENDDO
815 !
816 !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
817 !
818 
819 !... DERIVITIVE TENSOR OF HEIGHT
820 !
821 !$omp parallel do &
822 !$omp private (j,i,xnsum,xland,xfp,yfp,xfpyfp, &
823 !$omp xfp2,yfp2,lono,lato,jst,jen,ilist,numx,j1,i2,i1, &
824 !$omp loni,lati,i0,ip1,hi0,hip1,hj0,hjp1,ijax, &
825 !$omp hijax,hi1j1,lono_rad,lato_rad)
826  jloop : DO j=1,jm
827  iloop : DO i=1,im
828  hx2(i,j) = 0.0
829  hy2(i,j) = 0.0
830  hxy(i,j) = 0.0
831  xnsum = 0.0
832  xland = 0.0
833  xfp = 0.0
834  yfp = 0.0
835  xfpyfp = 0.0
836  xfp2 = 0.0
837  yfp2 = 0.0
838  hl(i,j) = 0.0
839  hk(i,j) = 0.0
840  hlprim(i,j) = 0.0
841  theta(i,j) = 0.0
842  gamma(i,j) = 0.
843  sigma2(i,j) = 0.
844  sigma(i,j) = 0.
845 
846  lono(1) = lon_c(i,j)
847  lono(2) = lon_c(i+1,j)
848  lono(3) = lon_c(i+1,j+1)
849  lono(4) = lon_c(i,j+1)
850  lato(1) = lat_c(i,j)
851  lato(2) = lat_c(i+1,j)
852  lato(3) = lat_c(i+1,j+1)
853  lato(4) = lat_c(i,j+1)
854  lato_rad = lato *d2r
855  lono_rad = lono *d2r
856  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
857 
858  do j1 = jst, jen; do i2 = 1, numx
859  i1 = ilist(i2)
860  loni = i1*delxn
861  lati = -90 + j1*delxn
862  inside : if(inside_a_polygon(loni*d2r,lati*d2r,4, &
863  lono_rad,lato_rad))then
864 
865 !=== set the rest of the indexs for ave: 2pt staggered derivitive
866 !
867  i0 = i1 - 1
868  if (i1 - 1 .le. 0 ) i0 = i0 + imn
869  if (i1 - 1 .gt. imn) i0 = i0 - imn
870 
871  ip1 = i1 + 1
872  if (i1 + 1 .le. 0 ) ip1 = ip1 + imn
873  if (i1 + 1 .gt. imn) ip1 = ip1 - imn
874 
875  xland = xland + float(zslm(i1,j1))
876  xnsum = xnsum + 1.
877 
878  hi0 = float(zavg(i0,j1))
879  hip1 = float(zavg(ip1,j1))
880 
881  if(hi0 .lt. -990.) hi0 = 0.0
882  if(hip1 .lt. -990.) hip1 = 0.0
883 !........ xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
884  xfp = 0.5 * ( hip1 - hi0 ) / deltax(j1)
885  xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/deltax(j1) )** 2
886 
887 ! --- not at boundaries
888 !RAB if ( J1 .ne. JST(1) .and. J1 .ne. JEN(JM) ) then
889  if ( j1 .ne. 1 .and. j1 .ne. jmn ) then
890  hj0 = float(zavg(i1,j1-1))
891  hjp1 = float(zavg(i1,j1+1))
892  if(hj0 .lt. -990.) hj0 = 0.0
893  if(hjp1 .lt. -990.) hjp1 = 0.0
894 !....... yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
895  yfp = 0.5 * ( hjp1 - hj0 ) / deltay
896  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/deltay )**2
897 !
898 !..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
899 ! === the NH pole: NB J1 goes from High at NP to Low toward SP
900 !
901 !RAB elseif ( J1 .eq. JST(1) ) then
902  elseif ( j1 .eq. 1 ) then
903  ijax = i1 + imn/2
904  if (ijax .le. 0 ) ijax = ijax + imn
905  if (ijax .gt. imn) ijax = ijax - imn
906 !..... at N pole we stay at the same latitude j1 but cross to opp side
907  hijax = float(zavg(ijax,j1))
908  hi1j1 = float(zavg(i1,j1))
909  if(hijax .lt. -990.) hijax = 0.0
910  if(hi1j1 .lt. -990.) hi1j1 = 0.0
911 !....... yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
912  yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/deltay
913  yfp2 = yfp2 + 0.25 * ( ( 0.5 * ( hijax - hi1j1) ) / deltay )**2
914 !
915 ! === the SH pole: NB J1 goes from High at NP to Low toward SP
916 !
917  elseif ( j1 .eq. jmn ) then
918  ijax = i1 + imn/2
919  if (ijax .le. 0 ) ijax = ijax + imn
920  if (ijax .gt. imn) ijax = ijax - imn
921  hijax = float(zavg(ijax,j1))
922  hi1j1 = float(zavg(i1,j1))
923  if(hijax .lt. -990.) hijax = 0.0
924  if(hi1j1 .lt. -990.) hi1j1 = 0.0
925 !..... yfp = yfp + 0.5 * (0.5 * ( hijax - hi1j1) )/DELTAY
926  yfp = 0.5 * (0.5 * ( hijax - hi1j1) )/deltay
927  yfp2 = yfp2 + 0.25 * ( (0.5 * (hijax - hi1j1) ) / deltay )**2
928  endif
929 !
930 ! === The above does an average across the pole for the bndry in j.
931 !
932  xfpyfp = xfpyfp + xfp * yfp
933  ENDIF inside
934 !
935 ! === average the HX2, HY2 and HXY
936 ! === This will be done over all land
937 !
938  ENDDO
939  ENDDO
940 !
941 ! === HTENSR
942 !
943  xnsum_gt_1 : IF(xnsum.GT.1.) THEN
944  IF(slm(i,j).NE.0.) THEN
945  IF (xland > 0) THEN
946  hx2(i,j) = xfp2 / xland
947  hy2(i,j) = yfp2 / xland
948  hxy(i,j) = xfpyfp / xland
949  ELSE
950  hx2(i,j) = xfp2 / xnsum
951  hy2(i,j) = yfp2 / xnsum
952  hxy(i,j) = xfpyfp / xnsum
953  ENDIF
954  ENDIF
955 !=== degub testing
956  if (debug) then
957  print *," I,J,i1,j1:", i,j,i1,j1,xland,slm(i,j)
958  print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
959  print *," HX2,HY2,HXY:",hx2(i,j),hy2(i,j),hxy(i,j)
960  ENDIF
961 !
962 ! === make the principal axes, theta, and the degree of anisotropy,
963 ! === and sigma2, the slope parameter
964 !
965  hk(i,j) = 0.5 * ( hx2(i,j) + hy2(i,j) )
966  hl(i,j) = 0.5 * ( hx2(i,j) - hy2(i,j) )
967  hlprim(i,j) = sqrt(hl(i,j)*hl(i,j) + hxy(i,j)*hxy(i,j))
968  IF( hl(i,j).NE. 0. .AND. slm(i,j) .NE. 0. ) THEN
969 
970  theta(i,j) = 0.5 * atan2(hxy(i,j),hl(i,j)) / d2r
971 ! === for testing print out in degrees
972 ! THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
973  ENDIF
974  sigma2(i,j) = ( hk(i,j) + hlprim(i,j) )
975  if ( sigma2(i,j) .GE. 0. ) then
976  sigma(i,j) = sqrt(sigma2(i,j) )
977  if (sigma2(i,j) .ne. 0. .and. &
978  hk(i,j) .GE. hlprim(i,j) ) &
979  gamma(i,j) = sqrt( (hk(i,j) - hlprim(i,j)) / sigma2(i,j) )
980  else
981  sigma(i,j)=0.
982  endif
983  ENDIF xnsum_gt_1
984  if (debug) then
985  print *," I,J,THETA,SIGMA,GAMMA,",i,j,theta(i,j),sigma(i,j),gamma(i,j)
986  print *," HK,HL,HLPRIM:",hk(i,j),hl(i,j),hlprim(i,j)
987  endif
988  ENDDO iloop
989  ENDDO jloop
990 !$omp end parallel do
991 
992  RETURN
993  END SUBROUTINE makepc2
994 
1023  SUBROUTINE makeoa2(ZAVG,zslm,VAR,OA4,OL,ELVMAX,ORO,&
1024  IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,dx,dy, &
1025  is_south_pole,is_north_pole )
1026 
1027  use orog_utils, only : get_lat_angle, get_lon_angle, &
1028  get_index, inside_a_polygon, &
1029  get_xnsum, get_xnsum2, &
1030  get_xnsum3
1031 
1032  implicit none
1033 
1034  integer, intent(in) :: im,jm,imn,jmn
1035  integer, intent(in) :: zavg(imn,jmn),zslm(imn,jmn)
1036 
1037  logical, intent(in) :: is_south_pole(im,jm), is_north_pole(im,jm)
1038 
1039  real, intent(in) :: dx(im,jm), dy(im,jm)
1040  real, intent(in) :: lon_c(im+1,jm+1), lat_c(im+1,jm+1)
1041  real, intent(in) :: lon_t(im,jm), lat_t(im,jm)
1042  real, intent(in) :: oro(im,jm), var(im,jm)
1043 
1044  real, intent(out) :: ol(im,jm,4),oa4(im,jm,4)
1045  real, intent(out) :: elvmax(im,jm)
1046 
1047  real, parameter :: MISSING_VALUE = -9999.
1048  real, parameter :: D2R = 3.14159265358979/180.
1049 
1050  integer :: i,j,ilist(imn),numx,i1,j1,ii1
1051  integer :: jst, jen, kwd
1052 
1053  real :: glat(jmn)
1054  real :: zmax(im,jm)
1055  real :: lono(4),lato(4),loni,lati
1056  real :: lono_rad(4), lato_rad(4)
1057  real :: delxn,hc,height,xnpu,xnpd,t
1058  real :: lon,lat,dlon,dlat,dlat_old
1059  real :: lon1,lat1,lon2,lat2
1060  real :: xnsum11,xnsum12,xnsum21,xnsum22
1061  real :: hc_11, hc_12, hc_21, hc_22
1062  real :: xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
1063  real :: xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
1064 
1065  print*,"- CREATE ASYMETRY AND LENGTH SCALE."
1066 !
1067 !---- GLOBAL XLAT AND XLON ( DEGREE )
1068 !
1069  delxn = 360./imn ! MOUNTAIN DATA RESOLUTION
1070 
1071  DO j=1,jmn
1072  glat(j) = -90. + (j-1) * delxn + delxn * 0.5
1073  ENDDO
1074  print*,'- IM=',im,' JM=',jm,' IMN=',imn,' JMN=',jmn
1075 !
1076 !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
1077 !
1078  DO j=1,jm
1079  DO i=1,im
1080  elvmax(i,j) = oro(i,j)
1081  zmax(i,j) = 0.0
1082 !---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
1083 ! IN A GRID BOX
1084  elvmax(i,j) = zmax(i,j)
1085  ENDDO
1086  ENDDO
1087 
1088 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
1089 ! --- to JM or to JM1
1090 !$omp parallel do &
1091 !$omp private (j,i,hc,lono,lato,jst,jen,ilist,numx,j1,ii1,i1,loni, &
1092 !$omp lati,height,lono_rad,lato_rad)
1093  DO j=1,jm
1094  DO i=1,im
1095  hc = 1116.2 - 0.878 * var(i,j)
1096  lono(1) = lon_c(i,j)
1097  lono(2) = lon_c(i+1,j)
1098  lono(3) = lon_c(i+1,j+1)
1099  lono(4) = lon_c(i,j+1)
1100  lato(1) = lat_c(i,j)
1101  lato(2) = lat_c(i+1,j)
1102  lato(3) = lat_c(i+1,j+1)
1103  lato(4) = lat_c(i,j+1)
1104  lono_rad = lono * d2r
1105  lato_rad = lato * d2r
1106  call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
1107  do j1 = jst, jen; do ii1 = 1, numx
1108  i1 = ilist(ii1)
1109  loni = i1*delxn
1110  lati = -90 + j1*delxn
1111  if(inside_a_polygon(loni*d2r,lati*d2r,4, &
1112  lono_rad,lato_rad))then
1113 
1114  height = float(zavg(i1,j1))
1115  IF(height.LT.-990.) height = 0.0
1116  IF ( height .gt. oro(i,j) ) then
1117  if ( height .gt. zmax(i,j) )zmax(i,j) = height
1118  ENDIF
1119  endif
1120  ENDDO ; ENDDO
1121  ENDDO
1122  ENDDO
1123 !$omp end parallel do
1124 
1125  DO j=1,jm
1126  DO i=1,im
1127  elvmax(i,j) = zmax(i,j)
1128  ENDDO
1129  ENDDO
1130 
1131  DO kwd = 1, 4
1132  DO j=1,jm
1133  DO i=1,im
1134  oa4(i,j,kwd) = 0.0
1135  ol(i,j,kwd) = 0.0
1136  ENDDO
1137  ENDDO
1138  ENDDO
1139  !
1140 ! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
1141 !
1142 !---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
1143 !---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
1144 ! (KWD = 1 2 3 4)
1145 ! ( WD = W S SW NW)
1146 
1147 !$omp parallel do &
1148 !$omp private (j,i,lon,lat,kwd,dlon,dlat,lon1,lon2,lat1,lat2, &
1149 !$omp xnsum11,xnsum12,xnsum21,xnsum22,xnpu,xnpd, &
1150 !$omp xnsum1_11,xnsum2_11,hc_11, xnsum1_12,xnsum2_12, &
1151 !$omp hc_12,xnsum1_21,xnsum2_21,hc_21, xnsum1_22, &
1152 !$omp xnsum2_22,hc_22)
1153  DO j=1,jm
1154  DO i=1,im
1155  lon = lon_t(i,j)
1156  lat = lat_t(i,j)
1157  !--- for around north pole, oa and ol are all 0
1158 
1159  if(is_north_pole(i,j)) then
1160  print*, "- SET OA1 = 0 AND OL=0 AT I,J=", i,j
1161  do kwd = 1, 4
1162  oa4(i,j,kwd) = 0.
1163  ol(i,j,kwd) = 0.
1164  enddo
1165  else if(is_south_pole(i,j)) then
1166  print*, "- SET OA1 = 0 AND OL=1 AT I,J=", i,j
1167  do kwd = 1, 4
1168  oa4(i,j,kwd) = 0.
1169  ol(i,j,kwd) = 1.
1170  enddo
1171  else
1172 
1173  !--- for each point, find a lat-lon grid box with same dx and dy as the cubic grid box
1174  dlon = get_lon_angle(dx(i,j), lat )
1175  dlat = get_lat_angle(dy(i,j))
1176  !--- adjust dlat if the points are close to pole.
1177  if( lat-dlat*0.5<-90.) then
1178  print*, "- AT I,J =", i,j, lat, dlat, lat-dlat*0.5
1179  print*, "FATAL ERROR: lat-dlat*0.5<-90."
1180  call errexit(4)
1181  endif
1182  if( lat+dlat*2 > 90.) then
1183  dlat_old = dlat
1184  dlat = (90-lat)*0.5
1185  print*, "- AT I,J=",i,j," ADJUST DLAT FROM ", &
1186  dlat_old, " TO ", dlat
1187  endif
1188  !--- lower left
1189  lon1 = lon-dlon*1.5
1190  lon2 = lon-dlon*0.5
1191  lat1 = lat-dlat*0.5
1192  lat2 = lat+dlat*0.5
1193 
1194  if(lat1<-90 .or. lat2>90) then
1195  print*, "- AT UPPER LEFT I=,J=", i, j, lat, dlat,lat1,lat2
1196  endif
1197  xnsum11 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1198  zavg,zslm,delxn)
1199 
1200  !--- upper left
1201  lon1 = lon-dlon*1.5
1202  lon2 = lon-dlon*0.5
1203  lat1 = lat+dlat*0.5
1204  lat2 = lat+dlat*1.5
1205  if(lat1<-90 .or. lat2>90) then
1206  print*, "- AT LOWER LEFT I=,J=", i, j, lat, dlat,lat1,lat2
1207  endif
1208  xnsum12 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1209  zavg,zslm,delxn)
1210 
1211  !--- lower right
1212  lon1 = lon-dlon*0.5
1213  lon2 = lon+dlon*0.5
1214  lat1 = lat-dlat*0.5
1215  lat2 = lat+dlat*0.5
1216  if(lat1<-90 .or. lat2>90) then
1217  print*, "- AT UPPER RIGHT I=,J=", i, j, lat, dlat,lat1,lat2
1218  endif
1219  xnsum21 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1220  zavg,zslm,delxn)
1221 
1222  !--- upper right
1223  lon1 = lon-dlon*0.5
1224  lon2 = lon+dlon*0.5
1225  lat1 = lat+dlat*0.5
1226  lat2 = lat+dlat*1.5
1227  if(lat1<-90 .or. lat2>90) then
1228  print*, "- AT LOWER RIGHT I=,J=", i, j, lat, dlat,lat1,lat2
1229  endif
1230 
1231  xnsum22 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1232  zavg,zslm,delxn)
1233 
1234  xnpu = xnsum11 + xnsum12
1235  xnpd = xnsum21 + xnsum22
1236  IF (xnpd .NE. xnpu) oa4(i,j,1) = 1. - xnpd / max(xnpu , 1.)
1237 
1238  xnpu = xnsum11 + xnsum21
1239  xnpd = xnsum12 + xnsum22
1240  IF (xnpd .NE. xnpu) oa4(i,j,2) = 1. - xnpd / max(xnpu , 1.)
1241 
1242  xnpu = xnsum11 + (xnsum12+xnsum21)*0.5
1243  xnpd = xnsum22 + (xnsum12+xnsum21)*0.5
1244  IF (xnpd .NE. xnpu) oa4(i,j,3) = 1. - xnpd / max(xnpu , 1.)
1245 
1246  xnpu = xnsum12 + (xnsum11+xnsum22)*0.5
1247  xnpd = xnsum21 + (xnsum11+xnsum22)*0.5
1248  IF (xnpd .NE. xnpu) oa4(i,j,4) = 1. - xnpd / max(xnpu , 1.)
1249 
1250 
1251  !--- calculate OL3 and OL4
1252  !--- lower left
1253  lon1 = lon-dlon*1.5
1254  lon2 = lon-dlon*0.5
1255  lat1 = lat-dlat*0.5
1256  lat2 = lat+dlat*0.5
1257  if(lat1<-90 .or. lat2>90) then
1258  print*, "- AT UPPER LEFT I=,J=", i, j, lat, dlat,lat1,lat2
1259  endif
1260  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1261  zavg,delxn, xnsum1_11, xnsum2_11, hc_11)
1262 
1263  !--- upper left
1264  lon1 = lon-dlon*1.5
1265  lon2 = lon-dlon*0.5
1266  lat1 = lat+dlat*0.5
1267  lat2 = lat+dlat*1.5
1268  if(lat1<-90 .or. lat2>90) then
1269  print*, "- AT LOWER LEFT I=,J=", i, j, lat, dlat,lat1,lat2
1270  endif
1271  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1272  zavg,delxn, xnsum1_12, xnsum2_12, hc_12)
1273 
1274  !--- lower right
1275  lon1 = lon-dlon*0.5
1276  lon2 = lon+dlon*0.5
1277  lat1 = lat-dlat*0.5
1278  lat2 = lat+dlat*0.5
1279  if(lat1<-90 .or. lat2>90) then
1280  print*, "- AT UPPER RIGHT I=,J=", i, j, lat, dlat,lat1,lat2
1281  endif
1282  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1283  zavg,delxn, xnsum1_21, xnsum2_21, hc_21)
1284 
1285  !--- upper right
1286  lon1 = lon-dlon*0.5
1287  lon2 = lon+dlon*0.5
1288  lat1 = lat+dlat*0.5
1289  lat2 = lat+dlat*1.5
1290  if(lat1<-90 .or. lat2>90) then
1291  print*, "- AT LOWER RIGHT I=,J=", i, j, lat, dlat,lat1,lat2
1292  endif
1293  call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1294  zavg,delxn, xnsum1_22, xnsum2_22, hc_22)
1295 
1296  ol(i,j,3) = (xnsum1_22+xnsum1_11)/(xnsum2_22+xnsum2_11)
1297  ol(i,j,4) = (xnsum1_12+xnsum1_21)/(xnsum2_12+xnsum2_21)
1298 
1299  !--- calculate OL1 and OL2
1300  !--- lower left
1301  lon1 = lon-dlon*2.0
1302  lon2 = lon-dlon
1303  lat1 = lat
1304  lat2 = lat+dlat
1305  if(lat1<-90 .or. lat2>90) then
1306  print*, "- AT UPPER LEFT I=,J=", i, j, lat, dlat,lat1,lat2
1307  endif
1308  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1309  zavg,delxn, xnsum1_11, xnsum2_11, hc_11)
1310 
1311  !--- upper left
1312  lon1 = lon-dlon*2.0
1313  lon2 = lon-dlon
1314  lat1 = lat+dlat
1315  lat2 = lat+dlat*2.0
1316  if(lat1<-90 .or. lat2>90) then
1317  print*, "- AT LOWER LEFT I=,J=", i, j, lat, dlat,lat1,lat2
1318  endif
1319 
1320  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1321  zavg,delxn, xnsum1_12, xnsum2_12, hc_12)
1322 
1323  !--- lower right
1324  lon1 = lon-dlon
1325  lon2 = lon
1326  lat1 = lat
1327  lat2 = lat+dlat
1328  if(lat1<-90 .or. lat2>90) then
1329  print*, "- AT UPPER RIGHT I=,J=", i, j, lat, dlat,lat1,lat2
1330  endif
1331  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1332  zavg,delxn, xnsum1_21, xnsum2_21, hc_21)
1333 
1334  !--- upper right
1335  lon1 = lon-dlon
1336  lon2 = lon
1337  lat1 = lat+dlat
1338  lat2 = lat+dlat*2.0
1339  if(lat1<-90 .or. lat2>90) then
1340  print*, "- AT LOWER RIGHT I=,J=", i, j, lat, dlat,lat1,lat2
1341  endif
1342 
1343  call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat, &
1344  zavg,delxn, xnsum1_22, xnsum2_22, hc_22)
1345 
1346  ol(i,j,1) = (xnsum1_11+xnsum1_21)/(xnsum2_11+xnsum2_21)
1347  ol(i,j,2) = (xnsum1_21+xnsum1_22)/(xnsum2_21+xnsum2_22)
1348  ENDIF
1349  ENDDO
1350  ENDDO
1351 !$omp end parallel do
1352  DO kwd=1,4
1353  DO j=1,jm
1354  DO i=1,im
1355  t = oa4(i,j,kwd)
1356  oa4(i,j,kwd) = sign( min( abs(t), 1. ), t )
1357  ENDDO
1358  ENDDO
1359  ENDDO
1360 
1361  RETURN
1362 
1363  END SUBROUTINE makeoa2
orog_utils::minmax
subroutine, public minmax(im, jm, a, title, imax, jmax)
Print out the maximum and minimum values of an array and optionally pass back the i/j location of the...
Definition: orog_utils.F90:50
tersub
subroutine tersub(IM, JM, EFAC, OUTGRID, MASK_ONLY, EXTERNAL_MASK_FILE)
Driver routine to compute terrain.
Definition: mtnlm7_oclsm.F90:119
makemt2
subroutine makemt2(ZAVG, ZSLM, ORO, SLM, VAR, VAR4, IM, JM, IMN, JMN, lon_c, lat_c, lake_frac, land_frac)
Create the orography, standard deviation of orography and the convexity on a model tile...
Definition: mtnlm7_oclsm.F90:531
orog_utils::get_xnsum2
subroutine, public 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: orog_utils.F90:904
makeoa2
subroutine makeoa2(ZAVG, zslm, VAR, OA4, OL, ELVMAX, ORO, 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.F90:1026
io_utils::qc_orog_by_ramp
subroutine, public qc_orog_by_ramp(imn, jmn, zavg, zslm)
Quality control the global orography and landmask data over Antarctica using RAMP data...
Definition: io_utils.F90:669
orog_utils::get_xnsum
real function, public 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: orog_utils.F90:793
io_utils::read_mdl_grid_file
subroutine, public read_mdl_grid_file(mdl_grid_file, im, jm, geolon, geolon_c, geolat, geolat_c, dx, dy, is_north_pole, is_south_pole)
Read the grid dimensions from the model &#39;grid&#39; file.
Definition: io_utils.F90:454
io_utils
Module containing utilities that read and write data.
Definition: io_utils.F90:9
orog_utils::get_index
subroutine, public 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: orog_utils.F90:691
orog_utils::inside_a_polygon
logical function, public inside_a_polygon(lon1, lat1, npts, lon2, lat2)
Check if a point is inside a polygon.
Definition: orog_utils.F90:318
io_utils::read_global_orog
subroutine, public read_global_orog(imn, jmn, glob)
Read input global 30-arc second orography data.
Definition: io_utils.F90:552
io_utils::read_mdl_dims
subroutine, public read_mdl_dims(mdl_grid_file, im, jm)
Read the grid dimensions from the model &#39;grid&#39; file.
Definition: io_utils.F90:402
io_utils::read_global_mask
subroutine, public read_global_mask(imn, jmn, mask)
Read input global 30-arc second land mask data.
Definition: io_utils.F90:611
orog_utils::get_lat_angle
real function, public get_lat_angle(dy)
Convert the &#39;y&#39; direction distance of a cubed-sphere grid point to the corresponding distance in lati...
Definition: orog_utils.F90:128
makepc2
subroutine makepc2(ZAVG, ZSLM, THETA, GAMMA, SIGMA, 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.F90:770
orog_utils::remove_isolated_pts
subroutine, public remove_isolated_pts(im, jm, slm, oro, var, var4, oa, ol)
Remove isolated model points.
Definition: orog_utils.F90:561
orog_utils
Module containing utilites used by the orog program.
Definition: orog_utils.F90:8
orog_utils::get_xnsum3
subroutine, public 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: orog_utils.F90:1006
make_mask
subroutine make_mask(zslm, slm, land_frac, im, jm, imn, jmn, lon_c, lat_c)
Create the land-mask, land fraction.
Definition: mtnlm7_oclsm.F90:409
io_utils::write_netcdf
subroutine, public write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolon, geolat, lon, lat)
Write out orography file in netcdf format.
Definition: io_utils.F90:42
orog_utils::timef
real function, public timef()
Get the date/time from the system clock.
Definition: orog_utils.F90:1064
orog_utils::get_lon_angle
real function, public get_lon_angle(dx, lat_in)
Convert the &#39;x&#39; direction distance of a cubed-sphere grid point to the corresponding distance in long...
Definition: orog_utils.F90:152
io_utils::read_mask
subroutine, public read_mask(merge_file, slm, land_frac, lake_frac, im, jm)
Read the land mask file.
Definition: io_utils.F90:357
io_utils::write_mask_netcdf
subroutine, public write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geolat)
Write the land mask file.
Definition: io_utils.F90:267