Ignore:
Timestamp:
Apr 30, 2008 1:41:13 PM (14 years ago)
Author:
letzel
Message:
  • NCL scripts in trunk/SCRIPTS/NCL updated
File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/NCL/cross_sections.ncl

    r161 r162  
    6161   if ( .not. isvar("sort") ) then                      ; sort of plots
    6262      sort = "time"
    63       if (parameter(37) .NE. "time") then
    64          sort = parameter(37
     63      if (parameter(51) .NE. "time") then
     64         sort = parameter(51
    6565      end if
    6666   end if
    6767   if ( .not. isvar("mode") ) then                      ; pattern of contour plots
    6868      mode = "Fill"
    69       if (parameter(39) .NE. "Fill") then
    70          mode = parameter(39) 
     69      if (parameter(49) .NE. "Fill") then
     70         mode = parameter(49) 
    7171      end if
    7272   end if
    7373   if ( .not. isvar("fill_mode") ) then                 ; pattern of filling
    7474      fill_mode = "AreaFill"
    75       if (parameter(41) .NE. "AreaFill") then
    76          mode = parameter(41
     75      if (parameter(53) .NE. "AreaFill") then
     76         mode = parameter(53
    7777      end if
    7878   end if
    7979   if ( .not. isvar("shape") ) then                     ; keeping of aspect ratio
    8080      shape = 1
    81       if (parameter(43) .NE. "1") then
    82          shape = stringtointeger(parameter(43))
    83          if (stringtointeger(parameter(43)) .NE. 0) then
     81      if (parameter(55) .NE. "1") then
     82         shape = stringtointeger(parameter(55))
     83         if (stringtointeger(parameter(55)) .NE. 0) then
    8484            print(" ")
    8585            print("Please set 'shape' to 0 or 1")
     
    9595   if ( .not. isvar("xyc") ) then                       ; turn xy cross-section on or off
    9696      xyc = 0
    97       if (parameter(45) .NE. "0") then
    98          xyc = stringtointeger(parameter(45))
    99          if (stringtointeger(parameter(45)) .NE. 1) then
     97      if (parameter(57) .NE. "0") then
     98         xyc = stringtointeger(parameter(57))
     99         if (stringtointeger(parameter(57)) .NE. 1) then
    100100            print(" ")
    101101            print("Please set 'xyc' to 0 or 1")
     
    107107   if ( .not. isvar("xzc") ) then                       ; turn xz cross-section on or off
    108108      xzc = 0
    109       if (parameter(47) .NE. "0") then
    110          xzc = stringtointeger(parameter(47))
    111          if (stringtointeger(parameter(47)) .NE. 1) then
     109      if (parameter(59) .NE. "0") then
     110         xzc = stringtointeger(parameter(59))
     111         if (stringtointeger(parameter(59)) .NE. 1) then
    112112            print(" ")
    113113            print("Please set 'xzc' to 0 or 1")
     
    119119   if ( .not. isvar("yzc") ) then                       ; turn yz cross-section on or off
    120120      yzc = 0
    121       if (parameter(49) .NE. "0") then
    122          yzc = stringtointeger(parameter(49))
    123          if (stringtointeger(parameter(49)) .NE. 1) then
     121      if (parameter(61) .NE. "0") then
     122         yzc = stringtointeger(parameter(61))
     123         if (stringtointeger(parameter(61)) .NE. 1) then
    124124            print(" ")
    125125            print("Please set 'yzc' to 0 or 1")
     
    129129      end if
    130130   end if
    131    if ( .not. isvar("xs") ) then                        ; start of x-coordinate range
    132       xs = 0
    133       if (parameter(51) .NE. "0") then
    134          xs = stringtointeger(parameter(51))
    135       end if
    136    end if
    137    if ( .not. isvar("ys") ) then                        ; start of y-coordinate range
    138       ys = 0
    139       if (parameter(55) .NE. "0") then
    140          ys = stringtointeger(parameter(55))
    141       end if
    142    end if
    143    if ( .not. isvar("zs") ) then                        ; start of z-coordinate range
    144       zs = 0
    145       if (parameter(59) .NE. "0") then
    146          zs = stringtointeger(parameter(59))
    147       end if
    148    end if 
    149131
    150132   if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
     
    180162   if ( .not. isvar("vector") ) then                    ; sort of plots
    181163      vector = 0
    182       if (parameter(63) .NE. "0") then
    183          vector = stringtointeger(parameter(63))
    184          if (stringtointeger(parameter(63)) .NE. 1) then
     164      if (parameter(39) .NE. "0") then
     165         vector = stringtointeger(parameter(39))
     166         if (stringtointeger(parameter(39)) .NE. 1) then
    185167            print(" ")
    186168            print("Please set 'vector' to 0 or 1")
     
    192174   if ( .not. isvar("ref_mag") ) then                           ; sort of plots
    193175      ref_mag = 0.05
    194       if (parameter(71) .NE. "0.05") then
    195          ref_mag = stringtofloat(parameter(71))   
     176      if (parameter(47) .NE. "0.05") then
     177         ref_mag = stringtofloat(parameter(47))   
    196178      end if
    197179   end if
     
    208190   print(" ")
    209191   dim   = dimsizes(vNam)
     192
     193   ; ***************************************************
     194   ; check for kind of input file
     195   ; ***************************************************
     196
     197   check_3d = True
     198   do varn=0,dim-1
     199      checkxy = isStrSubset( vNam(varn),"xy")
     200      checkxz = isStrSubset( vNam(varn),"xz")
     201      checkyz = isStrSubset( vNam(varn),"yz")
     202      if (checkxy .OR. checkxz .OR. checkyz) then
     203         check_3d = False
     204         break
     205      end if
     206   end do
     207   if (.not. check_3d) then
     208      if (xyc .EQ. 1 .AND. .not. checkxy) then
     209         print(" ")
     210         print("Your input file doesn't have values for xy cross-sections")
     211         if (checkxz)then
     212            print("Select another input file or xz cross-sections")
     213         else
     214            print("Select another input file or yz cross-sections")
     215         end if
     216         print(" ")
     217         exit
     218      else
     219         print(" ")
     220         print("Your input file contains xy data")
     221         print(" ")
     222      end if
     223      if (xzc .EQ. 1 .AND. .not. checkxz) then
     224         print(" ")
     225         print("Your input file doesn't have values for xz cross-sections")
     226         if (checkxy)then
     227            print("Select another input file or xy cross-sections")
     228         else
     229            print("Select another input file or yz cross-sections")
     230         end if
     231         print(" ")
     232         exit
     233      else
     234         print(" ")
     235         print("Your input file contains xz data")
     236         print(" ")
     237      end if
     238      if (yzc .EQ. 1 .AND. .not. checkyz) then
     239         print(" ")
     240         print("Your input file doesn't have values for yz cross-sections")
     241         if (checkxy)then
     242            print("Select another input file or xy cross-sections")
     243         else
     244            print("Select another input file or xz cross-sections")
     245         end if
     246         print(" ")
     247         exit
     248      else
     249         print(" ")
     250         print("Your input file contains yz data")
     251         print(" ")
     252      end if
     253   else
     254      print(" ")
     255      print("Your input file: '"+file_in+"'")
     256      print("contains 3d or other data")
     257      print(" ")
     258   end if
     259   
    210260   if (dim .EQ. 0) then
    211261      print(" ")
    212       print("There are no data on file")
     262      print("There is no data on file")
    213263      print(" ")
     264      exit
    214265   end if
    215266
     
    241292   cs_res@lgLabelFontHeightF     = .02
    242293   cs_res@txFontHeightF          = .02
     294   cs_res@cnLevelSelectionMode    = "ManualLevels"
     295
    243296   cs_resP = True
    244297   cs_resP@txString               = f@title
    245298   cs_resP@txFuncCode             = "~"                 ; necessary for title
    246299   cs_resP@txFontHeightF          = .02
     300   
    247301 
    248302   if ( mode .eq. "Fill" ) then
     
    268322   ; *********************************************
    269323
    270    t    = f->time
    271    nt = dimsizes(t)
    272      
    273    ; *********************************************     
     324   t_all = f->time
     325   nt = dimsizes(t_all)
     326   delta_t = t_all(nt-1)/nt
     327
     328   ; ****************************************************       
    274329   ; start of time step and different types of mistakes that could be done
    275    ; *********************************************
     330   ; ****************************************************
    276331
    277332   if ( .not. isvar("start_time_step") ) then           
    278       start_time_step = 0
    279       if (parameter(13) .NE. "1") then
    280          if (parameter(13) .LT. "1")
    281             print(" ")
    282             print("Begin with time step 1")
     333      start_time_step=t_all(0)/3600
     334      if (parameter(13) .NE. "t(0)") then
     335         if (stringtodouble(parameter(13)) .GT. t_all(nt-1)/3600)
     336            print(" ")
     337            print("'start_time_step' = "+ parameter(13) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
     338            print(" ")
     339            print("Please select another 'start_time_step'")
    283340            print(" ")
    284341            exit
    285342         end if
    286          if (stringtointeger(parameter(13)) .GT. nt)
    287             print(" ")
    288             print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt))
     343         if (stringtofloat(parameter(13)) .LT. t_all(0)/3600)
     344            print(" ")
     345            print("'start_time_step' = "+ parameter(13) +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
     346            print(" ")
     347            print("Please select another 'start_time_step'")
    289348            print(" ")
    290349            exit
    291350         end if
    292          start_time_step = stringtointeger(parameter(13))-1
     351         start_time_step=stringtodouble(parameter(13))
    293352      end if
    294353   else
    295       if (start_time_step .LE. 0)
    296          print(" ")
    297          print("Begin with time step 1")
     354      if (start_time_step .GT. t_all(nt-1)/3600)
     355         print(" ")
     356         print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
     357         print(" ")
     358         print("Please select another 'start_time_step'")
    298359         print(" ")
    299360         exit
    300361      end if
    301       if (start_time_step .GT. nt)
    302          print(" ")
    303          print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt))
     362      if (start_time_step .LT. t_all(0)/3600)
     363         print(" ")
     364         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
     365         print(" ")
     366         print("Please select another 'start_time_step'")
    304367         print(" ")
    305368         exit
    306369      end if
    307       start_time_step = start_time_step - 1
    308    end if
    309 
     370   end if
     371   start_time_step = start_time_step*3600
     372   do i=0,nt-1     
     373      if (start_time_step .GE. t_all(i)-delta_t/2 .AND. start_time_step .LT. t_all(i)+delta_t/2)then
     374         st=i
     375         break
     376      end if
     377   end do
     378       
    310379   ; ****************************************************
    311380   ; end of time step and different types of mistakes that could be done
     
    313382
    314383   if ( .not. isvar("end_time_step") ) then             
    315       end_time_step = nt-1
    316       if (parameter(15) .NE. "nt") then
    317          if (parameter(15) .LE. "0")
    318             print(" ")
    319             print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
     384      end_time_step = t_all(nt-1)/3600
     385      if (parameter(15) .NE. "t(end)") then
     386         if (stringtodouble(parameter(15)) .LT. t_all(0)/3600)
     387            print(" ")
     388            print("'end_time_step' = "+parameter(15)+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
     389            print(" ")
     390            print("Please select another 'end_time_step'")
    320391            print(" ")
    321392            exit
    322393         end if
    323          if (stringtointeger(parameter(15)) .GT. nt)
    324             print(" ")
    325             print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt))
     394         if (stringtodouble(parameter(15)) .GT. t_all(nt-1)/3600)
     395            print(" ")
     396            print("'end_time_step' = "+ parameter(15) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
     397            print(" ")
     398            print("Please select another 'end_time_step'") 
    326399            print(" ")
    327400            exit
    328401         end if
    329          if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
    330             print(" ")
    331             print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(15))
     402         if (stringtodouble(parameter(15)) .LT. start_time_step/3600)
     403            print(" ")
     404            print("'end_time_step' = "+ parameter(15) +"h is lower than 'start_time_step' = "+parameter(13)+"h")
     405            print(" ")
     406            print("Please select another 'start_time_step' or 'end_time_step'")
    332407            print(" ")
    333408            exit
    334409         end if
    335          end_time_step = stringtointeger(parameter(15))-1 
     410         end_time_step = stringtodouble(parameter(15))
    336411      end if   
    337412   else
    338       if (end_time_step .LE. 0)
    339          print(" ")
    340          print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
     413      if (end_time_step .LT. t_all(0)/3600)
     414         print(" ")
     415         print("'end_time_step' = "+end_time_step+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
     416         print(" ")
     417         print("Please select another 'end_time_step'")
    341418         print(" ")
    342419         exit
    343420      end if
    344       if (end_time_step .GT. nt)
    345          print(" ")
    346          print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt))
    347          print(" ")
     421      if (end_time_step .GT. t_all(nt-1)/3600)
     422         print(" ")
     423         print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
     424         print(" ")
     425         print("Please select another 'end_time_step'") 
     426         print(" ")
    348427         exit
    349428      end if
    350       if (end_time_step .LT. start_time_step)
    351          print(" ")
    352          print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
     429      if (end_time_step .LT. start_time_step/3600)
     430         print(" ")
     431         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
     432         print(" ")
     433         print("Please select another 'start_time_step' or 'end_time_step'")
    353434         print(" ")
    354435         exit
    355436      end if
    356       end_time_step = end_time_step-1
    357    end if
    358 
     437   end if
     438   end_time_step = end_time_step*3600
     439   do i=0,nt-1     
     440      if (end_time_step .GE. t_all(i)-delta_t/2 .AND. end_time_step .LT. t_all(i)+delta_t/2)then
     441         et=i
     442         break
     443      end if
     444   end do
     445 
     446   delete(start_time_step)
     447   start_time_step=round(st,3)
     448   delete(end_time_step)
     449   end_time_step=round(et,3)
     450
     451   print(" ")
     452   print("Output from "+t_all(start_time_step)/3600+"h = "+t_all(start_time_step)+"s => index = "+start_time_step)
     453   print("       till "+t_all(end_time_step)/3600+"h = "+t_all(end_time_step)+"s => index = "+end_time_step)
     454   print(" ")
     455 
    359456   no_time=(end_time_step-start_time_step)+1
    360457
     
    365462   if (vector .EQ. 1) then
    366463      if ( .not. isvar("plotvec") ) then
    367          plotvec = parameter(69)
     464         plotvec = parameter(45)
    368465      end if
    369466      if ( .not. isvar("vec1") ) then
    370          vec1 = parameter(65)
    371          if (parameter(65) .EQ. "vec1") then
     467         vec1 = parameter(41)
     468         if (parameter(41) .EQ. "vec1") then
    372469            print(" ")
    373470            print("Please indicate Vector 1 ('vec1') for Vector-Plot")
     
    377474      end if
    378475      if ( .not. isvar("vec2") ) then
    379          vec2 = parameter(67)
    380          if (parameter(67) .EQ. "vec2") then
     476         vec2 = parameter(43)
     477         if (parameter(43) .EQ. "vec2") then
    381478            print(" ")
    382479            print("Please indicate Vector 2 ('vec2') for Vector-Plot")
     
    395492   ; ****************************************************
    396493
    397    do varn=dim-1,0,1 
    398 
    399       if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x" .or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then
    400          check = False
     494   if (xyc .EQ. 1) then
     495      do varn=0,dim-1
     496         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x")then   
     497            x_d     = f->$vNam(varn)$
     498            xdim    = dimsizes(x_d)-1
     499            delta_x = x_d(1)-x_d(0)
     500            break
     501         end if
     502      end do
     503      do varn=0,dim-1
     504         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then   
     505            y_d     = f->$vNam(varn)$
     506            ydim    = dimsizes(y_d)-1
     507            delta_y = y_d(1)-y_d(0)
     508            break
     509         end if
     510      end do
     511      do varn=0,dim-1
     512         if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
     513            z_d     = f->$vNam(varn)$
     514            zdim    = dimsizes(z_d)-1
     515            delta_z = 0
     516            break
     517         else
     518            if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
     519               z_d     = f->$vNam(varn)$
     520               zdim    = dimsizes(z_d)-1
     521               delta_z = -1.d
     522               break
     523            else
     524               zdim    = 0
     525               delta_z = -1.d
     526            end if
     527         end if
     528      end do
     529   end if
     530   if (xzc .EQ. 1) then
     531      do varn=0,dim-1
     532         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then
     533            x_d     = f->$vNam(varn)$
     534            xdim    = dimsizes(x_d)-1
     535            delta_x = x_d(1)-x_d(0)
     536            break
     537         end if
     538      end do
     539      do varn=0,dim-1   
     540         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
     541            z_d     = f->$vNam(varn)$
     542            zdim    = dimsizes(z_d)-1
     543            delta_z = z_d(1)-z_d(0)
     544            break
     545         end if
     546      end do
     547      do varn=0,dim-1
     548         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
     549            y_d     = f->$vNam(varn)$
     550            ydim    = dimsizes(y_d)-1
     551            delta_y = y_d(1)-y_d(0)
     552            break
     553         else
     554            if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
     555               y_d     = f->$vNam(varn)$
     556               ydim    = dimsizes(y_d)-1
     557               delta_y = -1.d
     558               break
     559            else
     560               ydim    = 0
     561               delta_y = -1.d
     562            end if
     563         end if
     564      end do
     565   end if
     566   if (yzc .EQ. 1) then
     567      do varn=0,dim-1 
     568         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
     569            y_d     = f->$vNam(varn)$
     570            ydim    = dimsizes(y_d)-1
     571            delta_y = y_d(1)-y_d(0)
     572            break
     573         end if
     574      end do
     575      do varn=0,dim-1
     576         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
     577            z_d     = f->$vNam(varn)$
     578            zdim    = dimsizes(z_d)-1
     579            delta_z = z_d(1)-z_d(0)
     580            break
     581         end if
     582      end do
     583      do varn=0,dim-1
     584         if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
     585            x_d     = f->$vNam(varn)$
     586            xdim    = dimsizes(x_d)-1
     587            delta_x = x_d(1)-x_d(0)
     588            break
     589         else
     590            if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
     591               x_d     = f->$vNam(varn)$
     592               xdim    = dimsizes(x_d)-1
     593               delta_x = -1.d
     594               break
     595            else
     596               xdim    = 0
     597               delta_x = -1.d
     598            end if
     599         end if
     600      end do
     601   end if
     602   
     603   ; ****************************************************
     604   ; set up ranges of x-, y- and z-coordinates
     605   ; ****************************************************         
     606                   
     607   if ( .not. isvar("xs") ) then               
     608      xs     = 0.0d
     609      xstart = 0
     610      if (parameter(63) .NE. "x0") then
     611         if (delta_x .EQ. -1) then
     612            print(" ")
     613            print("You cannot choose a start value for x, there are preseted layers for x")
     614            print(" ")
     615            xstart=0
     616         else
     617            if (stringtodouble(parameter(63)) .LT. 0-delta_x/2) then
     618               print(" ")
     619               print("range start for x-coordinate = "+parameter(63)+"m is lower than first value x = "+0+"m or xu = "+(0-delta_x/2)+"m")
     620               print(" ")
     621               exit
     622            end if
     623            if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
     624               if (stringtodouble(parameter(63)) .GE. xdim*delta_x) then
     625                  print(" ")
     626                  print("range start for x-coordinate = "+parameter(63)+"m is equal or greater than last value = "+xdim*delta_x+"m")
     627                  print(" ")
     628                  exit
     629               end if
     630            else
     631               if (stringtodouble(parameter(63)) .GT. xdim*delta_x) then
     632                  print(" ")
     633                  print("range start for x-coordinate = "+parameter(63)+"m is greater than last value = "+xdim*delta_x+"m")
     634                  print(" ")
     635                  exit
     636               end if
     637            end if
     638            xs = stringtodouble(parameter(63))
     639         end if
     640      end if
     641   else
     642      if (delta_x .EQ. -1) then
     643         print(" ")
     644         print("You cannot choose a start value for x, there are preseted layers for x")
     645         print(" ")
     646         xstart=0
    401647      else
    402          check = True
    403       end if   
    404       if (  isvar("var") ) then 
    405          check = isStrSubset( var,","+vNam(varn)+",")
    406       end if
    407       if (parameter(21) .NE. "variables") then
    408          var = parameter(21)
    409          check = isStrSubset( var,","+vNam(varn)+"," )
    410       end if   
    411 
    412       if(check) then
    413          print(vNam(varn))
    414          data_all = f->$vNam(varn)$
    415 
    416          ; ****************************************************
    417          ; set up ranges of x-, y- and z-coordinates
    418          ; ****************************************************         
    419                    
    420          xdim = dimsizes(data_all(0,0,0,:)) - 1   
    421          ydim = dimsizes(data_all(0,0,:,0)) - 1       
    422          zdim = dimsizes(data_all(0,:,0,0)) - 1 
    423        
    424          if ( .not. isvar("xe")) then     ; output x-coordinate range end (in index)
    425             xe = xdim
    426             if (parameter(53) .NE. "xdim") then
    427                if (stringtointeger(parameter(53)) .GT. xdim) then
    428                   print(" ")
    429                   print("range end for x-coordinate = "+parameter(53)+" is higher than available dimensions = "+xdim)
    430                   print(" ")
     648         if (xs .LT. 0-delta_x/2) then
     649            print(" ")
     650            print("range start for x-coordinate = "+xs+"m is lower than first value = "+(0-delta_x/2)+"m")
     651            print(" ")
     652            exit
     653         end if
     654         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
     655            if (xs .GE. xdim*delta_x) then
     656               print(" ")
     657               print("range start for x-coordinate = "+xs+"m is equal or greater than last value = "+xdim*delta_x+"m")
     658               print(" ")
     659               exit
     660            end if
     661         else
     662            if (xs .GT. xdim*delta_x) then
     663               print(" ")
     664               print("range start for x-coordinate = "+xs+"m is greater than last value = "+xdim*delta_x+"m")
     665               print(" ")
     666               exit
     667            end if
     668         end if
     669      end if
     670   end if
     671
     672   if ( .not. isvar("ys") ) then   
     673      ys = 0.0d
     674      ystart = 0
     675      if (parameter(67) .NE. "y0") then
     676         if (delta_y .EQ. -1) then
     677            print(" ")
     678            print("You cannot choose a start value for y, there are preseted layers for y")
     679            print(" ")
     680            ystart=0
     681         else
     682            if (stringtodouble(parameter(67)) .LT. 0-delta_y/2) then
     683               print(" ")
     684               print("range start for y-coordinate = "+parameter(67)+"m is lower than first value = "+0-delta_y/2+"m")
     685               print(" ")
     686               exit
     687            end if
     688            if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
     689               if (stringtodouble(parameter(67)) .GE. ydim*delta_y) then
     690                  print(" ")
     691                  print("range start for y-coordinate = "+parameter(67)+"m is equal or greater than last value = "+ydim*delta_y+"m")
     692                  print(" ")
    431693                  exit
    432694               end if
    433                if (stringtointeger(parameter(53)) .LT. 0 .OR. stringtointeger(parameter(53)) .LT. xs) then
    434                   print(" ")
    435                   print("range end for x-coordinate = "+parameter(53)+" is too small")
     695            else
     696               if (stringtodouble(parameter(67)) .GT. ydim*delta_y) then
     697                  print(" ")
     698                  print("range start for y-coordinate = "+parameter(67)+"m is greater than last value = "+ydim*delta_y+"m")
    436699                  print(" ")
    437700                  exit
    438701               end if
    439                xe = stringtointeger(parameter(53))
    440             end if
    441          else
    442             if (xe .GT. xdim) then
    443                print(" ")
    444                print("range end for x-coordinate = "+xe+" is higher than available dimensions = "+xdim)
    445                print(" ")
    446                exit
    447             end if
    448             if (xe .LT. 0 .OR. xe .LT. xs) then
    449                print(" ")
    450                print("range end for x-coordinate = "+xe+" is too small")
    451                print(" ")
    452                exit
    453             end if                 
    454          end if
    455 
    456          if ( .not. isvar("ye")) then     ; output y-coordinate range end (in index)
    457             ye = ydim
    458             if (parameter(57) .NE. "ydim") then
    459                if (stringtointeger(parameter(57)) .GT. ydim) then
    460                   print(" ")
    461                   print("range end for y-coordinate = "+parameter(57)+" is higher than available dimensions = "+ydim)
     702            end if
     703            ys = stringtodouble(parameter(67))
     704         end if
     705      end if
     706   else
     707      if (delta_y .EQ. -1) then
     708         print(" ")
     709         print("You cannot choose a start value for y, there are preseted layers for y")
     710         print(" ")
     711         ystart=0
     712      else
     713         if (ys .LT. 0-delta_y/2) then
     714            print(" ")
     715            print("range start for y-coordinate = "+ys+"m is lower than first value = "+0-delta_y/2+"m")
     716            print(" ")
     717            exit
     718         end if
     719         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
     720            if (ys .GE. ydim*delta_y) then
     721               print(" ")
     722               print("range start for y-coordinate = "+ys+"m is equal or greater than last value = "+ydim*delta_y+"m")
     723               print(" ")
     724               exit
     725            end if
     726         else
     727            if (ys .GT. ydim*delta_y) then
     728               print(" ")
     729               print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m")
     730               print(" ")
     731               exit
     732            end if
     733         end if
     734      end if
     735   end if
     736 
     737   if ( .not. isvar("zs") ) then                        ; start of z-coordinate range
     738      zs = 0
     739      if (parameter(71) .NE. "z0") then
     740         if (delta_z .EQ. -1) then
     741            print(" ")
     742            print("You cannot choose a start value for z, there are preseted layers for z")
     743            print(" ")
     744         else
     745            print(" ")
     746            print("Please mind to indicate start and end ranges for the z-coordinate in")
     747            print("indices not in 'meters'. Corresponding index and meter:")
     748            print(" ")
     749            print(" = "+z_d+" m")
     750            print(" ")
     751            if (stringtointeger(parameter(71)) .LT. 0) then
     752               print(" ")
     753               print("range start for z-coordinate = "+parameter(71)+" is lower than first gridpoint = 0")
     754               print(" ")
     755               exit
     756            end if
     757            if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
     758               if (stringtointeger(parameter(71)) .GE. zdim) then
     759                  print(" ")
     760                  print("range start for z-coordinate = "+parameter(71)+" is equal or greater than last gridpoint = "+zdim)
    462761                  print(" ")
    463762                  exit
    464763               end if
    465                if (stringtointeger(parameter(57)) .LT. 0 .OR. stringtointeger(parameter(57)) .LT. ys) then
    466                   print(" ")
    467                   print("range end for y-coordinate = "+parameter(57)+" is too small")
     764            else
     765               if (stringtodouble(parameter(71)) .GT. zdim) then
     766                  print(" ")
     767                  print("range start for z-coordinate = "+parameter(71)+" is greater than last gridpoint = "+zdim)
    468768                  print(" ")
    469769                  exit
    470770               end if
    471                ye = stringtointeger(parameter(57))
     771            end if
     772            zs = stringtointeger(parameter(71))
     773         end if
     774      end if
     775   else
     776      if (delta_z .EQ. -1) then
     777         print(" ")
     778         print("You cannot choose a start value for z, there are preseted layers for z")
     779         print(" ")
     780      else
     781         if (zs .LT. 0) then
     782            print(" ")
     783            print("range start for z-coordinate = "+zs+" is lower than first gridpoint = 0")
     784            print(" ")
     785            exit
     786         end if
     787         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
     788            if (zs .GE. zdim) then
     789               print(" ")
     790               print("range start for z-coordinate = "+zs+" is equal or greater than last gridpoint = "+zdim)
     791               print(" ")
     792               exit
     793            end if
     794         else
     795            if (zs .GT. zdim) then
     796               print(" ")
     797               print("range start for z-coordinate = "+zs+" is greater than last gridpoint = "+zdim)
     798               print(" ")
     799               exit
     800            end if
     801         end if
     802      end if
     803   end if 
     804
     805   if ( .not. isvar("xe")) then   
     806      xe   = xdim*delta_x
     807      xend = xdim
     808      if (parameter(65) .NE. "xdim") then
     809         if (delta_x .EQ. -1) then
     810            print(" ")
     811            print("You cannot choose an end value for x, there are preseted layers for x")
     812            print(" ")
     813            xend=xdim           
     814         else
     815            if (stringtodouble(parameter(65)) .GT. xdim*delta_x) then
     816               print(" ")
     817               print("range end for x-coordinate = "+parameter(65)+"m is greater than last value = "+xdim*delta_x+"m")
     818               print(" ")
     819               exit
    472820            end if
    473          else
    474             if (ye .GT. ydim) then
    475                print(" ")
    476                print("range end for y-coordinate = "+ye+" is higher than available dimensions = "+ydim)
    477                print(" ")
    478                exit
    479             end if
    480             if (ye .LT. 0 .OR. ye .LT. ys) then
    481                print(" ")
    482                print("range end for y-coordinate = "+ye+" is too small")
    483                print(" ")
    484                exit
    485             end if
    486          end if
    487 
    488          if ( .not. isvar("ze")) then     ; output z-coordinate range end (in index)
    489             ze = zdim
    490             if (parameter(61) .NE. "zdim") then
    491                if (stringtointeger(parameter(61)) .GT. zdim) then
    492                   print(" ")
    493                   print("range end for z-coordinate = "+parameter(61)+" is higher than available dimensions = "+zdim)
     821            if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
     822               if (stringtodouble(parameter(65)) .LE. 0-delta_x/2)
     823                  print(" ")
     824                  print("range end for x-coordinate = "+parameter(65)+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
    494825                  print(" ")
    495826                  exit
    496827               end if
    497                if (stringtointeger(parameter(61)) .LT. 0 .OR. stringtointeger(parameter(61)) .LT. zs) then
    498                   print(" ")
    499                   print("range end for z-coordinate = "+parameter(61)+" is too small")
     828               if (stringtodouble(parameter(65)) .LE. xs) then
     829                  print(" ")
     830                  print("range end for x-coordinate = "+parameter(65)+"m is equal or lower than start range = "+xs+"m")
    500831                  print(" ")
    501832                  exit
    502833               end if
    503                ze = stringtointeger(parameter(61))
     834            else
     835               if (stringtodouble(parameter(65)) .LT. 0-delta_x/2)
     836                  print(" ")
     837                  print("range end for x-coordinate = "+parameter(65)+"m is lower than first value = "+(0-delta_x/2)+"m")
     838                  print(" ")
     839                  exit
     840               end if
     841               if (stringtodouble(parameter(65)) .LT. xs) then
     842                  print(" ")
     843                  print("range end for x-coordinate = "+parameter(65)+"m is lower than start range = "+xs+"m")
     844                  print(" ")
     845                  exit
     846               end if
     847            end if         
     848            xe = stringtodouble(parameter(65))
     849         end if
     850      end if
     851   else
     852      if (delta_x .EQ. -1) then
     853         print(" ")
     854         print("You cannot choose an end value for x, there are preseted layers for x")
     855         print(" ")
     856         xend=xdim
     857      else
     858         if (xe .GT. xdim*delta_x) then
     859            print(" ")
     860            print("range end for x-coordinate = "+xe+"m is greater than last value = "+xdim*delta_x+"m")
     861            print(" ")
     862            exit
     863         end if
     864         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
     865            if (xe .LE. 0-delta_x/2)
     866               print(" ")
     867               print("range end for x-coordinate = "+xe+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
     868               print(" ")
     869               exit
     870            end if
     871            if (xe .LE. xs) then
     872               print(" ")
     873               print("range end for x-coordinate = "+xe+"m is equal or lower than start range = "+xs+"m")
     874               print(" ")
     875               exit
     876            end if
     877            if (stringtodouble(xe .EQ. xs+1)) then
     878               print(" ")
     879               print("range end for x-coordinate = "+xe+"m must be at least two more gridpoints greater than start range = "+xs+"m")
     880               print(" ")
     881               exit
     882            end if
     883         else
     884            if (xe .LT. 0-delta_x/2)
     885               print(" ")
     886               print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m")
     887               print(" ")
     888               exit
     889            end if
     890            if (xe .LT. xs) then
     891               print(" ")
     892               print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m")
     893               print(" ")
     894               exit
     895            end if
     896         end if
     897      end if               
     898   end if
     899
     900   if ( .not. isvar("ye")) then 
     901      ye   = ydim*delta_y
     902      yend = ydim
     903      if (parameter(69) .NE. "ydim") then
     904         if (delta_y .EQ. -1) then
     905            print(" ")
     906            print("You cannot choose an end value for y, there are preseted layers for y")
     907            print(" ")
     908            yend=ydim
     909         else
     910            if (stringtodouble(parameter(69)) .GT. ydim*delta_y) then
     911               print(" ")
     912               print("range end for y-coordinate = "+parameter(69)+"m is greater than last value = "+ydim*delta_y+"m")
     913               print(" ")
     914               exit
    504915            end if
    505          else
    506             if (ze .GT. zdim) then
    507                print(" ")
    508                print("range end for z-coordinate = "+ze+" is higher than available dimensions = "+zdim)
     916            if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
     917               if (stringtodouble(parameter(69)) .LE. 0-delta_y/2)
     918                  print(" ")
     919                  print("range end for y-coordinate = "+parameter(69)+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
     920                  print(" ")
     921                  exit
     922               end if
     923               if (stringtodouble(parameter(69)) .LE. ys) then
     924                  print(" ")
     925                  print("range end for y-coordinate = "+parameter(69)+"m is equal or lower than start range = "+ys+"m")
     926                  print(" ")
     927                  exit
     928               end if
     929               if (stringtodouble(parameter(69)) .EQ. ys+1) then
     930                  print(" ")
     931                  print("range end for y-coordinate = "+parameter(69)+"m must be at least two more gridpoints greater than start range = "+ys+"m")
     932                  print(" ")
     933                  exit
     934               end if
     935            else
     936               if (stringtodouble(parameter(69)) .LT. 0-delta_y/2)
     937                  print(" ")
     938                  print("range end for y-coordinate = "+parameter(69)+"m is lower than first value = "+(0-delta_y/2)+"m")
     939                  print(" ")
     940                  exit
     941               end if
     942               if (stringtodouble(parameter(69)) .LT. ys) then
     943                  print(" ")
     944                  print("range end for y-coordinate = "+parameter(69)+"m is lower than start range = "+ys+"m")
     945                  print(" ")
     946                  exit
     947               end if
     948            end if         
     949            ye = stringtodouble(parameter(69))
     950         end if
     951      end if
     952   else
     953      if (delta_y .EQ. -1) then
     954         print(" ")
     955         print("You cannot choose an end value for y, there are preseted layers for y")
     956         print(" ")
     957         yend=ydim
     958      else
     959         if (ye .GT. ydim*delta_y) then
     960            print(" ")
     961            print("range end for y-coordinate = "+ye+"m is greater than last value = "+ydim*delta_y+"m")
     962            print(" ")
     963            exit
     964         end if
     965         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
     966            if (ye .LE. 0-delta_y/2)
     967               print(" ")
     968               print("range end for y-coordinate = "+ye+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
     969               print(" ")
     970               exit
     971            end if
     972            if (ye .LE. ys) then
     973               print(" ")
     974               print("range end for y-coordinate = "+ye+"m is equal or lower than start range = "+ys+"m")
     975               print(" ")
     976               exit
     977            end if
     978            if (ye .EQ. ys+1) then
     979               print(" ")
     980               print("range end for y-coordinate = "+ye+"m must be at least two more gridpoints greater than start range = "+ys+"m")
     981               print(" ")
     982               exit
     983            end if
     984         else
     985            if (ye .LT. 0-delta_y/2)
     986               print(" ")
     987               print("range end for y-coordinate = "+ye+"m is lower than first value = "+(0-delta_y/2)+"m")
     988               print(" ")
     989               exit
     990            end if
     991            if (ye .LT. ys) then
     992               print(" ")
     993               print("range end for y-coordinate = "+ye+"m is lower than start range = "+ys+"m")
     994               print(" ")
     995               exit
     996            end if
     997         end if
     998      end if
     999   end if
     1000 
     1001   if ( .not. isvar("ze")) then 
     1002      ze = zdim
     1003      if (parameter(73) .NE. "zdim") then
     1004         if (delta_z .EQ. -1) then
     1005            print(" ")
     1006            print("You cannot choose an end value for z, there are preseted layers for z")
     1007            print(" ")           
     1008         else
     1009            if (stringtointeger(parameter(73)) .GT. zdim) then
     1010               print(" ")
     1011               print("range end for z-coordinate = "+parameter(73)+" is greater than last gridpoint = "+zdim)
    5091012               print(" ")
    510                exit
    511             end if
    512             if (ze .LT. 0 .OR. ze .LT. zs) then
    513                print(" ")
    514                print("range end for z-coordinate = "+ze+" is too small")
    515                print(" ")
    516                exit
    517             end if
    518          end if   
    519          delete(data_all)
    520       end if 
    521    end do
    522 
     1013               exit
     1014            end if
     1015         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
     1016            if (stringtointeger(parameter(73)) .LE. 0)
     1017               print(" ")
     1018               print("range end for z-coordinate = "+parameter(73)+" is equal or lower than first gridpoint = 0")
     1019               print(" ")
     1020               exit
     1021            end if
     1022            if (stringtointeger(parameter(73)) .LE. zs) then
     1023               print(" ")
     1024               print("range end for z-coordinate = "+parameter(73)+" is equal or lower than start range = "+zs)
     1025               print(" ")
     1026               exit
     1027            end if
     1028            if (stringtodouble(parameter(73)) .EQ. zs+1) then
     1029               print(" ")
     1030               print("range end for z-coordinate = "+parameter(73)+" must be at least two more gridpoints greater than start range = "+zs)
     1031               print(" ")
     1032               exit
     1033            end if
     1034         else
     1035            if (stringtointeger(parameter(73)) .LT. 0)
     1036               print(" ")
     1037               print("range end for z-coordinate = "+parameter(73)+" is lower than first gridpoint = 0")
     1038               print(" ")
     1039               exit
     1040            end if
     1041            if (stringtointeger(parameter(73)) .LT. zs) then
     1042               print(" ")
     1043               print("range end for z-coordinate = "+parameter(73)+" is lower than start range = "+zs)
     1044               print(" ")
     1045               exit
     1046            end if
     1047         end if         
     1048         ze = stringtointeger(parameter(73))
     1049      end if
     1050      end if
     1051   else
     1052      if (delta_z .EQ. -1) then
     1053         print(" ")
     1054         print("You cannot choose an end value for z, there are preseted layers for z")
     1055         print(" ")
     1056         ze = zdim
     1057      else
     1058         if (ze .GT. zdim) then
     1059            print(" ")
     1060            print("range end for z-coordinate = "+ze+" is greater than last gridpoint = "+zdim)
     1061            print(" ")
     1062            exit
     1063         end if
     1064         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
     1065            if (ze .LE. 0)
     1066               print(" ")
     1067               print("range end for z-coordinate = "+ze+" is equal or lower than first gridpoint = 0")
     1068               print(" ")
     1069               exit
     1070            end if
     1071            if (ze .LE. zs) then
     1072               print(" ")
     1073               print("range end for z-coordinate = "+ze+" is equal or lower than start range = "+zs)
     1074               print(" ")
     1075               exit
     1076            end if
     1077            if (ze .EQ. zs+1) then
     1078               print(" ")
     1079               print("range end for z-coordinate = "+ze+" must be at least two more gridpoints greater than start range = "+zs)
     1080               print(" ")
     1081               exit
     1082            end if
     1083         else
     1084            if (ze .LT. 0)
     1085               print(" ")
     1086               print("range end for z-coordinate = "+ze+" is lower than first gridpoint = 0")
     1087               print(" ")
     1088               exit
     1089            end if
     1090            if (ze .LT. zs) then
     1091               print(" ")
     1092               print("range end for z-coordinate = "+ze+" is lower than start range = "+zs)
     1093               print(" ")
     1094               exit
     1095            end if
     1096         end if
     1097      end if
     1098   end if
     1099
     1100   if (delta_x .NE. -1) then
     1101      do i=0,xdim   
     1102         if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
     1103            xstart=i
     1104            break
     1105         end if
     1106      end do
     1107      do i=0,xdim   
     1108         if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
     1109            xend=i
     1110            break
     1111         end if
     1112      end do
     1113   end if
     1114   if (delta_y .NE. -1) then
     1115      do i=0,xdim   
     1116         if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
     1117            ystart=i
     1118            break
     1119         end if
     1120      end do
     1121      do i=0,xdim   
     1122         if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
     1123            yend=i
     1124            break
     1125         end if
     1126      end do
     1127   end if
     1128   
     1129   delete(xs)
     1130   delete(xe)   
     1131   delete(ys)
     1132   delete(ye)
     1133
     1134   xs=xstart
     1135   xe=xend
     1136   ys=ystart
     1137   ye=yend
     1138
     1139   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
     1140      if (xe .EQ. xs+1) then
     1141         print(" ")
     1142         print("range end for x-coordinate="+parameter(65)+"m(=>value on file="+xe*delta_x+"m) must be at least two")
     1143         print("more gridpoints(="+2*delta_x+"m) greater than start range="+parameter(63)+"m(=>value on file="+xs*delta_x+"m)")
     1144         print(" ")
     1145         exit
     1146      end if
     1147   end if
     1148   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
     1149      if (ye .EQ. ys+1) then
     1150         print(" ")
     1151         print("range end for y-coordinate="+parameter(69)+"m(=>value on file="+ye*delta_y+"m) must be at least two")
     1152         print("more gridpoints(="+2*delta_y+"m) greater than start range="+parameter(67)+"m(=>value on file="+ys*delta_y+"m)")
     1153         print(" ")
     1154         exit
     1155      end if
     1156   end if
     1157   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
     1158      if (ze .EQ. zs+1) then
     1159         print(" ")
     1160         print("range end for x-coordinate="+parameter(73)+"(=> value on file="+ze+") must be at least two")
     1161         print("more gridpoints greater than start range="+parameter(71)+"(=>value on file="+zs+")")
     1162         print(" ")
     1163         exit
     1164      end if
     1165   end if
     1166 
    5231167   if (xyc .EQ. 1) then
    5241168      no_layer = (ze-zs)+1
     1169      data = new((/dim,nt,zdim+1,(ye-ys)+1,(xe-xs)+1/),float)
    5251170   end if
    5261171   if (xzc .EQ. 1) then
    5271172      no_layer = (ye-ys)+1
     1173      data = new((/dim,nt,(ze-zs)+1,ydim+1,(xe-xs)+1/),float)
    5281174   end if
    5291175   if (yzc .EQ. 1) then
    5301176      no_layer = (xe-xs)+1
    531    end if
    532 
    533    data  = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
     1177      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,xdim+1/),float)
     1178   end if
     1179
     1180   MinVal = new(dim,float)
     1181   MaxVal = new(dim,float)
    5341182
    5351183   ; ****************************************************
     
    5711219     
    5721220      data_all = f->$vNam(varn)$
    573 
     1221   
    5741222      if (vector .EQ. 1) then   
    5751223         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
     
    5901238      end if 
    5911239
    592       if(check) then 
     1240      if(check) then
    5931241
    5941242         no_var=no_var+1
     
    6031251         end if
    6041252
    605          data(varn,:,:,:,:)=data_all(0:nt-1,zs:ze,ys:ye,xs:xe)
     1253         if (xyc .EQ. 1) then
     1254            data(varn,:,:,:,:)=data_all(:,:,ys:ye,xs:xe)
     1255         end if
     1256         if ( xzc .eq. 1 ) then
     1257            data(varn,:,:,:,:)=data_all(:,zs:ze,:,xs:xe)
     1258         end if
     1259         if ( yzc .eq. 1) then
     1260            data(varn,:,:,:,:)=data_all(:,zs:ze,ys:ye,:)
     1261         end if
     1262
    6061263         if (check_vec1) then
    6071264            vect1=data_all
     
    6151272         data!2 = "z"
    6161273         data!3 = "y"
    617          data!4 = "x"   
     1274         data!4 = "x" 
     1275
     1276         MinVal(varn) = min(data(varn,:,:,:,:))
     1277         MaxVal(varn) = max(data(varn,:,:,:,:))
    6181278
    6191279      end if
     
    6221282
    6231283   end do
    624  
    625    plot_panel=new((/no_time*no_layer*no_var/),graphic)
     1284
     1285   if (no_var .EQ. 0) then
     1286      print(" ")
     1287      print("Please select a variable 'var=' or use the default value")
     1288      print(" ")
     1289      print("Your selection '"+var+"' does not exist on the input file")
     1290      print(" ")
     1291      exit
     1292   end if
    6261293
    6271294   if (vector .EQ. 1) then
     
    6481315   gsn_define_colormap(wks_ps,"rainbow+white")
    6491316
    650    plot=new((/no_time*no_layer*no_var/),graphic)
    651 
     1317   plot=new((/no_time*no_layer*no_var*2/),graphic)
     1318 
    6521319   page = 0
    6531320   n=0
    654 
    655    print("plots sorted by '"+sort+"'")
     1321   print(" ")
     1322   print("Plots sorted by '"+sort+"'")
    6561323   print(" ")
    6571324
     
    6601327   ; ***************************************************
    6611328
     1329   if (vector .EQ. 1 .AND. parameter(45) .EQ. "plotvec") then
     1330      do lo = los, loe                                         
     1331         do li = lis, lie
     1332            level = "=" + data&z(lo) + "m"       
     1333            vecres                  = True            ; vector only resources
     1334            vecres@gsnDraw          = False           ; don't draw
     1335            vecres@gsnFrame         = False           ; don't advance frame
     1336            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1337            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1338            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
     1339            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
     1340            vecres@gsnLeftString = "t=" + data&t(li) +"s  z"+level
     1341            vecres@tiXAxisString    = " "                                               
     1342            plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
     1343            n=n+1
     1344         end do
     1345      end do
     1346   end if
     1347
    6621348   check = True
    663 
     1349 
    6641350   do varn=dim-1,0,1
    6651351
     
    6821368   
    6831369      if(check) then
    684    
     1370         print(vNam(varn))
     1371         space=(MaxVal(varn)-MinVal(varn))/24
     1372 
     1373         cs_res@cnMinLevelValF = MinVal(varn)
     1374         cs_res@cnMaxLevelValF = MaxVal(varn)
     1375
     1376         cs_res@cnLevelSpacingF = space
     1377     
    6851378         ; ****************************************************
    6861379         ; loops over time and layer
    6871380         ; ****************************************************
    6881381
     1382         
    6891383         do lo = los, loe                                       
    690             do li = lis, lie                           
    691                  
     1384            do li = lis, lie
     1385
    6921386               ; ****************************************************
    6931387               ; xy cross section
     
    6991393                  cs_res@tiYAxisString = "y [m]"
    7001394                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
    701      
     1395                 
    7021396                  if ( sort .eq. "time" ) then
    7031397                     if ( data&z(lo) .eq. -1 ) then
     
    7061400                        level = "=" + data&z(lo) + "m"
    7071401                     end if
    708                      cs_res@gsnCenterString = "time=" + data&t(li) +"s  z"+level
     1402                     cs_res@gsnCenterString = "t=" + data&t(li) +"s  z"+level
    7091403                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
    710                      if (vector .EQ. 1) then
    711                         if (check_vecp)
    712                            vecres                  = True            ; vector only resources
    713                            vecres@gsnDraw          = False           ; don't draw
    714                            vecres@gsnFrame         = False           ; don't advance frame
    715                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    716                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    717                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    718                            vecres@gsnRightString   = " "             ; turn off right string
    719                            vecres@gsnLeftString    = " "             ; turn off left string
    720                            vecres@tiXAxisString    = " "   
    721                            plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
    722                            overlay(plot(n), plot_vec)
    723                         else 
    724                            vecres                  = True            ; vector only resources
    725                            vecres@gsnDraw          = False           ; don't draw
    726                            vecres@gsnFrame         = False           ; don't advance frame
    727                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    728                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    729                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    730                            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
    731                            vecres@gsnLeftString    = " "             ; turn off left string
    732                            vecres@tiXAxisString    = " "   
    733                            plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
    734                            plot(n) = plot_vec
    735                         end if
     1404                     if (vector .EQ. 1 .AND. check_vecp) then
     1405                        vecres                  = True            ; vector only resources
     1406                        vecres@gsnDraw          = False           ; don't draw
     1407                        vecres@gsnFrame         = False           ; don't advance frame
     1408                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1409                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1410                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref 
     1411                        vecres@gsnRightString   = " "
     1412                        vecres@gsnLeftString    = " "
     1413                        vecres@tiXAxisString    = " "   
     1414                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
     1415                        overlay(plot(n), plot_vec)
    7361416                     end if                         
    7371417                  end if
    738 
     1418                 
    7391419                  if ( sort .eq. "layer" ) then
    7401420                     if ( data&z(li) .eq. -1 ) then
     
    7451425                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  z"+ level
    7461426                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
    747                      if (vector .EQ. 1) then
    748                         if (check_vecp)
    749                            vecres                  = True            ; vector only resources
    750                            vecres@gsnDraw          = False           ; don't draw
    751                            vecres@gsnFrame         = False           ; don't advance frame
    752                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    753                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    754                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    755                            vecres@gsnRightString   = " "             ; turn off right string
    756                            vecres@gsnLeftString    = " "             ; turn off left string
    757                            vecres@tiXAxisString    = " "   
    758                            plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
    759                            overlay(plot(n), plot_vec)
    760                         else 
    761                            vecres                  = True            ; vector only resources
    762                            vecres@gsnDraw          = False           ; don't draw
    763                            vecres@gsnFrame         = False           ; don't advance frame
    764                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    765                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    766                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    767                            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
    768                            vecres@gsnLeftString    = " "             ; turn off left string
    769                            vecres@tiXAxisString    = " "   
    770                            plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
    771                            plot(n) = plot_vec
    772                         end if
     1427                     if (vector .EQ. 1 .AND. check_vecp) then
     1428                        vecres                  = True            ; vector only resources
     1429                        vecres@gsnDraw          = False           ; don't draw
     1430                        vecres@gsnFrame         = False           ; don't advance frame
     1431                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1432                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1433                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
     1434                        vecres@gsnRightString   = " "             ; turn off right string
     1435                        vecres@gsnLeftString    = " "             ; turn off left string
     1436                        vecres@tiXAxisString    = " "   
     1437                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
     1438                        overlay(plot(n), plot_vec)
    7731439                     end if
    7741440                  end if
     
    7841450                  cs_res@tiYAxisString   = "z [m]"
    7851451                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
    786 
     1452               
    7871453                  if ( sort .eq. "time" ) then
    7881454                     if ( data&z(lo) .eq. -1 ) then
     
    7931459                     cs_res@gsnCenterString = "t=" + data&t(li) + "s  y"+ level
    7941460                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
    795                      if (vector .EQ. 1) then
    796                         if (check_vecp)
    797                            vecres                  = True            ; vector only resources
    798                            vecres@gsnDraw          = False           ; don't draw
    799                            vecres@gsnFrame         = False           ; don't advance frame
    800                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    801                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    802                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    803                            vecres@gsnRightString   = " "             ; turn off right string
    804                            vecres@gsnLeftString    = " "             ; turn off left string
    805                            vecres@tiXAxisString    = " "   
    806                            plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
    807                            overlay(plot(n), plot_vec)
    808                         else 
    809                            vecres                  = True            ; vector only resources
    810                            vecres@gsnDraw          = False           ; don't draw
    811                            vecres@gsnFrame         = False           ; don't advance frame
    812                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    813                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    814                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    815                            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
    816                            vecres@gsnLeftString    = " "             ; turn off left string
    817                            vecres@tiXAxisString    = " "   
    818                            plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
    819                            plot(n) = plot_vec
    820                         end if
     1461                     if (vector .EQ. 1 .AND. check_vecp) then
     1462                        vecres                  = True            ; vector only resources
     1463                        vecres@gsnDraw          = False           ; don't draw
     1464                        vecres@gsnFrame         = False           ; don't advance frame
     1465                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1466                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1467                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
     1468                        vecres@gsnRightString   = " "             ; turn off right string
     1469                        vecres@gsnLeftString    = " "             ; turn off left string
     1470                        vecres@tiXAxisString    = " "   
     1471                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
     1472                        overlay(plot(n), plot_vec)
    8211473                     end if
    8221474                  end if
     
    8301482                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  y"+ level
    8311483                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
    832                      if (vector .EQ. 1) then
    833                         if (check_vecp)
    834                            vecres                  = True            ; vector only resources
    835                            vecres@gsnDraw          = False           ; don't draw
    836                            vecres@gsnFrame         = False           ; don't advance frame
    837                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    838                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    839                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    840                            vecres@gsnRightString   = " "             ; turn off right string
    841                            vecres@gsnLeftString    = " "             ; turn off left string
    842                            vecres@tiXAxisString    = " "   
    843                            plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
    844                            overlay(plot(n), plot_vec)
    845                         else 
    846                            vecres                  = True            ; vector only resources
    847                            vecres@gsnDraw          = False           ; don't draw
    848                            vecres@gsnFrame         = False           ; don't advance frame
    849                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    850                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    851                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    852                            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
    853                            vecres@gsnLeftString    = " "             ; turn off left string
    854                            vecres@tiXAxisString    = " "   
    855                            plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
    856                            plot(n) = plot_vec
    857                         end if
     1484                     if (vector .EQ. 1 .AND. check_vecp) then
     1485                        vecres                  = True            ; vector only resources
     1486                        vecres@gsnDraw          = False           ; don't draw
     1487                        vecres@gsnFrame         = False           ; don't advance frame
     1488                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1489                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1490                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
     1491                        vecres@gsnRightString   = " "             ; turn off right string
     1492                        vecres@gsnLeftString    = " "             ; turn off left string
     1493                        vecres@tiXAxisString    = " "   
     1494                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
     1495                        overlay(plot(n), plot_vec)
    8581496                     end if
    8591497                  end if                 
     
    8781516                     cs_res@gsnCenterString = "t=" + data&t(li) + "s  x"+ level
    8791517                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
    880                      if (vector .EQ. 1) then
    881                         if (check_vecp)
    882                            vecres                  = True            ; vector only resources
    883                            vecres@gsnDraw          = False           ; don't draw
    884                            vecres@gsnFrame         = False           ; don't advance frame
    885                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    886                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    887                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    888                            vecres@gsnRightString   = " "             ; turn off right string
    889                            vecres@gsnLeftString    = " "             ; turn off left string
    890                            vecres@tiXAxisString    = " "   
    891                            plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
    892                            overlay(plot(n), plot_vec)
    893                         else 
    894                            vecres                  = True            ; vector only resources
    895                            vecres@gsnDraw          = False           ; don't draw
    896                            vecres@gsnFrame         = False           ; don't advance frame
    897                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    898                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    899                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    900                            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
    901                            vecres@gsnLeftString    = " "             ; turn off left string
    902                            vecres@tiXAxisString    = " "   
    903                            plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
    904                            plot(n) = plot_vec
    905                         end if
     1518                     if (vector .EQ. 1 .AND. check_vecp) then
     1519                        vecres                  = True            ; vector only resources
     1520                        vecres@gsnDraw          = False           ; don't draw
     1521                        vecres@gsnFrame         = False           ; don't advance frame
     1522                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1523                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1524                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
     1525                        vecres@gsnRightString   = " "             ; turn off right string
     1526                        vecres@gsnLeftString    = " "             ; turn off left string
     1527                        vecres@tiXAxisString    = " "   
     1528                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
     1529                        overlay(plot(n), plot_vec)
    9061530                     end if
    9071531                  end if
     
    9151539                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  x"+ level
    9161540                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
    917                      if (vector .EQ. 1) then
    918                         if (check_vecp)
    919                            vecres                  = True            ; vector only resources
    920                            vecres@gsnDraw          = False           ; don't draw
    921                            vecres@gsnFrame         = False           ; don't advance frame
    922                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    923                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    924                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    925                            vecres@gsnRightString   = " "             ; turn off right string
    926                            vecres@gsnLeftString    = " "             ; turn off left string
    927                            vecres@tiXAxisString    = " "   
    928                            plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
    929                            overlay(plot(n), plot_vec)
    930                         else 
    931                            vecres                  = True            ; vector only resources
    932                            vecres@gsnDraw          = False           ; don't draw
    933                            vecres@gsnFrame         = False           ; don't advance frame
    934                            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
    935                            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
    936                            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
    937                            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
    938                            vecres@gsnLeftString    = " "             ; turn off left string
    939                            vecres@tiXAxisString    = " "   
    940                            plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
    941                            plot(n) = plot_vec
    942                         end if
     1541                     if (vector .EQ. 1 .AND. check_vecp)then
     1542                        vecres                  = True            ; vector only resources
     1543                        vecres@gsnDraw          = False           ; don't draw
     1544                        vecres@gsnFrame         = False           ; don't advance frame
     1545                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
     1546                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
     1547                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
     1548                        vecres@gsnRightString   = " "             ; turn off right string
     1549                        vecres@gsnLeftString    = " "             ; turn off left string
     1550                        vecres@tiXAxisString    = " "   
     1551                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
     1552                        overlay(plot(n), plot_vec)
    9431553                     end if
    9441554                  end if
     
    9471557            end do     
    9481558         end do 
    949       end if 
     1559      end if
    9501560   end do
    9511561
     
    9531563   ; merge plots onto one page
    9541564   ; ***************************************************
    955              
    956    if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
    957       gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
    958       print(" ")
    959       print("Outputs to .eps or .epsi have only one frame")
    960       print(" ")
    961    else
    962       do np = 0,no_layer*no_time*no_var-1,no_lines*no_columns   
    963          if ( np + no_lines*no_columns .gt. (no_layer*no_time*no_var)-1) then
    964             gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_lines,no_columns/),cs_resP)
    965          else
    966             gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
    967          end if
    968       end do
     1565
     1566   if (vector .EQ. 1 .AND. parameter(45) .EQ. "plotvec")then
     1567      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
     1568         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
     1569         print(" ")
     1570         print("Outputs to .eps or .epsi have only one frame")
     1571         print(" ")
     1572      else
     1573         do np = 0,no_layer*no_time*(no_var+1)-1,no_lines*no_columns   
     1574            if ( np + no_lines*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then
     1575               gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-1),(/no_lines,no_columns/),cs_resP)
     1576            else
     1577               gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
     1578            end if
     1579         end do
     1580      end if
     1581   else         
     1582      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
     1583         gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
     1584         print(" ")
     1585         print("Outputs to .eps or .epsi have only one frame")
     1586         print(" ")
     1587      else
     1588         do np = 0,no_layer*no_time*no_var-1,no_lines*no_columns   
     1589            if ( np + no_lines*no_columns .gt. (no_layer*no_time*no_var)-1) then
     1590               gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_lines,no_columns/),cs_resP)
     1591            else
     1592               gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
     1593            end if
     1594         end do
     1595      end if
    9691596   end if
    9701597
Note: See TracChangeset for help on using the changeset viewer.