text:
Sorrento Apartments is just across the street from Safeway, Rite Aid, and Office Max. It is the second apartment complex after Oakshade Commons Apartments Oakshade on the left of Valdora St., and is accessible via the W bus line W and M bus line M bus line. While it is not as posh as the surrounding apartments such as Avalon Apartments Avalon and Tanglewood, it offers a comfortable place for both students and families. There is a nice mixture of children in the day and drunk drunken students in the night during weekends. Current management is thought to be better than preOctober 2010.

They offer five floor plans in single and two story models. Features include a pool, spa and exercise room.

Sorrento Apartments also manages Brisa Villas and Avalon Apartments in South Davis.

 2012/2013 Rental Information

2 Bedroom Flats:           $1,375  $1,410
2 Bedroom Townhouses: $1,510  $1,550
3 Bedroom Flats:           $1,880  $1,910
3 Bedroom Townhouses: $1,995  $2,045
4 Bedroom Townhouses: $2,315  $2,490

Pictures

To learn more about rental housing in Davis, check out our Housing Guide as well as our Apartments pages.

  /Pre 2009  Reviews before 2009 Property Management Company change

20090123 21:28:33 nbsp I have been living in Sorrento for the past year and it has been fairly nice. The park, maintenance, size of the apartments, and general location are all ideal. However; the price that was initially advertised was not the price that I was given. Also, the management will constantly put letters in your mail if they have a complaint, which made me wonder why Im paying money to someone who cant even speak to me like a human being. Theres plenty of amenities to take advantage of, but the area constantly smells like kim chi, which is probably the most foul scent I have ever had the misfortune of smelling in my life. This smell comes from the 95% Asian population within Sorrento, so if you are looking to live in a diverse, social, attractive area then Sorrento is not for you. It certainly wasnt right for me. Users/Yiddo9

20090221 11:56:05 nbsp I have enjoyed my stay at Sorrento.  I moved in at the beginning of sept 08.  There were some problems when we moved in.  Rat droppings, plugs not working, windows not opening, window screens falling off, etc.  HOWEVER, maintenance was amazing and fixed everything within the week, if not the next day.  If we ever had a problem involving needing maintenance, they have came and fixed it by the end of the day.  Both of the maintenance guys are friendly and nice. We are yet to see a rat.  The light in the oven was out, and they simply just replace the stove completely.  I havent ever really talked to Linda, but when I see her around the complex she always smiles.  The staff is always extremely helpful.  When any of my roommates and I walk into the office, Bryan or Ryan Im still not sure :/, always knows which apartment I belong to and usually can figure out why I am in there right away.  The outgoing mail has been broken into a few times, which would have been nice if we were notified about it.  However, when asked about mail not being received they admitted and told us they cant actually do anything with the mail boxes because it is US postal service property.  The park is right behind the complex, great for just a walk or relaxing.  The floorplans are well thought out, despite what some previous comments said.  I rarely smell anything bad.  Landscaping does come once a week and starts somewhat early if you like to sleep in during the week.  The landscapers do stare at all of the girls walk by.  I have even seen them stop what they are doing, and be quiet, just staring.  My roommates and I have just learned to ignore them completely though.  Parking is amazing.  I have never found any problem parking.  There are often guests who park inside, and there is still always parking.  The appliances you see in the pictures of the model are the EXACT same appliances you get, so there are no surprises.  Overall, this is great place to live.  Unfortunately, none of us resigned the lease because ther only bad thing is the cost.   It is way too expensive. Users/CaitlinBallengee

20090523 23:49:27 nbsp I dont know what everyone is complaining about. Ive never had any of the problems mentioned here.
Staff is wonderful. Theyve always been nice and responsive. The management feel like residents, not like a business. They all recognized me when just from when I originally came by to look at the model and have always been good to talk to, even if nothings up. Maintenance was pretty much perfect. All our requests were handled either the day of or the next day, During our stay we got a new stove (I didnt even notice anything wrong with the first one, but new one was nice), and two new sinks after reporting minor leaks. They always fix up anything we ask for, even if its not something major that theyre supposed to be responsible for. Management also accepts packages, which was very convenient.

The floor plan is very well thought out. Its spacious, open and wellproportioned. The living room is a good size, but not so large that the bedrooms are small for the square footage. The kitchen is also very roomy, and the vaulted ceilings make everything seem bigger. The upstairs bathroom (I live in a 4bed/2bath) is split into a sink/vanity area, a shower room and a toilet room, which makes having three girls up stairs quite manageable.

Parking is amazing. I have never seen the lot full, and most the time its more than half empty. Theres no need for assigned parking and having guests over is never a hassle. Around half the parking is covered, also.

Its mostly quiet. Occasionally there will be parties, but theyre never out of hand, too late or on weeknights. Lately the gardeners have been out earlier, but assuming youre getting up before 10 for class its not a problem. The school children are sometimes loud and a bit shrill, but theyre mostly only loud in the afternoon, so unless youre constantly napping with  your windows open (like me) then its no big deal.

The location is excellent. Across the street, maybe 3 minutes walk, from the M and W, and the P and Q are both just a few minutes further. Its across the street from SafeWay, OfficeMax, House of Chang Chinese food, Common grounds coffee, davis creamery ice cream, rite aid, round table pizza pizza, walnut park a huge park with tennis courts, volleyball, playgrounds... Basically its great, even if South Davis seems far away. Its also only a short walk to the post office and DMV.

Rents not that bad. For the 4 bedroom 2 bath starting fall 2007 it was $2440, +$30 for an inunit washer/dryer. It seems like a big number, but when looking for somewhere else to live (My roommates broke into factions. Wanting to live alone is my only reason to leave here.) I found nowhere comparable in terms of size, amenities, and general quality even only a little above the price here, much less below it. Security deposit was $1000, which spread amongst roommates wasnt a problem. Given my past interactions with the staff I dont foresee issues with getting the deposit back, though Ill report back after Ive moved.

Also, the official policy is no pets though Ive seen cats around as well as a couple small dogs in units. Either no pets is just the official policy to give management the upper hand in case of problems or people are just keeping them secretly.

I know this is quite long, but given my extreme satisfaction living here, I wanted to balance out the negative comments, especially those against the staff. The ease of dealing with them alone makes it that much more comfortable. Theyve always made it feel like a home for me, not like a business where theyre just hanging around to get my rent. The length of this comment is proportional to how much I like Sorrento.

tl;dr AAA+++ would live here again, full recommendation. Users/omgstarkitty

20091118 12:43:10 nbsp Reporting back re: security deposit:
They waived the cleaning charges, probably because needing light cleaning is normal wear and tear, but they charged $750, 3/4 of the deposit to repaint the entire apartment. However, they waived part of that bill and only charged us 614/1000. Perhaps they waived the price of what a touch up would have been and only charged us for excessive costs. They also did a full repainting on one of our neighbors apartments. Both were occupied by us for only one year. They said they charged for it because after one year the apartment should only need touch ups. I washed most of the walls, but one roommate spackled most of her room, so maybe they had to repaint that. So the security deposit refund wasnt great, but overall I found living there to be pleasant.

The only other downside I forgot to mention was that the fan in the shower room wasnt strong enough and the laundry isnt wellvented. Anytime the dryer was running the second floor would be hot and humid and windows would get condensation on them. There was a bit of mold, in the shower room and also on some downstairs windows.  Users/omgstarkitty

20091228 5:22:00 nbsp Moved in Fall 2009 and on lease till Summer 2010, but Im iffy about the management and place. Theyre usually polite, but sometimes they give the impression that my coming into the office (to pickup packages, report problems, etc.) is a bother. When something happens, they are very slow to admit that they are in the wrong. Case in point: My fiance was supposed to receive a package on Dec. 21st, and both USPS and Amazon reported it as delivered on that day. The office tells him they dont have it. Its most likely because of their lack of accurate records, because the package log is in disarray with info often omitted. The management called the post office to check, and its been delivered. They still suggest that my fiance might have somehow missed it (hes been checking the mail everyday and hes the only one doing so at our place). Or that the mailperson didnt drop it off. Thats not the only incident though. Theyve been unhelpful when we talked to them about our downstairs neighbors smoking 24/7 (we have asthmatics whose conditions are worsening because of the smoke). Theyve had issues recognizing that my fiance has signed the lease, and kept saying that he wasnt (someone didnt have all the records together). Things take more than a few days to get fixed, and like I said before, I dont always feel welcome going into the office for something.

