Changeset 226 for trunk/SRC/ToBeReviewed/GRILLE/domdef.pro
- Timestamp:
- 03/16/07 10:22:26 (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SRC/ToBeReviewed/GRILLE/domdef.pro
r163 r226 4 4 ;+ 5 5 ; 6 ; @file_comments 7 ; Allows to extract a sub-domain of study by providing parameters 6 ; @file_comments 7 ; Allows to extract a sub-domain of study by providing parameters 8 8 ; needed for drawings (see outputs) 9 9 ; … … 12 12 ; @param Z1 {in}{optional} 13 13 ; For a 3d domain whose the horizontal part cover all glam 14 ; 14 ; 15 15 ; @param Z2 {in}{optional} 16 16 ; For a 3d domain whose the horizontal part cover all gphi … … 38 38 ; 39 39 ; @keyword FINDALWAYS 40 ; Force to redefine a box even when none point is find in the box. 40 ; Force to redefine a box even when none point is find in the box. 41 41 ; In this case, we select all the grid. 42 42 ; 43 43 ; @keyword GRIDTYPE {type=string or vector} 44 ; It is a string or a vector of strings containing the grids's name 45 ; (Only determined by 'T','U','V','W','F') for which the calculation 46 ; must be done. 44 ; It is a string or a vector of strings containing the grids's name 45 ; (Only determined by 'T','U','V','W','F') for which the calculation 46 ; must be done. 47 47 ; For example: 'T' or ['T','U'] 48 48 ; … … 55 55 ; grids. 56 56 ; 57 ; @keyword INDEX 58 ; We activate it if we want that all elements passed in input of domdef 59 ; refer to indexes of glam, gphi and gdep arrays rather than to values 57 ; @keyword INDEX 58 ; We activate it if we want that all elements passed in input of domdef 59 ; refer to indexes of glam, gphi and gdep arrays rather than to values 60 60 ; of these arrays. 61 61 ; 62 ; @keyword XINDEX 63 ; We activate it if we want that all elements passed in input of domdef 64 ; and concerning the X dimension refer to indexes of glam arrays rather 62 ; @keyword XINDEX 63 ; We activate it if we want that all elements passed in input of domdef 64 ; and concerning the X dimension refer to indexes of glam arrays rather 65 65 ; than to values of these arrays. 66 66 ; 67 ; @keyword YINDEX 68 ; We activate it if we want that all elements passed in input of domdef 69 ; and concerning the X dimension refer to indexes of gphi arrays rather 67 ; @keyword YINDEX 68 ; We activate it if we want that all elements passed in input of domdef 69 ; and concerning the X dimension refer to indexes of gphi arrays rather 70 70 ; than to values of these arrays. 71 ; 72 ; @keyword ZINDEX 73 ; We activate it if we want that all elements passed in input of domdef 74 ; and concerning the X dimension refer to indexes of gdep arrays rather 71 ; 72 ; @keyword ZINDEX 73 ; We activate it if we want that all elements passed in input of domdef 74 ; and concerning the X dimension refer to indexes of gdep arrays rather 75 75 ; than to values of these arrays. 76 76 ; 77 77 ; @uses 78 78 ; common.pro 79 ;80 79 ; 81 80 ; @history … … 123 122 ; 124 123 IF keyword_set(endpoints) THEN BEGIN 125 IF NOT keyword_set(type) THEN BEGIN 124 IF NOT keyword_set(type) THEN BEGIN 126 125 dummy = report(['If domdef is used do find the box associated' $ 127 126 , 'to endpoints, you must also specify type keyword']) 128 127 return 129 ENDIF 128 ENDIF 130 129 CASE N_PARAMS() OF 131 130 0: … … 137 136 section, BOXZOOM = boxzoom, ENDPOINTS = endpoints, TYPE = type, /ONLYBOX 138 137 return 139 ENDIF 138 ENDIF 140 139 141 140 ;------------------------------------------------------------------- … … 153 152 END 154 153 ENDCASE 155 RETURN 154 RETURN 156 155 ENDIF 157 156 ;------------------------------------------------------------------- … … 170 169 vert1t = 99999. & vert2t = -99999. & vert1w = 99999. & vert2w = -99999. 171 170 ; 172 IF jpj EQ 1 THEN BEGIN 171 IF jpj EQ 1 THEN BEGIN 173 172 IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 174 173 glamt = reform(glamt, jpi, jpj, /over) 175 174 gphit = reform(gphit, jpi, jpj, /over) 176 ENDIF 175 ENDIF 177 176 IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 178 177 glamu = reform(glamu, jpi, jpj, /over) 179 178 gphiu = reform(gphiu, jpi, jpj, /over) 180 ENDIF 179 ENDIF 181 180 IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 182 181 glamv = reform(glamv, jpi, jpj, /over) 183 182 gphiv = reform(gphiv, jpi, jpj, /over) 184 ENDIF 183 ENDIF 185 184 IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 186 185 glamf = reform(glamf, jpi, jpj, /over) 187 186 gphif = reform(gphif, jpi, jpj, /over) 188 ENDIF 187 ENDIF 189 188 ENDIF 190 189 ; … … 229 228 AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN 230 229 IF N_PARAMS() EQ 0 THEN BEGIN 231 ; find lon1 and lon2 the longitudinal bou daries of the full domain230 ; find lon1 and lon2 the longitudinal boundaries of the full domain 232 231 IF (where(gridtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t) 233 232 IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t) … … 237 236 lon1 = min([lon1t, lon1u, lon1v, lon1f]) 238 237 lon2 = max([lon2t, lon2u, lon2v, lon2f]) 239 ; find lat1 and lat2 the latitudinal bou daries of the full domain238 ; find lat1 and lat2 the latitudinal boundaries of the full domain 240 239 IF (where(gridtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t) 241 240 IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t) … … 245 244 lat1 = min([lat1t, lat1u, lat1v, lat1f]) 246 245 lat2 = max([lat2t, lat2u, lat2v, lat2f]) 247 ENDIF ELSE BEGIN 246 ENDIF ELSE BEGIN 248 247 lon1 = min([x1, x2], max = lon2) 249 248 lat1 = min([y1, y2], max = lat2) … … 288 287 IF (dom[0] EQ -1) THEN BEGIN 289 288 IF keyword_set(findalways) THEN BEGIN 290 ; if t grid parameters alreday defined, we use them... 289 ; if t grid parameters alreday defined, we use them... 291 290 CASE 1 OF 292 291 (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1:BEGIN … … 337 336 firstxv = firstxt & lastxv = lastxt & nxv = nxt 338 337 firstyv = firstyt & lastyv = lastyt & nyv = nyt 339 END 338 END 340 339 (where(gridtype eq 'U'))[0] NE -1:BEGIN 341 340 print, 'WARNING, empty V points box... we use the same index as U points...' 342 341 firstxv = firstxu & lastxv = lastxu & nxv = nxu 343 342 firstyv = firstyu & lastyv = lastyu & nyv = nyu 344 END 343 END 345 344 ELSE:BEGIN 346 345 print, 'WARNING, empty V points box... we get the neighnors to define a new box...' … … 385 384 firstxf = firstxt & lastxf = lastxt & nxf = nxt 386 385 firstyf = firstyt & lastyf = lastyt & nyf = nyt 387 END 386 END 388 387 (where(gridtype eq 'U'))[0] NE -1:BEGIN 389 388 print, 'WARNING, empty F points box... we use the same index as U points...' 390 389 firstxf = firstxu & lastxf = lastxu & nxf = nxu 391 390 firstyf = firstyu & lastyf = lastyu & nyf = nyu 392 END 391 END 393 392 (where(gridtype eq 'V'))[0] NE -1:BEGIN 394 393 print, 'WARNING, empty F points box... we use the same index as V points...' 395 394 firstxf = firstxv & lastxf = lastxv & nxf = nxv 396 395 firstyf = firstyv & lastyf = lastyv & nyf = nyv 397 END 396 END 398 397 ELSE:BEGIN 399 398 print, 'WARNING, empty F points box... we get the neighnors to define a new box...' … … 406 405 RETURN 407 406 END 408 ENDCASE 407 ENDCASE 409 408 ENDIF ELSE BEGIN 410 409 ras = report('WARNING, The box does not contain any F points...') … … 469 468 firstyv = fsty & lastyv = lsty 470 469 nxv = nx & nyv = ny 471 ENDIF 470 ENDIF 472 471 ; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf 473 472 ; according to x1, x2, y1, y2 … … 537 536 print, 'WARNING, empty U points box... we use the same index as T points...' 538 537 firstyu = firstyt & lastyu = lastyt & nyu = nyt 539 END 538 END 540 539 ELSE:BEGIN 541 540 print, 'WARNING, empty U points box... we get the neighnors to define a new box...' … … 548 547 RETURN 549 548 END 550 ENDCASE 549 ENDCASE 551 550 ENDIF ELSE BEGIN 552 551 ras = report('WARNING, The box does not contain any U points...') … … 575 574 print, 'WARNING, empty V points box... we use the same index as T points...' 576 575 firstyv = firstyt & lastyv = lastyt & nyv = nyt 577 END 576 END 578 577 (where(gridtype eq 'U'))[0] NE -1:BEGIN 579 578 print, 'WARNING, empty V points box... we use the same index as U points...' 580 579 firstyv = firstyu & lastyv = lastyu & nyv = nyu 581 END 580 END 582 581 ELSE:BEGIN 583 582 print, 'WARNING, empty V points box... we get the neighnors to define a new box...' … … 589 588 ENDCASE 590 589 RETURN 591 END 592 ENDCASE 590 END 591 ENDCASE 593 592 ENDIF ELSE BEGIN 594 593 ras = report('WARNING, The box does not contain any V points...') … … 617 616 print, 'WARNING, empty F points box... we use the same index as T points...' 618 617 firstyf = firstyt & lastyf = lastyt & nyf = nyt 619 END 618 END 620 619 (where(gridtype eq 'U'))[0] NE -1:BEGIN 621 620 print, 'WARNING, empty F points box... we use the same index as U points...' 622 621 firstyf = firstyu & lastyf = lastyu & nyf = nyu 623 END 622 END 624 623 (where(gridtype eq 'V'))[0] NE -1:BEGIN 625 624 print, 'WARNING, empty F points box... we use the same index as V points...' 626 625 firstyf = firstyv & lastyf = lastyv & nyf = nyv 627 END 626 END 628 627 ELSE:BEGIN 629 628 print, 'WARNING, empty F points box... we get the neighnors to define a new box...' … … 636 635 RETURN 637 636 END 638 ENDCASE 637 ENDCASE 639 638 ENDIF ELSE BEGIN 640 639 ras = report('WARNING, The box does not contain any F points...') … … 705 704 print, 'WARNING, empty U points box... we use the same index as T points...' 706 705 firstxu = firstxt & lastxu = lastxt & nxu = nxt 707 END 706 END 708 707 ELSE:BEGIN 709 708 print, 'WARNING, empty U points box... we get the neighnors to define a new box...' … … 716 715 RETURN 717 716 END 718 ENDCASE 717 ENDCASE 719 718 ENDIF ELSE BEGIN 720 719 ras = report('WARNING, The box does not contain any U points...') … … 743 742 print, 'WARNING, empty V points box... we use the same index as T points...' 744 743 firstxv = firstxt & lastxv = lastxt & nxv = nxt 745 END 744 END 746 745 (where(gridtype eq 'U'))[0] NE -1:BEGIN 747 746 print, 'WARNING, empty V points box... we use the same index as U points...' 748 747 firstxv = firstxu & lastxv = lastxu & nxv = nxu 749 END 748 END 750 749 ELSE:BEGIN 751 750 print, 'WARNING, empty V points box... we get the neighnors to define a new box...' … … 757 756 ENDCASE 758 757 RETURN 759 END 760 ENDCASE 758 END 759 ENDCASE 761 760 ENDIF ELSE BEGIN 762 761 ras = report('WARNING, The box does not contain any V points...') … … 785 784 print, 'WARNING, empty F points box... we use the same index as T points...' 786 785 firstxf = firstxt & lastxf = lastxt & nxf = nxt 787 END 786 END 788 787 (where(gridtype eq 'U'))[0] NE -1:BEGIN 789 788 print, 'WARNING, empty F points box... we use the same index as U points...' 790 789 firstxf = firstxu & lastxf = lastxu & nxf = nxu 791 END 790 END 792 791 (where(gridtype eq 'V'))[0] NE -1:BEGIN 793 792 print, 'WARNING, empty F points box... we use the same index as V points...' 794 793 firstxf = firstxv & lastxf = lastxv & nxf = nxv 795 END 794 END 796 795 ELSE:BEGIN 797 796 print, 'WARNING, empty F points box... we get the neighnors to define a new box...' … … 803 802 ENDCASE 804 803 RETURN 805 END 806 ENDCASE 804 END 805 ENDCASE 807 806 ENDIF ELSE BEGIN 808 807 ras = report('WARNING, The box does not contain any F points...') … … 821 820 END 822 821 ENDCASE 823 ENDELSE 824 ;------------------------------------------------------------------- 825 ; The extracted domain is it regular or not? 822 ENDELSE 823 ;------------------------------------------------------------------- 824 ; The extracted domain is it regular or not? 826 825 ;------------------------------------------------------------------- 827 826 CASE 1 OF 828 827 ((where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN 829 ; to get faster, we first test the most basic cases befor testing the828 ; to get faster, we first test the most basic cases before testing the 830 829 ; full array. 831 830 CASE 0 OF … … 908 907 firstzt = domz[0] 909 908 lastzt = domz[nzt-1] 910 ENDIF ELSE BEGIN 909 ENDIF ELSE BEGIN 911 910 ras = report('WARNING, The box does not contain any T level...') 912 911 firstzt = -1 913 912 lastzt = -1 914 ENDELSE 913 ENDELSE 915 914 ENDIF 916 915 ; define firstzw, firstzw, nzw … … 920 919 lastzw = lastzt 921 920 nzw = nzt 922 ENDIF ELSE BEGIN 921 ENDIF ELSE BEGIN 923 922 domz = where(gdepw ge vert1 and gdepw le vert2, nzw) 924 923 IF nzw NE 0 THEN BEGIN 925 924 firstzw = domz[0] 926 925 lastzw = domz[nzw-1] 927 ENDIF ELSE BEGIN 926 ENDIF ELSE BEGIN 928 927 ras = report('WARNING, The box does not contain any W level...') 929 928 firstzw = -1 930 929 lastzw = -1 931 ENDELSE 932 ENDELSE 930 ENDELSE 931 ENDELSE 933 932 ENDIF 934 933 ;------------------------------------------------------------------- 935 934 ; vertical domain defined with the Z index 936 935 ;------------------------------------------------------------------- 937 ENDIF ELSE BEGIN 936 ENDIF ELSE BEGIN 938 937 CASE N_PARAMS() OF 939 938 2:fstz = min([x1, x2], max = lstz) … … 951 950 vert1t = min(gdept[fstz:lstz], max = vert2t) 952 951 firstzt = fstz & lastzt = lstz & nzt = nz 953 ENDIF 952 ENDIF 954 953 ; find vert1w, vert2w, firstzw, firstzw, nzw 955 954 ; according to (x1, x2) or (z1, z2) … … 957 956 vert1w = min(gdepw[fstz:lstz], max = vert2w) 958 957 firstzw = fstz & lastzw = lstz & nzw = nz 959 ENDIF 958 ENDIF 960 959 vert1 = min([vert1t, vert1w]) 961 960 vert2 = max([vert2t, vert2w]) … … 965 964 IF NOT keyword_set(key_forgetold) THEN BEGIN 966 965 @updateold 967 ENDIF 968 ;------------------------------------------------------------------- 969 if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun 966 ENDIF 967 ;------------------------------------------------------------------- 968 if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun 970 969 ;------------------------------------------------------------------- 971 970
Note: See TracChangeset
for help on using the changeset viewer.