orog_mask_tools 1.14.0
Loading...
Searching...
No Matches
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
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
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
tersub
subroutine tersub(im, jm, efac, outgrid, mask_only, external_mask_file)
Driver routine to compute terrain.
Definition mtnlm7_oclsm.F90:119
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
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
Module containing utilities that read and write data.
Definition io_utils.F90:9
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
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_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 'grid' file.
Definition io_utils.F90:454
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
io_utils::read_mdl_dims
subroutine, public read_mdl_dims(mdl_grid_file, im, jm)
Read the grid dimensions from the model 'grid' file.
Definition io_utils.F90:402
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
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
Module containing utilites used by the orog program.
Definition orog_utils.F90:8
orog_utils::get_lon_angle
real function, public get_lon_angle(dx, lat_in)
Convert the 'x' direction distance of a cubed-sphere grid point to the corresponding distance in long...
Definition orog_utils.F90:152
orog_utils::timef
real function, public timef()
Get the date/time from the system clock.
Definition orog_utils.F90:1064
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
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
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::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
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
orog_utils::get_lat_angle
real function, public get_lat_angle(dy)
Convert the 'y' direction distance of a cubed-sphere grid point to the corresponding distance in lati...
Definition orog_utils.F90:128
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
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