Perhaps I have been spoiled by the great management at Tanglewood. They worked with us through many troubles and were sincerely friendly. Sorrento isnt bad, and as mentioned above, theyre generally polite, but combined with my experiences with the management, the nonexistent gym, and overall quality of the complex, I definitely will not be staying here after my obligatory year. Users/SatomiNakajima

             I apologize that your experiences with our office have not been as pleasant as I would like. We pride ourselves in having an          open door policy and I hope that we can make you feel comfortable to come in to our office to share any of your concerns or comments.

            We accept parcels and packages in our office as a courtesy and convenience to our residents.  I believe that you have             received your misplaced package but admit that the delay was caused by our confusion with the lease.  It is now clear to          us that while your fiancés given name is on the lease, it is not the name he goes by and we will not assume that it will be the   name on future packages.

            Thank you for bringing your concerns to our attention so that we can improve.

            Lori Brenner
            Community Manager

20100123 23:53:07 nbsp The place is very nice, lots of parking and very quiet.
I am looking to sublease my apartment from preferably March 20th ish til Aug 31 when the lease ends. How ever I am flexible. Contact for more info!
The apartment is a 2 bd 1.5 bath townhouse. Users/Sulli

20100318 11:42:11 nbsp THERE IS AN APT AVAILABLE! a 2bd 1.5 townhouse! :) available for rent june 2010 until aug 2011! or if youd prefer.. sept 2010 until aug 2011! let me know :) email me at nighteffectx8k@yahoo.com thanks! this is a great place to live btw Users/sallythatgirl

20100902 23:55:28 nbsp I LOVED LIVING HERE! although my apt my have been a burden on sorrento.. they worked with us! <3 maintenance is AMAZING! i walked into the office mentioning that there was mold in the bathroom (fan not strong enough) and then 5 min after i got back into my apartment.. someone was there! mold GONE! and has not been back!

i was really sad because i JUST got a car and it was VANDALIZED here :( $2000 worth of damage.. i only had my car for 3 days too! and NO ONE knew i had a car.. beware of jealous people :(

but i want to end on a good note because i loved living here a lot :) people in the office are really friendly! annnnd.. theyre actually really good with packages.. got 3 of em and i was able to pick them up right when i saw the truck drive in! Users/sallythatgirl

20110415 21:33:24 nbsp been living here for my whole college life and i have no complaints whatsoever about my experience here. the workers of the apartment complex are super awesome and friendly. 0 complaints Users/pol

20110711 12:25:52 nbsp I miss this place. They have spacious rooms and maintenance was always on top of things.  Users/mechee

20110808 17:01:50 nbsp I lived here for a year, and now that Im moving out this month, it makes me realize how much Im going to miss this place. :( I snagged a 3BR apartment by luck because that unit wanted to switch to a 2BR or something, and Im so glad I got to live here. Its also pretty obvious what a nice place Sorrento is since all of the rooms fill up stupidly fast during leasing season.

I think my favorite thing about this apartment complex is the management. Baker, Brian, and all of the others are incredibly nice and helpful; they also have fresh cookies and coffee everyday in the main office in case you feeling like dropping by before class or something. :) They have an online rentpaying system so that you dont even have to write them checks if you dont want to, they are really organized and on top of mailing/packages (they also have an automatic email notification so you dont have to call the office yourself every day and check), and their maintenance is extremely fast in responding to work orders/requests/apartment issues (within one or two days). Everyone who works there is extremely flexible and understanding.

One time I went home for a weekend, came back Sunday night, and realized that I forgot my keys inside. I was outside of my apartment with a bunch of crap I brought back from home, and I just called the office and left a message, even though the office was closed. I think it serves as some sort of an emergency thing and someone must have been notified of a new message, because 5 minutes later, someone came by and let me in with a key. Our apartment unit also had some roaches during a super rainy, wet couple of weeks and still knowing that it wasnt really Sorrentos fault that the dampness attracts bugs and stuff, they set up an appointment for a spraying and pest control like, the week of. I also had problems with leases and stuff, and while its not really their business to help you get through that (if you sign, its your responsibility), they still offered support and helped me get through it.

As for the apartments, I know that our unit was very clean and pretty new looking, but I havent been in any of the others. Everything else, including the exercise room and pool, is nice, and the apartment is super close to Safeway, Rite Aid, Office Max, etc. and also to the bus stops. Whats not to love about Sorrento? :) The prices are reasonable, management is friendly and cooperative, and the location is perfect. I feel like even if rent was a little more expensive it would still be worth it to live here, simply because they make your experience of living here that much better.
 Users/susiebear

20121106 13:12:52 nbsp I have been living in Sorrento for about three years. I really enjoy living here. The maintenance guys are great they fixed our washer/dryer and the bathtub. I live with five other girls our bathtub its always clog. The maintenance guys never complain. Last week the maintenance guys and Thomas the manager came door to door to just check see if any houses need to fix anything. They gave out icecream, too. It was so nice of them ) I love living here b/c the bus is convenient and there is a Safeway right across the street from us.  Users/AppleSun

20130901 18:05:31 nbsp I lived in Avalon for the past three years, and I only have great things to say about the amazing staff. They have never failed to respond to our problems on the day of our maintenance requests, and they even replaced our washer and refrigerator during those years. The staff has always been friendly, patient, and fun to talk to, both before and after the change in management. Only a few days ago, we moved across the street to their Sorrento properties. Although the staff was swamped with the demands of other residents, they helped us a ton by making sure that we would have plenty of time carry things over, clean our previous house, and not be homeless for any amount of time. They are simply awesome! Users/wildcard418

20130903 18:12:24 nbsp Life in these most beautiful apartment homes has been a pleasure.  Well designed, quality in every aspect of the word...carpeting, appliances, fireplace, fixtures.  Management addressees concerns almost immediately, known for maintaining the property as functionally and aesthetically as possible.  Quite and safe community, proximal /within walking distance to grocery store, restaurants, coffee houses, office supply store and bank.  Bus line M & W, as well as several Yolo buses stop near the complex on Cowell Blvd.  Youll like it here if you choose to become a resident.   Users/dyerike


text:
Zensoku Tuning is a new tuner orientated automotive parts retail company based in Davis, CA.  Currently offering a variety of HID (High Intensity Discharge) light kits, headlamp assemblies, and suspension parts.  Soon to be partnered up with http://www.whiteline.com.au/ Whiteline Automotive to provide more quality performance parts.

The owners of the company are both into motors ports and love their cars.  Eddie owns a black 2003 Subaru WRX and Joseph owns a red 2005 Honda S2000.  If you see us driving around feel free to come by and talk to us, we are friendly people; and we always hook up our friends!

20080608 20:34:07 nbsp Where are you guys located at? I need some HID! Users/StevenDaubert

20080630 01:20:46 nbsp Hey Steven!  Were currently doing business through the internet and looking for a more permanent location.  Good news is that well deliver it to your door if you live in Davis!.  Visit our website for more information: www.zensokutuning.com Users/Eddie

20080630 01:21:29 nbsp Were doing special introductory pricing right now so please email me for MORE INFORMATION! Users/Eddie

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c     Subroutine to get vectors of data
c     AJ_Kettle, 16May2019
c     27Nov2019: subroutine modified for qff program
c     11Mar2020: used for USAF update program

      SUBROUTINE get_stnconfig_qff(s_directory_stnconfig,
     +  s_filename_stnconfig_new,
     +  l_scoutput_rgh,l_scoutput,ilen,
     +  l_numfield,s_vec_header,s_arch_fields)

      IMPLICIT NONE
c************************************************************************
c     Variables passed into program

      CHARACTER(LEN=300)  :: s_directory_stnconfig
      CHARACTER(LEN=300)  :: s_filename_stnconfig_new

      INTEGER             :: l_scoutput_rgh
      INTEGER             :: l_scoutput
      INTEGER             :: ilen

      CHARACTER(LEN=1000) :: s_header

      INTEGER             :: l_numfield
      CHARACTER(LEN=ilen) :: s_vec_header(50)
      CHARACTER(LEN=ilen) :: s_arch_fields(l_scoutput_rgh,50)
c*****
c     Variables used within program

      INTEGER             :: i,j,k,ii,jj,kk
      INTEGER             :: io

      CHARACTER(LEN=1000) :: s_pathandname
      CHARACTER(LEN=1000) :: s_linget

      CHARACTER(LEN=ilen)  :: s_vec_fields(50)

      INTEGER             :: i_len_header
      INTEGER             :: i_len_line
c************************************************************************
      print*,'just entered get_stnconfig_qff'

      s_pathandname=
     +  TRIM(s_directory_stnconfig)//TRIM(s_filename_stnconfig_new)

      ii=0
      jj=0

      OPEN(UNIT=5,
     +  FILE=TRIM(s_pathandname),
     +  FORM='formatted',
     +  STATUS='OLD',ACTION='READ')

      READ(5,1000,IOSTAT=io) s_header

      i_len_header=LEN_TRIM(s_header)

c     Call subroutine to get pipe-separated fields
      CALL get_elements_stnconfig_line(s_header,
     +  l_numfield,s_vec_header)

c      print*,'s_header=',TRIM(s_header)
c      print*,'i_len_header=',i_len_header
c      STOP 'get_stnconfig_write'

      DO
       READ(5,1000,IOSTAT=io) s_linget
 1000  FORMAT(a1000)

       IF (io .GT. 0) THEN
        WRITE(*,*) 'Check input. Something went wrong'
        GOTO 100
       ELSE IF (io .LT. 0) THEN
        WRITE(*,*) 'end of file reached'
        GOTO 100
       ELSE
        ii=ii+1

c       Crash condition if too many lines
        IF (ii.GT.l_scoutput_rgh) THEN
         print*,'emergency stop; stnconfig data lines exceeds limit'
         STOP 'get_stnconfig_qff'
        ENDIF

        i_len_line=LEN_TRIM(s_linget)

c       Call subroutine to get pipe-separated fields
        CALL get_elements_stnconfig_line(s_linget,
     +    l_numfield,s_vec_fields)

c       Archive fields
        DO j=1,l_numfield
         s_arch_fields(ii,j)=s_vec_fields(j)
        ENDDO

        IF (i_len_line.GT.1000) THEN
         print*,'emergency stop line length over range'
         STOP 'get_stnconfig_write'
        ENDIF

       ENDIF
      ENDDO
 100  CONTINUE
      CLOSE(UNIT=5)

      l_scoutput=ii

      print*,'l_numfield=',l_numfield
      print*,'s_vec_header=',(TRIM(s_vec_header(i))//'|',i=1,l_numfield)
      print*,'ii,l_scoutput=',ii,l_scoutput
      print*,'l_scoutput_rgh=',l_scoutput_rgh

c      DO i=1,l_numfield
c       print*,i,TRIM(s_vec_header(i))
c      ENDDO

      print*,'just leaving get_stnconfig_qff'

c      STOP 'get_stnconfig_qff'

      RETURN
      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


module moverlap
use constants
implicit none

contains

function overlap(n1,i1,j1,a1,k1,n2,i2,j2,a2,k2) result(ans)
 integer, intent(in) :: n1,i1,j1,n2,i2,j2
real(quad), intent(in) :: a1,a2,k1,k2
real(quad) :: ans
! - - - - - - - - - - - - - - - - - -
integer :: key
real(quad) :: ap, spp
real(quad) :: p1,p2,p3,p4,p5,p6,p7,p8,p9,p10
real(quad) :: e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12
! - - - - - - - - - - - - - - - - - -

key = n1 + 2*i1 + 4*j1 + 8*n2 + 16*i2 + 32*j2

select case (key)
  case (0)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    ans = spp
  case (1)
    ans = 0
  case (2)
    ans = 0
  case (3)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p3 = k1**2
    e2 = exp(-(p1*p3)/4.)
    ans = (e2*k1*p1*spp)/2.
  case (4)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p3 = k1**2
    e2 = exp(-(p1*p3)/4.)
    ans = e2*spp
  case (5)
    ans = 0
  case (6)
    ans = 0
  case (7)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p3 = k1**2
    e12 = exp(-(p1*p3))
    ans = (e12*k1*p1*spp)/2.
  case (8)
    ans = 0
  case (9)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    ans = (p1*spp)/2.
  case (10)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p3 = k1**2
    e2 = exp(-(p1*p3)/4.)
    ans = (e2*k1*p1*spp)/2.
  case (11)
    ans = 0
  case (12)
    ans = 0
  case (13)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p3 = k1**2
    p10 = ap**(-2)
    e2 = exp(-(p1*p3)/4.)
    ans = (e2*p10*(2*ap - p3)*spp)/4.
  case (14)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p3 = k1**2
    e12 = exp(-(p1*p3))
    ans = (e12*k1*p1*spp)/2.
  case (15)
    ans = 0
  case (16)
    ans = 0
  case (17)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p2 = k2**2
    e1 = exp(-(p1*p2)/4.)
    ans = (e1*k2*p1*spp)/2.
  case (18)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = ((e3 - e4)*spp)/2.
  case (19)
    ans = 0
  case (20)
    ans = 0
  case (21)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = ((e3*(-k1 + k2) + e4*(k1 + k2))*p1*spp)/4.
  case (22)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p8 = (-2*k1 + k2)**2
    p9 = (2*k1 + k2)**2
    e7 = exp(-(p1*p8)/4.)
    e8 = exp(-(p1*p9)/4.)
    ans = ((e7 - e8)*spp)/4.
  case (23)
    ans = 0
  case (24)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p2 = k2**2
    e1 = exp(-(p1*p2)/4.)
    ans = (e1*k2*p1*spp)/2.
  case (25)
    ans = 0
  case (26)
    ans = 0
  case (27)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    p10 = ap**(-2)
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = (p10*(e3*(2*ap - p4) + e4*(-2*ap + p5))*spp)/8.
  case (28)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = ((e3*(-k1 + k2) + e4*(k1 + k2))*p1*spp)/4.
  case (29)
    ans = 0
  case (30)
    ans = 0
  case (31)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p8 = (-2*k1 + k2)**2
    p9 = (2*k1 + k2)**2
    p10 = ap**(-2)
    e7 = exp(-(p1*p8)/4.)
    e8 = exp(-(p1*p9)/4.)
    ans = (p10*(e7*(2*ap - p8) + e8*(-2*ap + p9))*spp)/16.
  case (32)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p2 = k2**2
    e1 = exp(-(p1*p2)/4.)
    ans = e1*spp
  case (33)
    ans = 0
  case (34)
    ans = 0
  case (35)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = ((e3*(k1 - k2) + e4*(k1 + k2))*p1*spp)/4.
  case (36)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = ((e3 + e4)*spp)/2.
  case (37)
    ans = 0
  case (38)
    ans = 0
  case (39)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p8 = (-2*k1 + k2)**2
    p9 = (2*k1 + k2)**2
    e7 = exp(-(p1*p8)/4.)
    e8 = exp(-(p1*p9)/4.)
    ans = ((2*(e7 + e8)*k1 + (-e7 + e8)*k2)*p1*spp)/8.
  case (40)
    ans = 0
  case (41)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p2 = k2**2
    p10 = ap**(-2)
    e1 = exp(-(p1*p2)/4.)
    ans = (e1*p10*(2*ap - p2)*spp)/4.
  case (42)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = ((e3*(k1 - k2) + e4*(k1 + k2))*p1*spp)/4.
  case (43)
    ans = 0
  case (44)
    ans = 0
  case (45)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    p10 = ap**(-2)
    e3 = exp(-(p1*p4)/4.)
    e4 = exp(-(p1*p5)/4.)
    ans = (p10*(2*ap*(e3 + e4) - e3*p4 - e4*p5)*spp)/8.
  case (46)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p8 = (-2*k1 + k2)**2
    p9 = (2*k1 + k2)**2
    e7 = exp(-(p1*p8)/4.)
    e8 = exp(-(p1*p9)/4.)
    ans = ((2*(e7 + e8)*k1 + (-e7 + e8)*k2)*p1*spp)/8.
  case (47)
    ans = 0
  case (48)
    ans = 0
  case (49)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p2 = k2**2
    e11 = exp(-(p1*p2))
    ans = (e11*k2*p1*spp)/2.
  case (50)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p6 = (k1 - 2*k2)**2
    p7 = (k1 + 2*k2)**2
    e5 = exp(-(p1*p6)/4.)
    e6 = exp(-(p1*p7)/4.)
    ans = ((e5 - e6)*spp)/4.
  case (51)
    ans = 0
  case (52)
    ans = 0
  case (53)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p6 = (k1 - 2*k2)**2
    p7 = (k1 + 2*k2)**2
    e5 = exp(-(p1*p6)/4.)
    e6 = exp(-(p1*p7)/4.)
    ans = (((-e5 + e6)*k1 + 2*(e5 + e6)*k2)*p1*spp)/8.
  case (54)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    e9 = exp(-(p1*p4))
    e10 = exp(-(p1*p5))
    ans = ((-e10 + e9)*spp)/8.
  case (55)
    ans = 0
  case (56)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p2 = k2**2
    e11 = exp(-(p1*p2))
    ans = (e11*k2*p1*spp)/2.
  case (57)
    ans = 0
  case (58)
    ans = 0
  case (59)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p6 = (k1 - 2*k2)**2
    p7 = (k1 + 2*k2)**2
    p10 = ap**(-2)
    e5 = exp(-(p1*p6)/4.)
    e6 = exp(-(p1*p7)/4.)
    ans = (p10*(e5*(2*ap - p6) + e6*(-2*ap + p7))*spp)/16.
  case (60)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p6 = (k1 - 2*k2)**2
    p7 = (k1 + 2*k2)**2
    e5 = exp(-(p1*p6)/4.)
    e6 = exp(-(p1*p7)/4.)
    ans = (((-e5 + e6)*k1 + 2*(e5 + e6)*k2)*p1*spp)/8.
  case (61)
    ans = 0
  case (62)
    ans = 0
  case (63)
    ap = a1 + a2
    spp = sqrt(q_pi/ap)
    p1 = 1/ap
    p4 = (k1 - k2)**2
    p5 = (k1 + k2)**2
    p10 = ap**(-2)
    e9 = exp(-(p1*p4))
    e10 = exp(-(p1*p5))
    ans = (p10*(e9*(ap - 2*p4) - e10*(ap - 2*p5))*spp)/16.
end select
end function overlap
end module moverlap


text:
      SUBROUTINE ZSTEGR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
     $                   M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK,
     $                   LIWORK, INFO )
*
*  -- LAPACK computational routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1999
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, RANGE
      INTEGER            IL, INFO, IU, LDZ, LIWORK, LWORK, M, N
      DOUBLE PRECISION   ABSTOL, VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            ISUPPZ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
      COMPLEX*16         Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
* ZSTEGR computes selected eigenvalues and, optionally, eigenvectors
* of a real symmetric tridiagonal matrix T.  Eigenvalues and
* eigenvectors can be selected by specifying either a range of values
* or a range of indices for the desired eigenvalues. The eigenvalues
* are computed by the dqds algorithm, while orthogonal eigenvectors are
* computed from various ``good'' L D L^T representations (also known as
* Relatively Robust Representations). Gram-Schmidt orthogonalization is
* avoided as far as possible. More specifically, the various steps of
* the algorithm are as follows. For the i-th unreduced block of T,
*     (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T
*         is a relatively robust representation,
*     (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high
*         relative accuracy by the dqds algorithm,
*     (c) If there is a cluster of close eigenvalues, "choose" sigma_i
*         close to the cluster, and go to step (a),
*     (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T,
*         compute the corresponding eigenvector by forming a
*         rank-revealing twisted factorization.
*  The desired accuracy of the output can be specified by the input
*  parameter ABSTOL.
*
*  For more details, see "A new O(n^2) algorithm for the symmetric
*  tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon,
*  Computer Science Division Technical Report No. UCB/CSD-97-971,
*  UC Berkeley, May 1997.
*
*  Note 1 : Currently ZSTEGR is only set up to find ALL the n
*  eigenvalues and eigenvectors of T in O(n^2) time
*  Note 2 : Currently the routine ZSTEIN is called when an appropriate
*  sigma_i cannot be chosen in step (c) above. ZSTEIN invokes modified
*  Gram-Schmidt when eigenvalues are close.
*  Note 3 : ZSTEGR works only on machines which follow ieee-754
*  floating-point standard in their handling of infinities and NaNs.
*  Normal execution of ZSTEGR may create NaNs and infinities and hence
*  may abort due to a floating point exception in environments which
*  do not conform to the ieee standard.
*
*  Arguments
*  =========
*
*  JOBZ    (input) CHARACTER*1
*          = 'N':  Compute eigenvalues only;
*          = 'V':  Compute eigenvalues and eigenvectors.
*
*  RANGE   (input) CHARACTER*1
*          = 'A': all eigenvalues will be found.
*          = 'V': all eigenvalues in the half-open interval (VL,VU]
*                 will be found.
*          = 'I': the IL-th through IU-th eigenvalues will be found.
********** Only RANGE = 'A' is currently supported *********************
*
*  N       (input) INTEGER
*          The order of the matrix.  N >= 0.
*
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the n diagonal elements of the tridiagonal matrix
*          T. On exit, D is overwritten.
*
*  E       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the (n-1) subdiagonal elements of the tridiagonal
*          matrix T in elements 1 to N-1 of E; E(N) need not be set.
*          On exit, E is overwritten.
*
*  VL      (input) DOUBLE PRECISION
*  VU      (input) DOUBLE PRECISION
*          If RANGE='V', the lower and upper bounds of the interval to
*          be searched for eigenvalues. VL < VU.
*          Not referenced if RANGE = 'A' or 'I'.
*
*  IL      (input) INTEGER
*  IU      (input) INTEGER
*          If RANGE='I', the indices (in ascending order) of the
*          smallest and largest eigenvalues to be returned.
*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
*          Not referenced if RANGE = 'A' or 'V'.
*
*  ABSTOL  (input) DOUBLE PRECISION
*          The absolute error tolerance for the
*          eigenvalues/eigenvectors. IF JOBZ = 'V', the eigenvalues and
*          eigenvectors output have residual norms bounded by ABSTOL,
*          and the dot products between different eigenvectors are
*          bounded by ABSTOL. If ABSTOL is less than N*EPS*|T|, then
*          N*EPS*|T| will be used in its place, where EPS is the
*          machine precision and |T| is the 1-norm of the tridiagonal
*          matrix. The eigenvalues are computed to an accuracy of
*          EPS*|T| irrespective of ABSTOL. If high relative accuracy
*          is important, set ABSTOL to DLAMCH( 'Safe minimum' ).
*          See Barlow and Demmel "Computing Accurate Eigensystems of
*          Scaled Diagonally Dominant Matrices", LAPACK Working Note #7
*          for a discussion of which matrices define their eigenvalues
*          to high relative accuracy.
*
*  M       (output) INTEGER
*          The total number of eigenvalues found.  0 <= M <= N.
*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
*
*  W       (output) DOUBLE PRECISION array, dimension (N)
*          The first M elements contain the selected eigenvalues in
*          ascending order.
*
*  Z       (output) COMPLEX*16 array, dimension (LDZ, max(1,M) )
*          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
*          contain the orthonormal eigenvectors of the matrix T
*          corresponding to the selected eigenvalues, with the i-th
*          column of Z holding the eigenvector associated with W(i).
*          If JOBZ = 'N', then Z is not referenced.
*          Note: the user must ensure that at least max(1,M) columns are
*          supplied in the array Z; if RANGE = 'V', the exact value of M
*          is not known in advance and an upper bound must be used.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.  LDZ >= 1, and if
*          JOBZ = 'V', LDZ >= max(1,N).
*
*  ISUPPZ  (output) INTEGER ARRAY, dimension ( 2*max(1,M) )
*          The support of the eigenvectors in Z, i.e., the indices
*          indicating the nonzero elements in Z. The i-th eigenvector
*          is nonzero only in elements ISUPPZ( 2*i-1 ) through
*          ISUPPZ( 2*i ).
*
*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal
*          (and minimal) LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= max(1,18*N)
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  IWORK   (workspace/output) INTEGER array, dimension (LIWORK)
*          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
*
*  LIWORK  (input) INTEGER
*          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
*
*          If LIWORK = -1, then a workspace query is assumed; the
*          routine only calculates the optimal size of the IWORK array,
*          returns this value as the first entry of the IWORK array, and
*          no error message related to LIWORK is issued by XERBLA.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = 1, internal error in DLARRE,
*                if INFO = 2, internal error in ZLARRV.
*
*  Further Details
*  ===============
*
*  Based on contributions by
*     Inderjit Dhillon, IBM Almaden, USA
*     Osni Marques, LBNL/NERSC, USA
*     Ken Stanley, Computer Science Division, University of
*       California at Berkeley, USA
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
      COMPLEX*16         CZERO
      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ
      INTEGER            I, IBEGIN, IEND, IINDBL, IINDWK, IINFO, IINSPL,
     $                   INDGRS, INDWOF, INDWRK, ITMP, J, JJ, LIWMIN,
     $                   LWMIN, NSPLIT
      DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SCALE, SMLNUM,
     $                   THRESH, TMP, TNRM, TOL
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           LSAME, DLAMCH, DLANST
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLARRE, DSCAL, XERBLA, ZLARRV, ZLASET, ZSWAP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DBLE, MAX, MIN, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      WANTZ = LSAME( JOBZ, 'V' )
      ALLEIG = LSAME( RANGE, 'A' )
      VALEIG = LSAME( RANGE, 'V' )
      INDEIG = LSAME( RANGE, 'I' )
*
      LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) )
      LWMIN = 18*N
      LIWMIN = 10*N
*
      INFO = 0
      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
         INFO = -2
*
*     The following two lines need to be removed once the
*     RANGE = 'V' and RANGE = 'I' options are provided.
*
      ELSE IF( VALEIG .OR. INDEIG ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( VALEIG .AND. N.GT.0 .AND. VU.LE.VL ) THEN
         INFO = -7
      ELSE IF( INDEIG .AND. IL.LT.1 ) THEN
         INFO = -8
*     The following change should be made in DSTEVX also, otherwise
*     IL can be specified as N+1 and IU as N.
*     ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN
      ELSE IF( INDEIG .AND. ( IU.LT.IL .OR. IU.GT.N ) ) THEN
         INFO = -9
      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
         INFO = -14
      ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -17
      ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
         INFO = -19
      END IF
      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = LWMIN
         IWORK( 1 ) = LIWMIN
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'ZSTEGR', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      M = 0
      IF( N.EQ.0 )
     $   RETURN
*
      IF( N.EQ.1 ) THEN
         IF( ALLEIG .OR. INDEIG ) THEN
            M = 1
            W( 1 ) = D( 1 )
         ELSE
            IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
               M = 1
               W( 1 ) = D( 1 )
            END IF
         END IF
         IF( WANTZ )
     $      Z( 1, 1 ) = ONE
         RETURN
      END IF
*
*     Get machine constants.
*
      SAFMIN = DLAMCH( 'Safe minimum' )
      EPS = DLAMCH( 'Precision' )
      SMLNUM = SAFMIN / EPS
      BIGNUM = ONE / SMLNUM
      RMIN = SQRT( SMLNUM )
      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
*
*     Scale matrix to allowable range, if necessary.
*
      SCALE = ONE
      TNRM = DLANST( 'M', N, D, E )
      IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
         SCALE = RMIN / TNRM
      ELSE IF( TNRM.GT.RMAX ) THEN
         SCALE = RMAX / TNRM
      END IF
      IF( SCALE.NE.ONE ) THEN
         CALL DSCAL( N, SCALE, D, 1 )
         CALL DSCAL( N-1, SCALE, E, 1 )
         TNRM = TNRM*SCALE
      END IF
      INDGRS = 1
      INDWOF = 2*N + 1
      INDWRK = 3*N + 1
*
      IINSPL = 1
      IINDBL = N + 1
      IINDWK = 2*N + 1
*
      CALL ZLASET( 'Full', N, N, CZERO, CZERO, Z, LDZ )
*
*     Compute the desired eigenvalues of the tridiagonal after splitting
*     into smaller subblocks if the corresponding of-diagonal elements
*     are small
*
      THRESH = EPS*TNRM
      CALL DLARRE( N, D, E, THRESH, NSPLIT, IWORK( IINSPL ), M, W,
     $             WORK( INDWOF ), WORK( INDGRS ), WORK( INDWRK ),
     $             IINFO )
      IF( IINFO.NE.0 ) THEN
         INFO = 1
         RETURN
      END IF
*
      IF( WANTZ ) THEN
*
*        Compute the desired eigenvectors corresponding to the computed
*        eigenvalues
*
         TOL = MAX( ABSTOL, DBLE( N )*THRESH )
         IBEGIN = 1
         DO 20 I = 1, NSPLIT
            IEND = IWORK( IINSPL+I-1 )
            DO 10 J = IBEGIN, IEND
               IWORK( IINDBL+J-1 ) = I
   10       CONTINUE
            IBEGIN = IEND + 1
   20    CONTINUE
*
         CALL ZLARRV( N, D, E, IWORK( IINSPL ), M, W, IWORK( IINDBL ),
     $                WORK( INDGRS ), TOL, Z, LDZ, ISUPPZ,
     $                WORK( INDWRK ), IWORK( IINDWK ), IINFO )
         IF( IINFO.NE.0 ) THEN
            INFO = 2
            RETURN
         END IF
*
      END IF
*
      IBEGIN = 1
      DO 40 I = 1, NSPLIT
         IEND = IWORK( IINSPL+I-1 )
         DO 30 J = IBEGIN, IEND
            W( J ) = W( J ) + WORK( INDWOF+I-1 )
   30    CONTINUE
         IBEGIN = IEND + 1
   40 CONTINUE
*
*     If matrix was scaled, then rescale eigenvalues appropriately.
*
      IF( SCALE.NE.ONE ) THEN
         CALL DSCAL( M, ONE / SCALE, W, 1 )
      END IF
*
*     If eigenvalues are not in order, then sort them, along with
*     eigenvectors.
*
      IF( NSPLIT.GT.1 ) THEN
         DO 60 J = 1, M - 1
            I = 0
            TMP = W( J )
            DO 50 JJ = J + 1, M
               IF( W( JJ ).LT.TMP ) THEN
                  I = JJ
                  TMP = W( JJ )
               END IF
   50       CONTINUE
            IF( I.NE.0 ) THEN
               W( I ) = W( J )
               W( J ) = TMP
               IF( WANTZ ) THEN
                  CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
                  ITMP = ISUPPZ( 2*I-1 )
                  ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 )
                  ISUPPZ( 2*J-1 ) = ITMP
                  ITMP = ISUPPZ( 2*I )
                  ISUPPZ( 2*I ) = ISUPPZ( 2*J )
                  ISUPPZ( 2*J ) = ITMP
               END IF
            END IF
   60    CONTINUE
      END IF
*
      WORK( 1 ) = LWMIN
      IWORK( 1 ) = LIWMIN
      RETURN
*
*     End of ZSTEGR
*
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:59:12 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE GET_EDGE__genmod
          INTERFACE
            SUBROUTINE GET_EDGE(NGNOD,N,ITYPE,E1,E2,IXMIN,IXMAX,IZMIN,  &
     &IZMAX,SENS)
              INTEGER(KIND=4), INTENT(IN) :: NGNOD
              INTEGER(KIND=4), INTENT(IN) :: N(NGNOD)
              INTEGER(KIND=4), INTENT(IN) :: ITYPE
              INTEGER(KIND=4), INTENT(IN) :: E1
              INTEGER(KIND=4), INTENT(IN) :: E2
              INTEGER(KIND=4), INTENT(OUT) :: IXMIN
              INTEGER(KIND=4), INTENT(OUT) :: IXMAX
              INTEGER(KIND=4), INTENT(OUT) :: IZMIN
              INTEGER(KIND=4), INTENT(OUT) :: IZMAX
              INTEGER(KIND=4), INTENT(OUT) :: SENS
            END SUBROUTINE GET_EDGE
          END INTERFACE
        END MODULE GET_EDGE__genmod


text:
!>-----------------------------------------
!! Main Program
!!
!! Initialize options and memory in init_model
!! Read initial conditions in bc_init (from a restart file if requested)
!! initialize physics packages in init_physics (e.g. tiedke and thompson if used)
!! If this run is a restart run, then set start to the restart timestep
!!      in otherwords, ntimesteps is the number of BC updates from the beginning of the entire model
!!      run, not just from the begining of this restart run
!! calculate model time in seconds based on the time between BC updates (in_dt)
!! Calculate the next model output time from current model time + output time delta (out_dt)
!!
!! Finally, loop until ntimesteps are reached updating boundary conditions and stepping the model forward
!!
!!  @author
!!  Ethan Gutmann (gutmann@ucar.edu)
!!
!!-----------------------------------------
program icar

    use iso_fortran_env,    only : output_unit
    use options_interface,  only : options_t
    use domain_interface,   only : domain_t
    use boundary_interface, only : boundary_t
    use output_interface,   only : output_t
    use time_step,          only : step                               ! Advance the model forward in time
    use initialization,     only : init_model
    use timer_interface,    only : timer_t
    use time_object,        only : Time_type
    use time_delta_object,  only : time_delta_t
    use wind,               only : update_winds
    use restart_interface,  only : restart_model
    use icar_constants,     only : kVARS

    implicit none

    type(options_t) :: options
    type(domain_t)  :: domain
    type(boundary_t):: boundary
    type(output_t)  :: dataset
    ! type(output_t)  :: surface_dataset
    type(timer_t)   :: initialization_timer, total_timer, input_timer, output_timer, physics_timer
    type(Time_type) :: next_output
    type(time_delta_t) :: small_time_delta

    character(len=1024) :: file_name
    character(len=49)   :: file_date_format = '(I4,"-",I0.2,"-",I0.2,"_",I0.2,"-",I0.2,"-",I0.2)'
    integer :: i
    integer :: output_vars(18)

    call small_time_delta%set(1)

    call total_timer%start()
    call initialization_timer%start()
    !-----------------------------------------
    !  Model Initialization
    !
    ! Reads config options and initializes domain and boundary conditions
    call init_model(options, domain, boundary)

    if (this_image()==1) write(*,*) "Setting up output files"
    ! should be combined into a single setup_output call
    call dataset%set_domain(domain)
    call dataset%add_variables(options%vars_for_restart, domain)

    output_vars = [kVARS%precipitation, kVARS%snowfall, kVARS%graupel, kVARS%cloud_fraction, kVARS%shortwave, kVARS%longwave, &
                   kVARS%sensible_heat, kVARS%latent_heat, kVARS%u_10m, kVARS%v_10m, kVARS%temperature_2m, kVARS%humidity_2m, &
                   kVARS%surface_pressure, kVARS%soil_totalmoisture, kVARS%snow_water_equivalent, kVARS%skin_temperature,     &
                   kVARS%latitude, kVARS%longitude]
    do i=1,size(output_vars)
        options%vars_for_output(output_vars)=1
    enddo
    ! call surface_dataset%set_domain(domain)
    ! call surface_dataset%add_variables(options%vars_for_output, domain)

    if (options%parameters%restart) then
        if (this_image()==1) write(*,*) "Reading restart data"
        call restart_model(domain, dataset, options)
    endif

    !-----------------------------------------
    !-----------------------------------------
    !  Time Loop
    !
    !   note that a timestep here is a forcing input timestep O(1-3hr), not a physics timestep O(20-100s)
    write(file_name, '(A,I6.6,"_",A,".nc")')      &
            trim(options%parameters%output_file), &
            this_image(),                         &
            trim(domain%model_time%as_string(file_date_format))

    call initialization_timer%stop()

    i=1
    call output_timer%start()
    call dataset%save_file(trim(file_name), i, domain%model_time)
    ! call surface_dataset%save_file("surface_"//trim(file_name), i, domain%model_time)
    next_output = domain%model_time + options%parameters%output_dt
    call output_timer%stop()
    i = i + 1

    if (this_image()==1) write(*,*) "Initialization complete, beginning physics integration."

    do while (domain%model_time < options%parameters%end_time)
        call flush(output_unit)

        ! -----------------------------------------------------
        !
        !  Read input data if necessary
        !
        ! -----------------------------------------------------
        call input_timer%start()
        if (boundary%current_time <= (domain%model_time + small_time_delta) ) then
            if (this_image()==1) write(*,*) ""
            if (this_image()==1) write(*,*) " ----------------------------------------------------------------------"
            if (this_image()==1) write(*,*) "Updating Boundary conditions"
            call boundary%update_forcing(options)
            call domain%interpolate_forcing(boundary, update=.True.)
            call update_winds(domain, options)

            ! Make the boundary condition dXdt values into units of [X]/s
            call domain%update_delta_fields(boundary%current_time - domain%model_time)
        endif
        call input_timer%stop()

        ! -----------------------------------------------------
        !
        !  Integrate physics forward in time
        !
        ! -----------------------------------------------------
        if (this_image()==1) write(*,*) "Running Physics"
        if (this_image()==1) write(*,*) "  Model time = ", trim(domain%model_time%as_string())
        if (this_image()==1) write(*,*) "   End  time = ", trim(options%parameters%end_time%as_string())
        if (this_image()==1) write(*,*) "  Next Input = ", trim(boundary%current_time%as_string())
        if (this_image()==1) write(*,*) "  Next Output= ", trim(next_output%as_string())

        ! this is the meat of the model physics, run all the physics for the current time step looping over internal timesteps
        call physics_timer%start()
        call step(domain, step_end(boundary%current_time, next_output), options)
        call physics_timer%stop()

        ! -----------------------------------------------------
        !
        !  Write output data if it is time
        !
        ! -----------------------------------------------------
        ! This is a bit of a hack until the output object is set up better to handle files with specified number of steps per file (or months or...)
        ! ideally this will just become timer_start, save_file()...
        ! the output object needs a pointer to model_time and will have to know how to create output file names
        call output_timer%start()
        if ((domain%model_time + small_time_delta) >= next_output) then
            if (this_image()==1) write(*,*) "Writing output file"
            if (i>24) then
                write(file_name, '(A,I6.6,"_",A,".nc")')    &
                    trim(options%parameters%output_file),   &
                    this_image(),                           &
                    trim(domain%model_time%as_string(file_date_format))
                i = 1
            endif

            call dataset%save_file(trim(file_name), i, next_output)
            ! call surface_dataset%save_file("surface_"//trim(file_name), i, domain%model_time)

            next_output = next_output + options%parameters%output_dt

            i = i + 1
        endif

        call output_timer%stop()

    end do
    !
    !-----------------------------------------
    call total_timer%stop()

    if (this_image()==1) then
        write(*,*) ""
        write(*,*) "Model run from : ",trim(options%parameters%start_time%as_string())
        write(*,*) "           to  : ",trim(options%parameters%end_time%as_string())
        write(*,*) "Domain : ",trim(options%parameters%init_conditions_file)
        write(*,*) "Number of images:",num_images()
        write(*,*) ""
        write(*,*) "First image timing:"
        write(*,*) "total   : ", trim(total_timer%as_string())
        write(*,*) "init    : ", trim(initialization_timer%as_string())
        write(*,*) "input   : ", trim(input_timer%as_string())
        write(*,*) "output  : ", trim(output_timer%as_string())
        write(*,*) "physics : ", trim(physics_timer%as_string())
    endif

contains

    function step_end(time1, time2) result(min_time)
        implicit none
        type(Time_type), intent(in) :: time1
        type(Time_type), intent(in) :: time2
        type(Time_type) :: min_time

        if (time1 <= time2 ) then
            min_time = time1
        else
            min_time = time2
        endif

    end function

end program

! This is the Doxygen mainpage documentation.  This should be moved to another file at some point.

!>------------------------------------------
!!  @mainpage
!!
!!  @section Introduction
!!  ICAR is a simplified atmospheric model designed primarily for climate downscaling, atmospheric sensitivity tests,
!!  and hopefully educational uses. At this early stage, the model is still undergoing rapid development, and users
!!  are encouraged to get updates frequently.
!!
!!  @section Running_ICAR
!!  To run the model 3D time-varying atmospheric data are required, though an ideal test case can be generated for
!!  simple simulations as well. There are some sample python scripts to help make input forcing files, but the WRF
!!  pre-processing system can also be used. Low-resolution WRF output files can be used directly, various reanalysis
!!  and GCM output files can be used with minimal pre-processing (just get all the variables in the same netcdf file.)
!!  In addition, a high-resolution netCDF topography file is required. This will define the grid that ICAR will run on.
!!  Finally an ICAR options file is used to specify various parameters for the model. A sample options file is provided
!!  in the run/ directory.
!!
!!  @section Developing
!!  This document provides the primary API and code structure documentation. The code is based on github.com/NCAR/icar
!!  Developers are encouraged to fork the main git repository and maintain their own git repository from which to
!!  issue pull requests.
!!
!!  @section Reference
!!  Gutmann, E. D., I. Barstad, M. P. Clark, J. R. Arnold, and R. M. Rasmussen (2016),
!!  The Intermediate Complexity Atmospheric Research Model, J. Hydrometeor, doi:10.1175/JHM-D-15-0155.1.
!!
!!------------------------------------------

max_stars_repo_name: d-reynolds/HICAR

max_issues_repo_name: d-reynolds/HICAR

max_forks_repo_name: d-reynolds/HICAR


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine rscmp1                       *
c     *                                                              *
c     *                       written by : bh                        *
c     *                                                              *
c     *                   last modified : 06/30/91                   *
c     *                                                              *
c     *     this subroutine computes the polar decompostion of the   *
c     *     deformation gradient into the right stretch tensor.      *
c     *     the computations are for a gauss point of a q3disop      *
c     *     or l3disop element.                                      *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine rscmp1( span, f, u )
      implicit integer (a-z)
      include 'param_def'
c
c                      parameter declarations
c
      double precision
     & f(mxvl,ndim,*), u(mxvl,*)
c
c                      locally allocated arrays & constants
c
      double precision
     &  c(mxvl,nstr), cc(mxvl,nstr), iu(mxvl), iiu(mxvl),
     &  iiiu(mxvl), a1(mxvl), b1(mxvl), c1(mxvl), one
      data one / 1.0 /
c
c                       u is in symmetric upper triangular form.
c
c
c                       compute the invariants of the right
c                       stretch tensor, the metric tensor, and
c                       its square.
c
      call ivcmp1( span, f, c, cc, iu, iiu, iiiu )
c
c                       compute multipliers.
c
      do i = 1, span
         a1(i)= one/(iu(i)*iiu(i)-iiiu(i))
         b1(i)= iu(i)*iiiu(i)
         c1(i)= iu(i)*iu(i)-iiu(i)
      end do
c
c                       compute the right stretch tensor.
c
      do i = 1, span
         u(i,1)= a1(i) * ( b1(i) + c1(i)*c(i,1) - cc(i,1) )
         u(i,2)= a1(i) * (         c1(i)*c(i,2) - cc(i,2) )
         u(i,3)= a1(i) * ( b1(i) + c1(i)*c(i,3) - cc(i,3) )
         u(i,4)= a1(i) * (         c1(i)*c(i,4) - cc(i,4) )
         u(i,5)= a1(i) * (         c1(i)*c(i,5) - cc(i,5) )
         u(i,6)= a1(i) * ( b1(i) + c1(i)*c(i,6) - cc(i,6) )
      end do
c
      return
      end


text:
subroutine lev_CutElementsAndListByLayers(a)
!DESCRIPTION
!   This subroutine get the element status: cut=1 and non cut=0
!-----------------------------------------------------------------------
   use typre
   use def_parame
   use Mod_Memor
   use Mod_Mesh
   use Mod_Element
   use Mod_CutMesh
   use Mod_LevelSet
   implicit none
   class(LevelSetProblem), target :: a

   interface
      subroutine lev_LayerOut(a, NelemWithAtLeastOneauxPointDone,ElementDone,PointDone,ElementListByLayers)
         use typre
         use Mod_LevelSet
         use Mod_Element
         implicit none
         class(LevelSetProblem) :: a
         integer(ip) :: NelemWithAtLeastOneauxPointDone
         logical :: ElementDone(*),PointDone(*)
         integer(ip) :: ElementListByLayers(*)
      end subroutine

   end interface

   class(FiniteElement), pointer :: e => NULL()

   integer(ip) :: icomp,ipoin,jpoin,npoin,ielem,nelem,nnode,inode,elemi,pelpo,helem,jelem
   integer(ip) :: elemh,hpoin,elemj

   !auxiliar variables
   integer(ip) :: auxlelemcut,auxnnode1,auxnnode2,kelem,existin,existout
   logical, allocatable     :: ElementDone(:),PointDone(:)
   !surface subroutine
   integer(ip) :: ndime,auxAdd,fnode,fpoin,gpoin,gnode
   real(rp)    :: dsurf,hclen,rdinv,cutcriteria
   real(rp), allocatable :: ellev(:,:)

   integer(ip) :: NelemWithAtLeastOneauxPointDone, ElementStatus
   real(rp), pointer :: level(:)=> Null()

   level => a%level(:,1)

   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNelem(nelem)
   call a%Mesh%GetNdime(ndime)

   rdinv=1.0_rp/real(ndime)

   !Definition of arrays needed by the level set

   !allocate
   call a%Mesh%ElementAlloc(e,a%Memor,'DefaultRule','lev_CutElementsAndListByLayers')
   call a%Memor%alloc(e%mnode,1,ellev,'ellev','lev_CutElementsAndListByLayers')
   call a%Memor%alloc(nelem,ElementDone,'ElementDone','lev_CutElementsAndListByLayers')
   call a%Memor%alloc(npoin,PointDone,'PointDone','lev_CutElementsAndListByLayers')

   !Initialization
   ElementDone = .false.
   PointDone=.false.

   call a%CutMesh%ComputeIntersectionPoints(level)

   !Fluid definition
   !List of point cuts by the intreface and point status to define the sign of then
   call a%CutMesh%ComputeElementType(level)
   call a%CutMesh%ComputePoinType(level)
   !Only testing
   call a%CutMesh%ComputeSubelements(ndime)
   call a%CutMesh%ComputeEnrichBubble(ndime)

   kelem=0_ip

   !to check if the element has two fluids
   ielem=1
   do ielem=1,nelem
      !Load Element
      call a%Mesh%ElementLoad(ielem,e)

      call a%CutMesh%GetElementType(ielem,ElementStatus)

      if (ElementStatus == 0.0_rp) then
         !obtaining surface and intersection points
         dsurf=0.0_rp
         call a%CutMesh%GetSurfaceIntersection(ielem,e,dsurf)

         !Derivative (and detjm) at the center of gravity
         call e%elmdcg

         !Element length
         hclen=(e%weicg*e%detjm)**rdinv
         hclen=hclen/100000_rp

         if(dsurf>0.0_rp)then

            !List of cut elements initially not full
            kelem=kelem + 1

            a%ElementListByLayers(kelem) = ielem
            ElementDone(ielem) = .true.

            !Load Element
            call a%Mesh%ElementLoad(ielem,e)

            do gnode=1,e%pnode
               gpoin=e%lnods(gnode)
               PointDone(gpoin) = .true.
            end do
          end if

      end if
   end do

   NelemWithAtLeastOneauxPointDone = kelem

   call lev_LayerOut(a,NelemWithAtLeastOneauxPointDone,ElementDone,PointDone,a%ElementListByLayers)

   !deallocate
   call a%Memor%dealloc(e%mnode,1,ellev,'ellev','lev_CutElementsAndListByLayers')
   call a%Memor%dealloc(nelem,ElementDone,'ElementDone','lev_CutElementsAndListByLayers')

   call a%Memor%dealloc(npoin,PointDone,'PointDone','lev_CutElementsAndListByLayers')
   call a%Mesh%ElementDeAlloc(e,a%Memor,'DefaultRule','lev_CutElementsAndListByLayers')

end subroutine

subroutine lev_LayerOut(a, NelemWithAtLeastOneauxPointDone,ElementDone,PointDone,ElementListByLayers)
   use typre
   use Mod_LevelSet
   use Mod_Element
   use Mod_HangingNodes
   implicit none
   class(LevelSetProblem) :: a
   integer(ip) :: NelemWithAtLeastOneauxPointDone
   logical :: ElementDone(*),PointDone(*)
   integer(ip) :: ElementListByLayers(*)

   class(FiniteElement), pointer :: e => NULL()
   logical, allocatable :: auxPointDone(:)
   real(rp), allocatable :: ellev(:,:)
   integer(ip), pointer     :: lelpo(:) => NULL()

   integer(ip) :: icomp,ipoin,jpoin,npoin,ielem,nelem,nnode,inode,elemi,pelpo,helem,jelem
   integer(ip) :: elemh,hpoin,elemj

   integer(ip) :: ndime,auxAdd,fnode,fpoin,gpoin,gnode,npoinInterface,kelem
   integer(ip), allocatable :: auxelementlist(:)
   integer(ip) :: kfl_Hanging

   interface
      subroutine CheckLogicalParents(a,ipoin,array,TrueOrFalse)
      use typre
      use Mod_Mesh
      implicit none
      class(FemMesh) :: a
      integer(ip) :: ipoin
      logical :: array(*)
      logical :: TrueOrFalse
      end subroutine

   end interface

   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNelem(nelem)
   call a%Mesh%GetNdime(ndime)

   allocate(auxelementlist(nelem))
   auxelementlist= 0

   call a%Mesh%ElementAlloc(e,a%Memor,'DefaultRule','lev_CutElementsAndListByLayers')
   call a%Memor%alloc(e%mnode,1,ellev,'ellev','lev_CutElementsAndListByLayers')

   call a%Memor%alloc(npoin,auxPointDone,'auxPointDone','lev_CutElementsAndListByLayers')
   call a%Mesh%GetHanging(kfl_Hanging)

   !To complete the list of element layers
   !Do not do it if no element is cut
   if (NelemWithAtLeastOneauxPointDone > 0) then
      kelem = NelemWithAtLeastOneauxPointDone
      auxPointDone = .false.
      do elemh=1,nelem
         helem = ElementListByLayers(elemh)
         if (helem < 1 .or. helem > nelem) then
            return
            !We will complete it in the following passes

!             do ielem = 1,nelem
!                if (auxelementlist(ielem) == 0) then
!                   write(*,*) 'asdfasdf'
!                endif
!             enddo
!             call runend('Error in lev_layerOut')
         endif
         auxelementlist(helem) = 1
         !Load Element
         call a%Mesh%ElementLoad(helem,e)

         do inode = 1, e%pnode
            hpoin=e%lnods(inode)

            if(auxPointDone(hpoin) .eqv. .false.)then
               auxPointDone(hpoin)=.true.

               call a%Mesh%GetLelpo(hpoin,pelpo,lelpo)

               do elemj=1,pelpo
               jelem = lelpo(elemj)

                  auxAdd=0
                  call a%Mesh%ElementLoad(jelem,e)
                  do fnode =1,e%pnode
                     fpoin =e%lnods(fnode)

                     !Check for hanging nodes
                     if (kfl_Hanging == 1) then
                        if (PointDone(fpoin) .eqv. .false.) then
                           call CheckLogicalParents(a%Mesh,fpoin,PointDone,PointDone(fpoin))
                        endif
                     endif

                     if(PointDone(fpoin) .eqv. .true.)then
                        auxAdd = auxAdd +1
                     end if
                  end do

                  ! in linear elements we add elements when only one node is not calculated
                  if((ElementDone(jelem) .eqv. .false.) .and. (auxAdd>=ndime))then

                     ElementDone(jelem) = .true.
                     kelem=kelem+1
                     ElementListByLayers(kelem)=jelem

                     do gnode=1,e%pnode
                        gpoin=e%lnods(gnode)
                        PointDone(gpoin) = .true.
                     end do

                  end if
               end do
            end if
         end do
      end do
   endif

   call a%Memor%dealloc(npoin,auxPointDone,'auxPointDone','lev_CutElementsAndListByLayers')
   call a%Memor%dealloc(e%mnode,1,ellev,'ellev','lev_CutElementsAndListByLayers')

   call a%Mesh%ElementDeAlloc(e,a%Memor,'DefaultRule','lev_CutElementsAndListByLayers')

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      SUBROUTINE PSUBV(A,LDA,NA,V,IND,NI,I)
      IMPLICIT NONE

*
*  -- user defined driver routine --
*
*
*     .. Scalar Arguments ..
      INTEGER            NA, NI, I, LDA
*     ..
*     .. Array Arguments ..
      INTEGER            IND(*)
      DOUBLE PRECISION   A(LDA,*), V(*)
*     ..
*
*  Purpose
*  =======
*
*  PSUBV computes the I-th column vector V of A. The elements of V
*  are reduced to the actual index set IND.
*
*
*  Arguments
*  =========
*
*  A        (input) DOUBLE PRECISION array of DIMENSION (NA,NA).
*           Before entry, the array A must contain matrix of
*           covariances.
*           Unchanged on exit.
*
*  NA       (input) INTEGER.
*           On entry, NA specifies the number of rows/columns of A
*           as declared in the calling (sub) program. NA >= 0.
*           Unchanged on exit.
*
*  V        (output) DOUBLE PRECISION vector of DIMENSION (NI).
*           On exit, the I-th column vector of A reduced to the
*           actual index set IND.
*
*  IND      (input) INTEGER vector, DIMENSION (NI).
*           On entry, IND specifies the indices forced into
*           the actual solution.
*           Unchanged on exit.
*
*  NI       (input) INTEGER.
*           On entry, NI specifies the dimension of IND. 0 <= NI <= NA.
*           Unchanged on exit.
*
*
*  =====================================================================
*
*     ..
*     .. Parameters ..
*     ..
*     .. Local Scalars ..
      INTEGER            J
*     ..
*     .. Local Arrays ..
*     ..
*     .. External Functions ..
*     ..
*     .. External Subroutines ..
*     ..
*     .. Intrinsic Functions ..
*     ..
*     .. Executable Statements ..
*
*
      DO 10 J=1,NI
         V(J)=0
         V(J)=A(IND(J),I)
  10   CONTINUE
      END

max_stars_repo_name: cran/edesign

max_issues_repo_name: cran/edesign

max_forks_repo_name: cran/edesign


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program barot_obs_file

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

! Generates a regular set, currently hardcoded to 40 by 30,
! of observation file locations for the forced barot model.

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

write(*, *) 1200
do i = 1, 40
   do j = 1, 30
      lon = ((i - 0.5) / 40.0) * 360.0
! Don't want to have to interpolate lats
      lat = -86.0 + (j / 30.0) * 170.0
      write(*, *) lon, lat, (1e6)**2
   end do
end do
end program barot_obs_file

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess