text:
One of Davis newer Bars, the G Street WunderBar is just one more addition to the Davis Night Time Nightlife.

The former G Street Pub has now reopened as G Street WunderBar.  It is owned by the former head bartender of G Street Pub Chris Armanini.  It opened on October 26, 2011.  Before it opened, the location underwent intensive remodeling including new bathrooms, a new commercial kitchen and an operable sewer system.  They now offer extensive food choices.  They  have four pool tables, vintage video games, dart boards, a full bar in the pool room, and free blazingfast WiFi.

They enforce some much stricter policies than the old G Street Pub, including no smoking outside the patio entrance of the building, no skateboards inside, and no messenger bags, large totes or backpacks inside after 9:30PM (people are more than welcome to study during the day and take full advantage of their free WiFi).  The backpack policy can vary depending on the nightly entertainment, when the kitchen closes, and when the 2nd Security Shift phases on.

On Mondays, the bar hosts a great Trivia Nights Trivia Night from 911 PM

As always, remember all of the Bar Taboos and etiquette regarding Tipping.

Happy Hour MF 4:00pm7:00pm
$1 off all draft
$1 PBR
$3 Well Drinks
$4 Pitchers of PBR
$4.95 Pound Basket of Wings

In its history, the location has transitioned from Davis Meat Market & Grocery to a very fancy restaurant (A. J. Bumps), to a biker bars biker bar, to the wellworn, allpurpose drinking spot that it now is.

 WunderBar in process

20111026 15:56:17 nbsp Wunderbar seems an appropriate name for this establishment, since racial fascism akin to Nazi Germany was evident in my experience there last night. My friends and I (a colorful and multiethnic group of seven) had just begun a relatively tame boys night by playing eight ball and drinking cheap beer, in a throwback to the G Street Pub days. During the second round, a bouncer approached the darkest member of our group and informed him that he was intoxicated and needed to leave the property or else the police would become involved. I interjected myself into the conversation to find out exactly what my friend had done wrong (at this point I was told I may have to leave as well since I was wearing a wool beanie and only hats with brims were allowed). The bouncer was only following orders and so we soon took up the discussion with the management. Before long the police showed up. One of the officers seemed in cahoots with the management, enjoying an impromptu tour of the new space and entertaining the suggestion by a manager that my friend be charged with a Drunk in Public offense. The other officer, slightly of a darker complexion and with more syllables in his name, appeared to understand the situation and apologized that this was a private establishment and that they had the right to order anyone off the property if they so desired. We agreed (happily) to leave. The conversation took place outside the bar and when we decided to call it quits we were denied reentry to claim our sweaters and other belongings. My friend was told that he was not allowed in the bar from here on out, with threat of arrest if seen on the premises again. This was at the suggestion from the fraternalminded and lighteyed officer. To my surprise, I was given the same sentence: entry never again or else arrest. This was my punishment for asking my friends crime and for explaining things to the officers. Needless to say, my friends and I (including our money, which in part pays the salary of all employees there) will gladly stay away. If you choose to give this bromantic getaway your business be careful, especially if you are not white or dont belong to the network. Users/Andresito

  You are, of course, referring to your friend, Sellou, who has been 86d from almost every bar in Davis, including easygoing Sophias by Kevin himself(!), for getting belligerently drunk and pulling the Race Card when any staff member refuses to serve him further or let him regain entrance to an establishment.  Good riddance.  Your friend has a HUGE chip on his entitlementminded, elitist shoulder and should give up drinking in Public, because he will be arrested if he decides to step foot in this bar again.  Users/HerrWunderbar

  Knowing each of the staff members there, with the exception of two door personnel, I can personally vouch for everyone there that there is absolutely no racism in the establishment. If there was a problem with any member of your group I am sure it had absolutely nothing to do with race and everything to do with the level of intoxication of the group/individuals.  Users/WesP

20111028 09:41:17 nbsp Im willing to give the bar the benefit of the doubt. As for the police, I had a friend who was with the department for a few months and he said it was the most racist place hed ever worked (racial profiling, jokes, etc.)  Users/ScottSchiller
  I hate to say it, but that names just about every bar in this town, at least as far as profiling and jokes go. Ask any security people who causes the most problems and many wil say it is the townies, ask them who uses the most fake IDs and it will be women, usually of Asian or Middle Eastern descent... ask who tips the least and it will be Southeast Asians and Northern Europeans. Its not so much profiling as it is trends shared by specific groups of individuals. As for fights, it is almost always either drunk college jocks (regardless of color or ethnicity) or it is Asian/African American individuals. Its just what I have seen since working in the bar scene. I cant speak for the police, but as far as the bars go, it is just the trend of how things usually go in this town. I am quite sure this trend is different in every other city and even from bar to bar.  Users/WesP
   Ask whos ridiculously douchey and arrogant toward staff and its law students. Users/TomGarberson

20111028 15:46:45 nbsp My experience is that the food/drink service industry itself is largely racist.  Ask anyone whos ever worked waitstaff and they will tell you stories.   Users/OliviaY

     Nothing says nuance like a blanket statement backed by allegorical evidence. Users/HerrWunderbar

20111028 22:20:42 nbsp So this is where all of the locals hang out... Granted, it is one of the first dive bars Ive gone to. They have a nobeanies policy, they serve jaeger on tap, and there are automatic hand sanitizer machines in case you feel extra dirty. Users/argyle

20111108 14:37:50 nbsp anone know if there still do karoake tuesday nights?  Users/RobertSinclair

      Yes, from 9:30 til Close on Tuesdays.  Users/HerrWunderbar

20111123 16:29:29 nbsp They dont allow messenger bags or backpacks nor do they allow beanies. Can anyone explain to me the reasons for these policies and why no other bar in Davis enforces them? Also, like the old G st, they still have a $15 min for credit cards which is 50% more than any other bar in Davis. Users/JakeJames
   Bags take up space/trip people/could carry in outside booze?  Beanies: no hipsters allowed.

20111208 09:05:38 nbsp Since when did G St raise their standards? I remember the days hobos were allowed into the bar. Now they reject people simply for stumbling over. No big deal. Except now they get a negative comment on their daviswiki page. You suck new Wunderbar!!!!

PS I agree everyone that works there is a racist. But I could care less about that part. Users/msavrass

      No bar in town has so lenient of a door policy that hobos and stumbledrunks are allowed in. Why do you expect WunderBar to be any different, especially since its in no way affiliated to the old Pub?  Users/HerrWunderbar

20111217 17:23:29 nbsp The lights are way too bright for a bar or a restaurant for that fact.The bar/eaterys ambiance seems very sterile,kind of reminds me of a hospital waiting room especially with the bright florescent lights.Im kind of disappointed with the newGStreet barand did i mention the lights are too bright??
   Users/HallandOats

20120131 08:46:21 nbsp We visited this place for the first time on 01.23.12. The place is a dive. Not to mention a cowboy dress code only. The bar has more lights than a stadium. Youll need your sunglasses for this dive. While we sat at an unoccupied table a cowboy came up to us to tells us to vacate the table because he and his dyed blonde were there first. However, no one was at the table. I and my friend told him to take a hike and that made the little cowboy pissed. He threaten us with harm but once the shithead understood that we had attorneys (relatives) that could rip him a new AHole. He backed off. This place has a NeoNazi environment mixed in with the wanna be cowboys. Dont waste your time patronizing this dive. Spend your hard earned money elsewhere. Lets show this shithole of a place that we dont want it in Davis,Ca. Our town deserves better.  Vic .L  Users/VicttorLopez

    So you get into a baby bitchfest with another customer and its the bars fault (even though the bar backed up your side of the story)?  Im sure your incredibly great sense of judgement will be sorely missed there, especially with your drunk asses being thrown to the ground at Froggys later that night.  How was the ride in the Davis PD patrol car, comfortable? It sounds like you get the full use of our Citys wonderful services for which the rest of us have to pay.  Youre right, our town does deserve better: stay home.   Users/HerrWunderbar

20120312 12:13:10 nbsp Just went here for the first time since it reopened! Reminded of the old undergrad days. Pitchers of beer, awesome games of darts. You guys did a great job, glad to see this place back in action! Users/ShahinVahdat

20120619 00:36:03 nbsp I enjoy the hokcey game. Greatly. Still wishing for lower credit card minimums. Sometimes you just want a pint, not three or four. I know, I know, cash. Sigh. First world problems.  Users/JakeJames

20120709 15:05:34 nbsp The owner Chris was rude as a bartender (as were most behind the bar at G Street) so Im not patronizing the place now thats hes an owner. Waste of bandwidth to even elaborate. Users/MaryCarroll

     You are clearly a very patient and understanding customer to have sat through numerous interactions, on multiple occasions, to have formed such an opinion about the majority of the frontofhouse staff at the old venue; therefore, if you could just rechannel that sympathetic energy, were sure that you would find the new bar to be leaps and bounds above the old one in terms of ambiance, product, and service.  Please, come in and introduce yourself to me and see the improvements for yourself, and, as a token, Ill be happy to buy your first libation.  Users/HerrWunderbar

20120817 13:46:41 nbsp wunderbar is alright. Yes they still employ almost all of the racist/aggresive/nonintelligent aholes that swear they arent racist but you really have to be a drunken fool to think you will ever get your way at a bar vs bars owners/staff. Especially in a college town where it is assumed all kids arguing are just belligerant.

Now like everything in this world there are exceptions. For one, bartenders drew and another skinny, shorthaired guy, lighter skinned (looks like he just turned 21, lol) are really aware of whos tipping and whos not. Naturally, they will provide me with outstanding service because i tip them well. Secondly, i really like that our shoes dont stick to the floor and that girls we go out with arent twisting their ankles because of it (lol). Bathrooms were a MAJOR improvement.

What would make this an even better establishment is a better background check on bouncers. This is where the most complaints of racism/mistreatment comes from. If more than a person cries wolf, im going to believe theres a wolf. If anything maybe providing a more proactive approach to customer service training/awareness (because we are indeed customers and not animals to be treated the way bouncers believe they have the authority to do so)would cancel the need for extensive background checks (i hope you know this excludes a psychological evaluation which should be mandatory for anyone manhandling another humanbeing) and i bet all the allegations against this bar would dissipate in a matter of time.

Good work so far, but it could be even better if you stepped outside of your comfort zone and didnt hire from the more than circulated and neverrenewing bouncer pool of davis. Users/Ponchelas

20120928 11:16:52 nbsp Hmm, lets see:
racial fascism akin to Nazi Germany
the most racist place hed ever worked
everyone that works there is a racist
The place is a dive to which the comeback response included a baby bitchfest by the Wunderbars bitch and
the call for better background checks on the bouncers.

Well I became curious (before reading all of those comments) so I showed up.  Man, I wish I had read those comments before stepping into this place.

But, to be fair, you should see for yourself.  You too can be manhandled by the pathetic group of mindless muscles who call themselves bouncers. Users/DavidBarnum

20130420 20:01:37 nbsp Wow, Im actually kinda shocked at these comments. Ive been going to this bar for almost ten years. Ive never had any issues with racism. Sure there are a few assholes on staff, but thats anywhere you go. It is a dive bar. Its not that bad. Granted I usually go on Saturdays when they were having metal nights, Im friends with Mortal Atrocity and Wolfsblood members. As well as on occasion in the early evening to have a few drinks and play pool or foozhockey thing in the corner with friends. Never have had any trouble with the staff, only shitheads Ive had issues with have been the college kids. Stupid frat guys.

As far as the beanies and dress code shindig, maybe it is because Im a local and been going for so long I have never gotten any grief over the way Ive dressed. But nor have I seen anyone get harassed for it besides one time. Swear this guy had his pants to his shins. Never understood that style.

A few things I do like better with wunderbar. Better food, and bathrooms are way cleaner. Everything else pretty much feels the same to me. Its still my favorite bar in davis.  Users/Dwajon

20130424 17:35:30 nbsp My favorite bar in Davis, period. Best beer selection, excellent food for a bar (the burgers are magical!), pool tables that work, and a staff that is always warm, friendly and attentive. I really cant express in words how glad I am it reopened. I missed it. Now that its back, every other bar in Davis is difficult to compare. Users/BrianPakpour

20130526 15:39:55 nbsp Hello, I am 44 year old Chilean woman that does not have the stereotypical latin look, with a family, work, and a pretty normal Davis life. I hardly ever go out to Bars in Davis eventhough I have lived in Davis for 12 years. I feel the need to explain the following because I read the other comments on Daviswiki. To my surprise they have mainly to do with racisiim and above all the responses seem ridiculous and somewhat sarcastic. In my opinion someone that does not have the guts to use there name is not someone that should be answering the comments. But it is true that there is racisim in this place because last night, May 25, 2013, I went out with two girl friends to distract ourselves from life a while and by pure chance we ended up at the G Street Pub. We arrived around 9:00 pm. A bit after we arrived a man approached us with a red shirt and sat at our table and asked us from where we are from. My friend that speaks perfect english ( unlike myself that unfortunately does not speak it so well, but I do understand quite a lot ) says I am American but my family is Japanese, she is Chilean, and she is German. Then the man looks at me and says, OH, then I am from Guatemala and I responded great then you will be able to understand me, and I began to speak with him in Spanish and he looks at me and says, Oh I am just joking and I responded that I was not up for jokes and for him to leave my table. This is when I realized that the people with the red shirts are employees of the bar, possibly the door men or bouncers. It called attention, until this point this does not matter much but it does make sense with later events, because it was the most humiliating, embarrassing, and discriminating thing I have ever lived in my entire life.
When I arrived I ordered a Tequila Screwdriver with Don Julio, which they serve in a short glass. At around 11:30 pm we moved to the pool table area to be able to talk because the music was to loud. This is when I ordered my second Tequila Screwdriver with Don Julio of the night ( 2 hours later ). At 1:00 am my one friend wanted to leave so my other friend offered to take her home and come back because I did not want to leave just yet because I was having a good time with some other patrons of the bar that were teaching me to play pool. I was having lots of beginners luck and made the ball in the hole 2 times in a row which made me very happy and excited. I went to ask for my third drink of the night and asked for the same as before but the bar man left to poor the drink before I could say that it was with Don Julio, but since I had ordered two previous drinks like that I figured he left quickly because he knew what I was drinking since like any good barman. When he returned I asked politely if the drink contained Don Julio and he said, No you should have specified before I served it. I said I was sorry,( even though he was the one that did not let me finish my sentence before walking away ) and asked if he could change it for Don Julio because cheap tequila gives me a head ache, and he responded, No, you have to pay for this drink, so I again said I am sorry for the mistake and said that I no longer wanted anything. I went to sit on a chair to wait for my friend to pick me up. A few seconds later to men came up to me and asked me to leave, I asked why thinking it might be closing time. They told me it was because I did not want ot pay for a drink. I tried to explain with my little English what had happened and I noticed that they did not understand so I asked to barrow a phone. At first they said no, that I must leave. I said I know I have to leave but I need a phone to call my friend or husband to come pick me up, I expressed that I did not want to be alone standing outside the bar. Finally I was given a phone to use, my friend did not answer so I called my husband. I asked if I could wait there until he arrives and was told No, and a kind patron of the bar accompanied me outside. When my husband arrived in about 8 minutes because we live close by, I was not allowed inside to help clear up the misunderstanding and he was basically told a bunch of lies. This is the important part for me, I was told to leave because I did not pay for a drink. When my husband went inside to understand what happened he was told I was asked to leave for STEALING and besides this reason, that I had had to much to drink. And the most incredible part is that the manager told him that he had spoken to me in spanish and that he would never discriminate his people, and I believe him, he is discriminating me for not being his people. Yeh Right like they say her in the US. Lets use some logic, if he would have spoken to me in Spanish would we have had the misunderstanding over the tequila. And if you really do speak spanish, what a shame because you do not know how to use it. Please!!!, you are racist and a liar because all you said about me is not true, you have zero etiquette for your patrons. I understand and respect that we all have the right to ask someone to leave our private properties, but just be truthful, you guys ask people to leave because of a problem with RACISM. So if this is not racism what is it? You do not need to discredit citizens of Davis like myself, saying that I am a thief, or to drunk, or did not order the drink appropriately, or that I behaved inappropriately to hide your true reason.

P.S. For he that answers ( which who knows if it is the owner ), since your so cool and have an answer for everything, you can put your answer where it fits you best. It is I that will never go to your establishment again. It is really sad to me that this bar dirties the beauty and tranquility of this city because the people of Davis, since I have arrived have been the nicest people I have meet since being outside my country. In my opinion this Antro is a blemish on Davis.

  Users/ximemurillo

20130729 17:08:25 nbsp Love all the bartenders, have never had a negative interaction with any of them, nor the security guards in the past. Was there last week though, and a guy in a red shirt that Id not noticed before was running his mouth more like a patron than an employee. He was going on about how someone knew all the words to the Backdoor Boys and making negative gayshaming comments about Lance Bass. Its one thing if youre there drinking and want to make homophobic comments, but when you are representing G St. as an employee, should you really be making bad jokes that could / did offend customers? Users/OhHey

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module mg_intergrids

  use mg_mpi
  use mg_tictoc
  use mg_namelist
  use mg_grids
  use mg_mpi_exchange
  use mg_gather

  implicit none

contains
  !---------------!
  !- FINE2COARSE -! fine to coarse grid
  !------------------------------------------------------------
  subroutine fine2coarse(lev)

    ! coarsen grid(lev)%r to grid(lev+1)%b and set grid(lev+1)%p=0

    integer(kind=ip), intent(in) :: lev

    real(kind=rp),dimension(:,:,:),pointer :: r
    real(kind=rp),dimension(:,:,:),pointer :: b

    integer(kind=ip) :: nx, ny, nz

    nx = grid(lev+1)%nx
    ny = grid(lev+1)%ny
    nz = grid(lev+1)%nz

    r => grid(lev)%r

    if (grid(lev+1)%gather == 1) then
       !- b here is sub level rhs, we will gather at the end
       !- in the level+1 rhs
       b => grid(lev+1)%dummy3
       !- ngx and ngy are the number of processes being gathered in x and y
       !- values are 1 or 2
       nx = grid(lev+1)%nx / grid(lev+1)%ngx
       ny = grid(lev+1)%ny / grid(lev+1)%ngy
!       if(myrank == 0) write(*,*)"gather lev=",lev+1,"nx,ny,nz=",nx,ny,nz
    else
       b => grid(lev+1)%b
!       if(myrank == 0) write(*,*)"F2C   lev=",lev+1,"nx,ny,nz=",nx,ny,nz
    endif

    if ((aggressive).and.(lev==1)) then
       call fine2coarse_aggressive(r,b,nx,ny,nz)

    elseif (grid(lev)%nz == 1) then
       call fine2coarse_2D(r,b,nx,ny)

    else

       call fine2coarse_3D(r,b,nx,ny,nz)

    end if

    if (grid(lev+1)%gather == 1) then
!       if(myrank == 0) write(*,*)" *** dummy3(1,1,1)=",b(1,1:ny,1:nx)
       r => grid(lev+1)%dummy3
       b => grid(lev+1)%b
       call gather(lev+1,r,b)
!       if(myrank == 0) write(*,*)" *** after gather **** b(1,1,1)=",grid(lev+1)%b(1,1,1)
    endif

    call fill_halo(lev+1,grid(lev+1)%b)

    grid(lev+1)%p = 0._rp

  end subroutine fine2coarse

  !----------------------------------------
  subroutine fine2coarse_aggressive(x,y,nx,ny,nz)

    real(kind=rp)   , dimension(:,:,:), intent(in)    :: x !fine
    real(kind=rp)   , dimension(:,:,:), intent(inout) :: y ! coarse
    integer(kind=ip)                  , intent(in)    :: nx, ny, nz

    ! local
    integer(kind=ip):: i,j,k,k2

    do k=1,nz
       k2=(k-1)/8+1
       if(mod(k,8).eq.1)then
          do j=1,ny
             do i=1,nx
                y(i,j,k2) = x(i,j,k)*0.125_rp
             enddo
          enddo
       else
          do j=1,ny
             do i=1,nx
                y(i,j,k2) = y(i,j,k2)+x(i,j,k)*0.125_rp
             enddo
          enddo
       endif
    enddo

  end subroutine fine2coarse_aggressive

  !------------------------------------------------------------
  subroutine fine2coarse_2D(x,y,nx,ny)

    real(kind=rp),dimension(:,:,:),pointer,intent(in) :: x
    real(kind=rp),dimension(:,:,:),pointer,intent(out) :: y
    integer(kind=ip), intent(in) :: nx, ny

    integer(kind=ip) :: i,j,i2,j2
    integer(kind=ip) :: d

    d = size(x,1) ! vertical dimension of the fine level, can 2 or 1

   if(d==1)then
       ! x was already 2D
       do i2=1,nx
          i=2*i2-1
          do j2=1,ny
             j=2*j2-1
             y(1,j2,i2) = (x(1,j,i)+x(1,j,i+1)+x(1,j+1,i)+x(1,j+1,i+1))
          enddo
       enddo
    else
       ! x was 3D
       do i2=1,nx
          i=2*i2-1
          do j2=1,ny
             j=2*j2-1
             y(1,j2,i2) = (x(1,j,i)+x(1,j,i+1)+x(1,j+1,i)+x(1,j+1,i+1)&
                          +x(2,j,i)+x(2,j,i+1)+x(2,j+1,i)+x(2,j+1,i+1))
          enddo
       enddo
    endif

  end subroutine fine2coarse_2D

  !------------------------------------------------------------
  subroutine fine2coarse_3D(x,y,nx,ny,nz)
    !
    ! Fine2coarse 'x' from fine level l1 to 'y' on coarse level l2=l1+1
    real(kind=rp),dimension(:,:,:),pointer,intent(in) :: x
    real(kind=rp),dimension(:,:,:),pointer,intent(out) :: y
    integer(kind=ip), intent(in) :: nx, ny, nz
    ! local
    integer(kind=ip) :: i,j,k,i2,j2,k2
    real(kind=rp):: z

    do i2=1,nx
       i=2*i2-1
       do j2=1,ny
          j=2*j2-1
          do k2=1,nz
             k=2*k2-1
             z = x(k,j,i)  +x(k,j,i+1)  +x(k,j+1,i)  +x(k,j+1,i+1) &
                  + x(k+1,j,i)+x(k+1,j,i+1)+x(k+1,j+1,i)+x(k+1,j+1,i+1)
             y(k2,j2,i2) = z
          enddo
       enddo
    enddo

  end subroutine fine2coarse_3D

  !---------------!
  !- COARSE2FINE -! coarse to fine grid
  !------------------------------------------------------------
  subroutine coarse2fine(lev)

    ! interpolate grid(lev+1)%p to grid(lev)%r and add it to grid(lev)%p

    integer(kind=ip), intent(in) :: lev

    real(kind=rp),dimension(:,:,:),pointer :: rf
    real(kind=rp),dimension(:,:,:),pointer :: pc

    integer(kind=ip) :: nxc, nyc, nzc

    nxc = grid(lev+1)%nx
    nyc = grid(lev+1)%ny
    nzc = grid(lev+1)%nz

    pc => grid(lev+1)%p

    if (grid(lev+1)%gather == 1) then
       !       if(myrank==0)write(*,*)"SPLIT!!!"
       rf => grid(lev+1)%dummy3
       call split(lev+1,pc,rf)
       pc => grid(lev+1)%dummy3
       nxc = grid(lev+1)%nx / grid(lev+1)%ngx
       nyc = grid(lev+1)%ny / grid(lev+1)%ngy
       nzc = grid(lev+1)%nz
    endif

    rf => grid(lev)%r

    if (trim(interp_type)=='nearest')  then

       if ((aggressive).and.(lev==1)) then

          call coarse2fine_aggressive(rf,pc,nxc,nyc,nzc)

       elseif (grid(lev)%nz == 1) then
          call coarse2fine_2D_nearest(rf,pc,nxc,nyc)
       else
          call coarse2fine_3D_nearest(rf,pc,nxc,nyc,nzc)
       end if

    elseif ( trim(interp_type)=='linear')then

       if ((aggressive).and.(lev==1)) then

          call coarse2fine_aggressive(rf,pc,nxc,nyc,nzc)

       elseif (grid(lev)%nz == 1) then

          call coarse2fine_2D_linear(rf,pc,nxc,nyc)
       else
          call coarse2fine_3D_linear(rf,pc,nxc,nyc,nzc)

       end if

    endif

    call fill_halo(lev,grid(lev)%r)

    grid(lev)%p = grid(lev)%p + grid(lev)%r

  end subroutine coarse2fine

  !------------------------------------------------------------
  subroutine coarse2fine_aggressive(x,y,nx,ny,nz)

    real(kind=rp),dimension(:,:,:),pointer,intent(in)  :: x
    real(kind=rp),dimension(:,:,:),pointer,intent(out) :: y
    integer(kind=ip),intent(in) :: nx, ny, nz

    !TODO
    integer(kind=ip) ::idum ! line to remove
    idum = nx               ! line to remove
    idum = ny               ! line to remove
    idum = nz               ! line to remove
    y = x                   ! line to remove
    write(*,*)'Error:  coarse2fine_aggressive not available yet !'
    stop -1
    !TODO

  end subroutine coarse2fine_aggressive

  !------------------------------------------------------------
  subroutine coarse2fine_2D_nearest(xf,xc,nx,ny)

    real(kind=rp),dimension(:,:,:),pointer,intent(out):: xf
    real(kind=rp),dimension(:,:,:),pointer,intent(in) :: xc
    integer(kind=ip),intent(in) :: nx, ny

    ! local
    integer(kind=ip) :: i,j,i2,j2
    integer(kind=ip) :: d

    d = size(xf,1) ! vertical dimension of the fine level, can 2 or 1

    if(d==1)then
       ! xf is also 2D
       do i2=1,nx
          i=2*i2-1
          do j2=1,ny
             j=2*j2-1
             xf(1,j  ,i  ) = xc(1,j2,i2)
             xf(1,j+1,i  ) = xc(1,j2,i2)
             xf(1,j  ,i+1) = xc(1,j2,i2)
             xf(1,j+1,i+1) = xc(1,j2,i2)
          enddo
       enddo
    else
       ! xf is 3D
       do i2=1,nx
          i=2*i2-1
          do j2=1,ny
             j=2*j2-1
             xf(1,j  ,i  ) = xc(1,j2,i2)
             xf(2,j  ,i  ) = xc(1,j2,i2)
             xf(1,j+1,i  ) = xc(1,j2,i2)
             xf(2,j+1,i  ) = xc(1,j2,i2)
             xf(1,j  ,i+1) = xc(1,j2,i2)
             xf(2,j  ,i+1) = xc(1,j2,i2)
             xf(1,j+1,i+1) = xc(1,j2,i2)
             xf(2,j+1,i+1) = xc(1,j2,i2)
          enddo
       enddo
    endif

  end subroutine coarse2fine_2D_nearest

  !------------------------------------------------------------
  subroutine coarse2fine_2D_linear(xf,xc,nx,ny)

    real(kind=rp),dimension(:,:,:),pointer,intent(out) :: xf
    real(kind=rp),dimension(:,:,:),pointer,intent(in)  :: xc
    integer(kind=ip),intent(in) :: nx, ny

    ! local
  integer(kind=ip) :: i,j,i2,j2,k,k2
  real(kind=rp) :: a,b,c
  !
  ! weights for bilinear in (i,j)
  a = 9._rp / 16._rp
  b = 3._rp / 16._rp
  c = 1._rp / 16._rp

  k  = 1
  k2 = 1

  do i2=1,nx
     i=2*i2-1
     do j2=1,ny
        j=2*j2-1
        xf(k  ,j  ,i  ) =  &
             + a * xc(k2,j2  ,i2) + c * xc(k2,j2-1,i2-1) &
             + b * xc(k2,j2-1,i2) + b * xc(k2,j2  ,i2-1)
        xf(k  ,j+1,i  ) =  &
             + a * xc(k2,j2  ,i2) + c * xc(k2,j2+1,i2-1) &
             + b * xc(k2,j2+1,i2) + b * xc(k2,j2  ,i2-1)
        xf(k  ,j  ,i+1) =  &
             + a * xc(k2,j2  ,i2) + c * xc(k2,j2-1,i2+1) &
             + b * xc(k2,j2-1,i2) + b * xc(k2,j2  ,i2+1)
        xf(k  ,j+1,i+1) =  &
             + a * xc(k2,j2  ,i2) + c * xc(k2,j2+1,i2+1) &
             + b * xc(k2,j2+1,i2) + b * xc(k2,j2  ,i2+1)
     enddo
  enddo

  end subroutine coarse2fine_2D_linear

  !------------------------------------------------------------
  subroutine coarse2fine_3D_nearest(xf,xc,nx,ny,nz)

    real(kind=rp),dimension(:,:,:),pointer,intent(out) :: xf
    real(kind=rp),dimension(:,:,:),pointer,intent(in)  :: xc
    integer(kind=ip),intent(in) :: nx, ny, nz

    ! local
    integer(kind=ip) :: i,j,k,i2,j2,k2
    !
    do i2=1,nx
       i=2*i2-1
       do j2=1,ny
          j=2*j2-1
          do k2=1,nz
             k=2*k2-1
             xf(k  ,j  ,i  ) = xc(k2,j2,i2)
             xf(k+1,j  ,i  ) = xc(k2,j2,i2)
             xf(k  ,j+1,i  ) = xc(k2,j2,i2)
             xf(k+1,j+1,i  ) = xc(k2,j2,i2)
             xf(k  ,j  ,i+1) = xc(k2,j2,i2)
             xf(k+1,j  ,i+1) = xc(k2,j2,i2)
             xf(k  ,j+1,i+1) = xc(k2,j2,i2)
             xf(k+1,j+1,i+1) = xc(k2,j2,i2)
          enddo
       enddo
    enddo

  end subroutine coarse2fine_3D_nearest

  !------------------------------------------------------------
  subroutine coarse2fine_3D_linear(xf,xc,nx,ny,nz)

    real(kind=rp),dimension(:,:,:),pointer,intent(out) :: xf
    real(kind=rp),dimension(:,:,:),pointer,intent(in)  :: xc
    integer(kind=ip),intent(in) :: nx, ny, nz

    ! local
    integer(kind=ip) :: i,j,k,i2,j2,k2,kp
    real(kind=rp) :: a,b,c,d,e,f,g
    !
    ! weights for bilinear in (i,j), nearest in k
    a = 9._rp / 16._rp
    b = 3._rp / 16._rp
    c = 1._rp / 16._rp
    !
    ! weights for trilinear in (i,j,k)
    d = 27._rp / 64._rp
    e =  9._rp / 64._rp
    f =  3._rp / 64._rp
    g =  1._rp / 64._rp
    !
    do i2=1,nx
       i=2*i2-1
       do j2=1,ny
          j=2*j2-1
          ! bottom level
          k  = 1
          k2 = 1
          xf(k  ,j  ,i  ) =  &
               + a * xc(k2,j2  ,i2) + c * xc(k2,j2-1,i2-1) &
               + b * xc(k2,j2-1,i2) + b * xc(k2,j2  ,i2-1)
          xf(k  ,j+1,i  ) =  &
               + a * xc(k2,j2  ,i2) + c * xc(k2,j2+1,i2-1) &
               + b * xc(k2,j2+1,i2) + b * xc(k2,j2  ,i2-1)
          xf(k  ,j  ,i+1) =  &
               + a * xc(k2,j2  ,i2) + c * xc(k2,j2-1,i2+1) &
               + b * xc(k2,j2-1,i2) + b * xc(k2,j2  ,i2+1)
          xf(k  ,j+1,i+1) =  &
               + a * xc(k2,j2  ,i2) + c * xc(k2,j2+1,i2+1) &
               + b * xc(k2,j2+1,i2) + b * xc(k2,j2  ,i2+1)
          ! interior level
          do k=2,nz*2-1
             ! kp = k2+1 for k=2,4 ..
             ! kp = k2-1 for k=3,5 ..
             k2 = ((k+1)/2)
             kp = k2-(mod(k,2)*2-1)
             xf(k  ,j  ,i  ) =  &
                  + d * xc(k2,j2  ,i2) + f * xc(k2,j2-1,i2-1) &
                  + e * xc(k2,j2-1,i2) + e * xc(k2,j2  ,i2-1) &
                  + e * xc(kp,j2  ,i2) + g * xc(kp,j2-1,i2-1) &
                  + f * xc(kp,j2-1,i2) + f * xc(kp,j2  ,i2-1)
             xf(k  ,j+1,i  ) =  &
                  + d * xc(k2,j2  ,i2) + f * xc(k2,j2+1,i2-1) &
                  + e * xc(k2,j2+1,i2) + e * xc(k2,j2  ,i2-1) &
                  + e * xc(kp,j2  ,i2) + g * xc(kp,j2+1,i2-1) &
                  + f * xc(kp,j2+1,i2) + f * xc(kp,j2  ,i2-1)
             xf(k  ,j  ,i+1) = &
                  + d * xc(k2,j2  ,i2) + f * xc(k2,j2-1,i2+1) &
                  + e * xc(k2,j2-1,i2) + e * xc(k2,j2  ,i2+1) &
                  + e * xc(kp,j2  ,i2) + g * xc(kp,j2-1,i2+1) &
                  + f * xc(kp,j2-1,i2) + f * xc(kp,j2  ,i2+1)
             xf(k  ,j+1,i+1) = &
                  + d * xc(k2,j2  ,i2) + f * xc(k2,j2+1,i2+1) &
                  + e * xc(k2,j2+1,i2) + e * xc(k2,j2  ,i2+1) &
                  + e * xc(kp,j2  ,i2) + g * xc(kp,j2+1,i2+1) &
                  + f * xc(kp,j2+1,i2) + f * xc(kp,j2  ,i2+1)
          enddo
          ! top level
          k = nz*2
          xf(k  ,j  ,i  ) =  0.5_rp                      * ( &
               a * xc(k2,j2  ,i2) + c * xc(k2,j2-1,i2-1) +   &
               b * xc(k2,j2-1,i2) + b * xc(k2,j2  ,i2-1)   )
          xf(k  ,j+1,i  ) =  0.5_rp                      * ( &
               a * xc(k2,j2  ,i2) + c * xc(k2,j2+1,i2-1) +   &
               b * xc(k2,j2+1,i2) + b * xc(k2,j2  ,i2-1)   )
          xf(k  ,j  ,i+1) =  0.5_rp                      * ( &
               a * xc(k2,j2  ,i2) + c * xc(k2,j2-1,i2+1) +   &
               b * xc(k2,j2-1,i2) + b * xc(k2,j2  ,i2+1)   )
          xf(k  ,j+1,i+1) =  0.5_rp                      * ( &
               a * xc(k2,j2  ,i2) + c * xc(k2,j2+1,i2+1) +   &
               b * xc(k2,j2+1,i2) + b * xc(k2,j2  ,i2+1)   )
       enddo
    enddo

  end subroutine coarse2fine_3D_linear

end module mg_intergrids

max_stars_repo_name: nmolem/mgroms

max_issues_repo_name: nmolem/mgroms

max_forks_repo_name: nmolem/mgroms


text:
MODULE Nelder_Mead

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
PRIVATE
PUBLIC :: minim

CONTAINS

SUBROUTINE minim(p, step, nop, func, maxfn, iprint, stopcr, nloop, iquad,  &
                 simp, var, functn, ifault)

!     A PROGRAM FOR FUNCTION MINIMIZATION USING THE SIMPLEX METHOD.

!     FOR DETAILS, SEE NELDER & MEAD, THE COMPUTER JOURNAL, JANUARY 1965

!     PROGRAMMED BY D.E.SHAW,
!     CSIRO, DIVISION OF MATHEMATICS & STATISTICS
!     P.O. BOX 218, LINDFIELD, N.S.W. 2070

!     WITH AMENDMENTS BY R.W.M.WEDDERBURN
!     ROTHAMSTED EXPERIMENTAL STATION
!     HARPENDEN, HERTFORDSHIRE, ENGLAND

!     Further amended by Alan Miller
!     CSIRO Division of Mathematical & Information Sciences
!     Private Bag 10, CLAYTON, VIC. 3169

!     Fortran 90 conversion by Alan Miller, June 1995
!     Alan.Miller @ vic.cmis.csiro.au
!     Latest revision - 5 December 1999

!     ARGUMENTS:-
!     P()     = INPUT, STARTING VALUES OF PARAMETERS
!               OUTPUT, FINAL VALUES OF PARAMETERS
!     STEP()  = INPUT, INITIAL STEP SIZES
!     NOP     = INPUT, NO. OF PARAMETERS, INCL. ANY TO BE HELD FIXED
!     FUNC    = OUTPUT, THE FUNCTION VALUE CORRESPONDING TO THE FINAL
!                 PARAMETER VALUES.
!     maxfn     = INPUT, THE MAXIMUM NO. OF FUNCTION EVALUATIONS ALLOWED.
!               Say, 20 times the number of parameters, NOP.
!     IPRINT  = INPUT, PRINT CONTROL PARAMETER
!                 < 0 NO PRINTING
!                 = 0 PRINTING OF PARAMETER VALUES AND THE FUNCTION
!                     VALUE AFTER INITIAL EVIDENCE OF CONVERGENCE.
!                 > 0 AS FOR IPRINT = 0 PLUS PROGRESS REPORTS AFTER EVERY
!                     IPRINT EVALUATIONS, PLUS PRINTING FOR THE INITIAL SIMPLEX.
!     STOPCR  = INPUT, STOPPING CRITERION.
!               The criterion is applied to the standard deviation of
!               the values of FUNC at the points of the simplex.
!     NLOOP   = INPUT, THE STOPPING RULE IS APPLIED AFTER EVERY NLOOP
!               FUNCTION EVALUATIONS.   Normally NLOOP should be slightly
!               greater than NOP, say NLOOP = 2*NOP.
!     IQUAD   = INPUT, = 1 IF FITTING OF A QUADRATIC SURFACE IS REQUIRED
!                      = 0 IF NOT
!               N.B. The fitting of a quadratic surface is strongly
!               recommended, provided that the fitted function is
!               continuous in the vicinity of the minimum.   It is often
!               a good indicator of whether a premature termination of
!               the search has occurred.
!     SIMP    = INPUT, CRITERION FOR EXPANDING THE SIMPLEX TO OVERCOME
!               ROUNDING ERRORS BEFORE FITTING THE QUADRATIC SURFACE.
!               The simplex is expanded so that the function values at
!               the points of the simplex exceed those at the supposed
!               minimum by at least an amount SIMP.
!     VAR()   = OUTPUT, CONTAINS THE DIAGONAL ELEMENTS OF THE INVERSE OF
!               THE INFORMATION MATRIX.
!     FUNCTN  = INPUT, NAME OF THE USER'S SUBROUTINE - ARGUMENTS (P,FUNC)
!               WHICH RETURNS THE FUNCTION VALUE FOR A GIVEN SET OF
!               PARAMETER VALUES IN ARRAY P.
!****     FUNCTN MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
!     IFAULT  = OUTPUT, = 0 FOR SUCCESSFUL TERMINATION
!                 = 1 IF MAXIMUM NO. OF FUNCTION EVALUATIONS EXCEEDED
!                 = 2 IF INFORMATION MATRIX IS NOT +VE SEMI-DEFINITE
!                 = 3 IF NOP < 1
!                 = 4 IF NLOOP < 1

!     N.B. P, STEP AND VAR (IF IQUAD = 1) MUST HAVE DIMENSION AT LEAST NOP
!          IN THE CALLING PROGRAM.

!*****************************************************************************

INTEGER, INTENT(IN)        :: nop, maxfn, iprint, nloop, iquad
INTEGER, INTENT(OUT)       :: ifault
REAL (dp), INTENT(IN)      :: stopcr, simp
REAL (dp), INTENT(IN OUT)  :: p(:), step(:)
REAL (dp), INTENT(OUT)     :: var(:), func

INTERFACE
  SUBROUTINE functn(p, func)
    IMPLICIT NONE
    INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)  :: p(:)
    REAL (dp), INTENT(OUT) :: func
  END SUBROUTINE functn
END INTERFACE

!     Local variables

REAL (dp)   :: g(nop+1,nop), h(nop+1), pbar(nop), pstar(nop), pstst(nop), &
               aval(nop), pmin(nop), temp(nop), bmat(nop*(nop+1)/2),  &
               vc(nop*(nop+1)/2), ymin, rmax, hstst, a0, hmin, test, hmean, &
               hstd, hstar, hmax, savemn

REAL (dp), PARAMETER :: zero = 0._dp, half = 0.5_dp, one = 1._dp, two = 2._dp
INTEGER     :: i, i1, i2, iflag, ii, ij, imax, imin, irank, irow, j, j1, jj, &
               k, l, loop, nap, neval, nmore, np1, nullty

!     A = REFLECTION COEFFICIENT, B = CONTRACTION COEFFICIENT, AND
!     C = EXPANSION COEFFICIENT.

REAL (dp), PARAMETER :: a = 1._dp, b = 0.5_dp, c = 2._dp

!     SET LOUT = LOGICAL UNIT NO. FOR OUTPUT

INTEGER, PARAMETER :: lout = 6

!     IF PROGRESS REPORTS HAVE BEEN REQUESTED, PRINT HEADING

IF (iprint > 0) WRITE (lout,5000) iprint

!     CHECK INPUT ARGUMENTS

ifault = 0
IF (nop <= 0) ifault = 3
IF (nloop <= 0) ifault = 4
IF (ifault /= 0) RETURN

!     SET NAP = NO. OF PARAMETERS TO BE VARIED, I.E. WITH STEP /= 0

nap = COUNT(step /= zero)
neval = 0
loop = 0
iflag = 0

!     IF NAP = 0 EVALUATE FUNCTION AT THE STARTING POINT AND RETURN

IF (nap <= 0) THEN
  CALL functn(p,func)
  RETURN
END IF

!     SET UP THE INITIAL SIMPLEX

20 g(1,:) = p
irow = 2
DO i = 1, nop
  IF (step(i) /= zero) THEN
    g(irow,:) = p
    g(irow,i) = p(i) + step(i)
    irow = irow + 1
  END IF
END DO

np1 = nap + 1
DO i = 1, np1
  p = g(i,:)
  CALL functn(p,h(i))
  neval = neval + 1
  IF (iprint > 0) THEN
    WRITE (lout,5100) neval, h(i), p
  END IF
END DO

!     START OF MAIN CYCLE.

!     FIND MAX. & MIN. VALUES FOR CURRENT SIMPLEX (HMAX & HMIN).

Main_loop: DO
  loop = loop + 1
  imax = 1
  imin = 1
  hmax = h(1)
  hmin = h(1)
  DO i = 2, np1
    IF (h(i) > hmax) THEN
      imax = i
      hmax = h(i)
    ELSE
      IF (h(i) < hmin) THEN
        imin = i
        hmin = h(i)
      END IF
    END IF
  END DO

!     FIND THE CENTROID OF THE VERTICES OTHER THAN P(IMAX)

  pbar = zero
  DO i = 1, np1
    IF (i /= imax) THEN
      pbar = pbar + g(i,:)
    END IF
  END DO
  pbar = pbar / nap

!     REFLECT MAXIMUM THROUGH PBAR TO PSTAR,
!     HSTAR = FUNCTION VALUE AT PSTAR.

  pstar = a * (pbar - g(imax,:)) + pbar
  CALL functn(pstar,hstar)
  neval = neval + 1
  IF (iprint > 0) THEN
    IF (MOD(neval,iprint) == 0) WRITE (lout,5100) neval, hstar, pstar
  END IF

!     IF HSTAR < HMIN, REFLECT PBAR THROUGH PSTAR,
!     HSTST = FUNCTION VALUE AT PSTST.

  IF (hstar < hmin) THEN
    pstst = c * (pstar - pbar) + pbar
    CALL functn(pstst,hstst)
    neval = neval + 1
    IF (iprint > 0) THEN
      IF (MOD(neval,iprint) == 0) WRITE (lout,5100) neval, hstst, pstst
    END IF

!     IF HSTST < HMIN REPLACE CURRENT MAXIMUM POINT BY PSTST AND
!     HMAX BY HSTST, THEN TEST FOR CONVERGENCE.

    IF (hstst >= hmin) THEN   ! REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR.
      g(imax,:) = pstar
      h(imax) = hstar
    ELSE
      g(imax,:) = pstst
      h(imax) = hstst
    END IF
    GO TO 250
  END IF

!     HSTAR IS NOT < HMIN.
!     TEST WHETHER IT IS < FUNCTION VALUE AT SOME POINT OTHER THAN
!     P(IMAX).   IF IT IS REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR.

  DO i = 1, np1
    IF (i /= imax) THEN
      IF (hstar < h(i)) THEN  ! REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR.
        g(imax,:) = pstar
        h(imax) = hstar
        GO TO 250
      END IF
    END IF
  END DO

!     HSTAR > ALL FUNCTION VALUES EXCEPT POSSIBLY HMAX.
!     IF HSTAR <= HMAX, REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR.

  IF (hstar <= hmax) THEN
    g(imax,:) = pstar
    hmax = hstar
    h(imax) = hstar
  END IF

!     CONTRACTED STEP TO THE POINT PSTST,
!     HSTST = FUNCTION VALUE AT PSTST.

  pstst = b * g(imax,:) + (one-b) * pbar
  CALL functn(pstst,hstst)
  neval = neval + 1
  IF (iprint > 0) THEN
    IF (MOD(neval,iprint) == 0) WRITE (lout,5100) neval, hstst, pstst
  END IF

!     IF HSTST < HMAX REPLACE P(IMAX) BY PSTST & HMAX BY HSTST.

  IF (hstst <= hmax) THEN
    g(imax,:) = pstst
    h(imax) = hstst
    GO TO 250
  END IF

!     HSTST > HMAX.
!     SHRINK THE SIMPLEX BY REPLACING EACH POINT, OTHER THAN THE CURRENT
!     MINIMUM, BY A POINT MID-WAY BETWEEN ITS CURRENT POSITION AND THE
!     MINIMUM.

  DO i = 1, np1
    IF (i /= imin) THEN
      DO j = 1, nop
        IF (step(j) /= zero) g(i,j) = (g(i,j) + g(imin,j)) * half
        p(j) = g(i,j)
      END DO
      CALL functn(p,h(i))
      neval = neval + 1
      IF (iprint > 0) THEN
        IF (MOD(neval,iprint) == 0) WRITE (lout,5100) neval, h(i), p
      END IF
    END IF
  END DO

!     IF LOOP = NLOOP TEST FOR CONVERGENCE, OTHERWISE REPEAT MAIN CYCLE.

  250 IF (loop < nloop) CYCLE Main_loop

!     CALCULATE MEAN & STANDARD DEVIATION OF FUNCTION VALUES FOR THE
!     CURRENT SIMPLEX.

  hmean = SUM( h(1:np1) ) / np1
  hstd = SUM( (h(1:np1) - hmean) ** 2 )
  hstd = SQRT(hstd / np1)

!     IF THE RMS > STOPCR, SET IFLAG & LOOP TO ZERO AND GO TO THE
!     START OF THE MAIN CYCLE AGAIN.

  IF (hstd > stopcr .AND. neval <= maxfn) THEN
    iflag = 0
    loop = 0
    CYCLE Main_loop
  END IF

!     FIND THE CENTROID OF THE CURRENT SIMPLEX AND THE FUNCTION VALUE THERE.

  DO i = 1, nop
    IF (step(i) /= zero) THEN
      p(i) = SUM( g(1:np1,i) ) / np1
    END IF
  END DO
  CALL functn(p,func)
  neval = neval + 1
  IF (iprint > 0) THEN
    IF (MOD(neval,iprint) == 0) WRITE (lout,5100) neval, func, p
  END IF

!     TEST WHETHER THE NO. OF FUNCTION VALUES ALLOWED, maxfn, HAS BEEN
!     OVERRUN; IF SO, EXIT WITH IFAULT = 1.

  IF (neval > maxfn) THEN
    ifault = 1
    IF (iprint < 0) RETURN
    WRITE (lout,5200) maxfn
    WRITE (lout,5300) hstd
    WRITE (lout,5400) p
    WRITE (lout,5500) func
    RETURN
  END IF

!     CONVERGENCE CRITERION SATISFIED.
!     IF IFLAG = 0, SET IFLAG & SAVE HMEAN.
!     IF IFLAG = 1 & CHANGE IN HMEAN <= STOPCR THEN SEARCH IS COMPLETE.

  IF (iprint >= 0) THEN
    WRITE (lout,5600)
    WRITE (lout,5400) p
    WRITE (lout,5500) func
  END IF

  IF (iflag == 0 .OR. ABS(savemn-hmean) >= stopcr) THEN
    iflag = 1
    savemn = hmean
    loop = 0
  ELSE
    EXIT Main_loop
  END IF

END DO Main_loop

IF (iprint >= 0) THEN
  WRITE (lout,5700) neval
  WRITE (lout,5800) p
  WRITE (lout,5900) func
END IF
IF (iquad <= 0) RETURN

!------------------------------------------------------------------

!     QUADRATIC SURFACE FITTING

IF (iprint >= 0) WRITE (lout,6000)

!     EXPAND THE FINAL SIMPLEX, IF NECESSARY, TO OVERCOME ROUNDING
!     ERRORS.

hmin = func
nmore = 0
DO i = 1, np1
  DO
    test = ABS(h(i)-func)
    IF (test < simp) THEN
      DO j = 1, nop
        IF (step(j) /= zero) g(i,j) = (g(i,j)-p(j)) + g(i,j)
        pstst(j) = g(i,j)
      END DO
      CALL functn(pstst,h(i))
      nmore = nmore + 1
      neval = neval + 1
      IF (h(i) >= hmin) CYCLE
      hmin = h(i)
      IF (iprint >= 0) WRITE (lout,5100) neval, hmin, pstst
    ELSE
      EXIT
    END IF
  END DO
END DO

!     FUNCTION VALUES ARE CALCULATED AT AN ADDITIONAL NAP POINTS.

DO i = 1, nap
  i1 = i + 1
  pstar = (g(1,:) + g(i1,:)) * half
  CALL functn(pstar,aval(i))
  nmore = nmore + 1
  neval = neval + 1
END DO

!     THE MATRIX OF ESTIMATED SECOND DERIVATIVES IS CALCULATED AND ITS
!     LOWER TRIANGLE STORED IN BMAT.

a0 = h(1)
DO i = 1, nap
  i1 = i - 1
  i2 = i + 1
  DO j = 1, i1
    j1 = j + 1
    pstst = (g(i2,:) + g(j1,:)) * half
    CALL functn(pstst,hstst)
    nmore = nmore + 1
    neval = neval + 1
    l = i * (i-1) / 2 + j
    bmat(l) = two * (hstst + a0 - aval(i) - aval(j))
  END DO
END DO

l = 0
DO i = 1, nap
  i1 = i + 1
  l = l + i
  bmat(l) = two * (h(i1) + a0 - two*aval(i))
END DO

!     THE VECTOR OF ESTIMATED FIRST DERIVATIVES IS CALCULATED AND
!     STORED IN AVAL.

DO i = 1, nap
  i1 = i + 1
  aval(i) = two * aval(i) - (h(i1) + 3._dp*a0) * half
END DO

!     THE MATRIX Q OF NELDER & MEAD IS CALCULATED AND STORED IN G.

pmin = g(1,:)
DO i = 1, nap
  i1 = i + 1
  g(i1,:) = g(i1,:) - g(1,:)
END DO

DO i = 1, nap
  i1 = i + 1
  g(i,:) = g(i1,:)
END DO

!     INVERT BMAT

CALL syminv(bmat, nap, bmat, temp, nullty, ifault, rmax)
IF (ifault == 0) THEN
  irank = nap - nullty
ELSE                                 ! BMAT not +ve definite
                                     ! Resume search for the minimum
  IF (iprint >= 0) WRITE (lout,6100)
  ifault = 2
  IF (neval > maxfn) RETURN
  WRITE (lout,6200)
  step = half * step
  GO TO 20
END IF

!     BMAT*A/2 IS CALCULATED AND STORED IN H.

DO i = 1, nap
  h(i) = zero
  DO j = 1, nap
    IF (j <= i) THEN
      l = i * (i-1) / 2 + j
    ELSE
      l = j * (j-1) / 2 + i
    END IF
    h(i) = h(i) + bmat(l) * aval(j)
  END DO
END DO

!     FIND THE POSITION, PMIN, & VALUE, YMIN, OF THE MINIMUM OF THE
!     QUADRATIC.

ymin = DOT_PRODUCT( h(1:nap), aval(1:nap) )
ymin = a0 - ymin
DO i = 1, nop
  pstst(i) = DOT_PRODUCT( h(1:nap), g(1:nap,i) )
END DO
pmin = pmin - pstst
IF (iprint >= 0) THEN
  WRITE (lout,6300) ymin, pmin
  WRITE (lout,6400)
END IF

!     Q*BMAT*Q'/2 IS CALCULATED & ITS LOWER TRIANGLE STORED IN VC

DO i = 1, nop
  DO j = 1, nap
    h(j) = zero
    DO k = 1, nap
      IF (k <= j) THEN
        l = j * (j-1) / 2 + k
      ELSE
        l = k * (k-1) / 2 + j
      END IF
      h(j) = h(j) + bmat(l) * g(k,i) * half
    END DO
  END DO

  DO j = i, nop
    l = j * (j-1) / 2 + i
    vc(l) = DOT_PRODUCT( h(1:nap), g(1:nap,j) )
  END DO
END DO

!     THE DIAGONAL ELEMENTS OF VC ARE COPIED INTO VAR.

j = 0
DO i = 1, nop
  j = j + i
  var(i) = vc(j)
END DO
IF (iprint < 0) RETURN
WRITE (lout,6500) irank
CALL print_tri_matrix(vc, nop, lout)

WRITE (lout,6600)
CALL syminv(vc, nap, bmat, temp, nullty, ifault, rmax)

!     BMAT NOW CONTAINS THE INFORMATION MATRIX

WRITE (lout,6700)
CALL print_tri_matrix(bmat, nop, lout)

ii = 0
ij = 0
DO i = 1, nop
  ii = ii + i
  IF (vc(ii) > zero) THEN
    vc(ii) = one / SQRT(vc(ii))
  ELSE
    vc(ii) = zero
  END IF
  jj = 0
  DO j = 1, i - 1
    jj = jj + j
    ij = ij + 1
    vc(ij) = vc(ij) * vc(ii) * vc(jj)
  END DO
  ij = ij + 1
END DO

WRITE (lout,6800)
ii = 0
DO i = 1, nop
  ii = ii + i
  IF (vc(ii) /= zero) vc(ii) = one
END DO
CALL print_tri_matrix(vc, nop, lout)

!     Exit, on successful termination.

WRITE (lout,6900) nmore
RETURN

5000 FORMAT (' Progress Report every',i4,' function evaluations'/  &
             ' EVAL.   FUNC.VALUE.          PARAMETER VALUES')
5100 FORMAT (/' ', i4, '  ', g12.5, '  ', 5G11.4, 3(/t22, 5G11.4))
5200 FORMAT (' No. of function evaluations > ',i5)
5300 FORMAT (' RMS of function values of last simplex =', g14.6)
5400 FORMAT (' Centroid of last simplex =',4(/' ', 6G13.5))
5500 FORMAT (' Function value at centroid =', g14.6)
5600 FORMAT (/' EVIDENCE OF CONVERGENCE')
5700 FORMAT (/' Minimum found after',i5,' function evaluations')
5800 FORMAT (' Minimum at',4(/' ', 6G13.6))
5900 FORMAT (' Function value at minimum =', g14.6)
6000 FORMAT (/' Fitting quadratic surface about supposed minimum'/)
6100 FORMAT (/' MATRIX OF ESTIMATED SECOND DERIVATIVES NOT +VE DEFN.'/ &
             ' MINIMUM PROBABLY NOT FOUND'/)
6200 FORMAT (/t11, 'Search restarting'/)
6300 FORMAT (' Minimum of quadratic surface =',g14.6,' at',4(/' ', 6G13.5))
6400 FORMAT (' IF THIS DIFFERS BY MUCH FROM THE MINIMUM ESTIMATED ',      &
             'FROM THE MINIMIZATION,'/' THE MINIMUM MAY BE FALSE &/OR THE ',  &
             'INFORMATION MATRIX MAY BE INACCURATE'/)
6500 FORMAT (' Rank of information matrix =',i3/ &
             ' Inverse of information matrix:-')
6600 FORMAT (/' If the function minimized was -LOG(LIKELIHOOD),'/         &
             ' this is the covariance matrix of the parameters.'/         &
             ' If the function was a sum of squares of residuals,'/       &
             ' this matrix must be multiplied by twice the estimated ',   &
             'residual variance'/' to obtain the covariance matrix.'/)
6700 FORMAT (' INFORMATION MATRIX:-'/)
6800 FORMAT (/' CORRELATION MATRIX:-')
6900 FORMAT (/' A further',i4,' function evaluations have been used'/)

END SUBROUTINE minim

SUBROUTINE syminv(a, n, c, w, nullty, ifault, rmax)

!     ALGORITHM AS7, APPLIED STATISTICS, VOL.17, 1968.

!     ARGUMENTS:-
!     A()    = INPUT, THE SYMMETRIC MATRIX TO BE INVERTED, STORED IN
!                LOWER TRIANGULAR FORM
!     N      = INPUT, ORDER OF THE MATRIX
!     C()    = OUTPUT, THE INVERSE OF A (A GENERALIZED INVERSE IF C IS
!                SINGULAR), ALSO STORED IN LOWER TRIANGULAR FORM.
!                C AND A MAY OCCUPY THE SAME LOCATIONS.
!     W()    = WORKSPACE, DIMENSION AT LEAST N.
!     NULLTY = OUTPUT, THE RANK DEFICIENCY OF A.
!     IFAULT = OUTPUT, ERROR INDICATOR
!                 = 1 IF N < 1
!                 = 2 IF A IS NOT +VE SEMI-DEFINITE
!                 = 0 OTHERWISE
!     RMAX   = OUTPUT, APPROXIMATE BOUND ON THE ACCURACY OF THE DIAGONAL
!                ELEMENTS OF C.  E.G. IF RMAX = 1.E-04 THEN THE DIAGONAL
!                ELEMENTS OF C WILL BE ACCURATE TO ABOUT 4 DEC. DIGITS.

!     LATEST REVISION - 1 April 1985

!***************************************************************************

REAL (dp), INTENT(IN OUT) :: a(:), c(:), w(:)
INTEGER, INTENT(IN)       :: n
INTEGER, INTENT(OUT)      :: nullty, ifault
REAL (dp), INTENT(OUT)    :: rmax

REAL (dp), PARAMETER :: zero = 0._dp, one = 1._dp
INTEGER              :: i, icol, irow, j, jcol, k, l, mdiag, ndiag, nn, nrow
REAL (dp)            :: x

nrow = n
ifault = 1
IF (nrow > 0) THEN
  ifault = 0

!     CHOLESKY FACTORIZATION OF A, RESULT IN C

  CALL chola(a, nrow, c, nullty, ifault, rmax, w)
  IF (ifault == 0) THEN

!     INVERT C & FORM THE PRODUCT (CINV)'*CINV, WHERE CINV IS THE INVERSE
!     OF C, ROW BY ROW STARTING WITH THE LAST ROW.
!     IROW = THE ROW NUMBER, NDIAG = LOCATION OF LAST ELEMENT IN THE ROW.

    nn = nrow * (nrow+1) / 2
    irow = nrow
    ndiag = nn
    10 IF (c(ndiag) /= zero) THEN
      l = ndiag
      DO i = irow, nrow
        w(i) = c(l)
        l = l + i
      END DO
      icol = nrow
      jcol = nn
      mdiag = nn

      30 l = jcol
      x = zero
      IF (icol == irow) x = one / w(irow)
      k = nrow
      40 IF (k /= irow) THEN
        x = x - w(k) * c(l)
        k = k - 1
        l = l - 1
        IF (l > mdiag) l = l - k + 1
        GO TO 40
      END IF

      c(l) = x / w(irow)
      IF (icol == irow) GO TO 60
      mdiag = mdiag - icol
      icol = icol - 1
      jcol = jcol - 1
      GO TO 30
    END IF ! (c(ndiag) /= zero)

    l = ndiag
    DO j = irow, nrow
      c(l) = zero
      l = l + j
    END DO

    60 ndiag = ndiag - irow
    irow = irow - 1
    IF (irow /= 0) GO TO 10
  END IF
END IF
RETURN

END SUBROUTINE syminv

SUBROUTINE chola(a, n, u, nullty, ifault, rmax, r)

!     ALGORITHM AS6, APPLIED STATISTICS, VOL.17, 1968, WITH
!     MODIFICATIONS BY A.J.MILLER

!     ARGUMENTS:-
!     A()    = INPUT, A +VE DEFINITE MATRIX STORED IN LOWER-TRIANGULAR
!                FORM.
!     N      = INPUT, THE ORDER OF A
!     U()    = OUTPUT, A LOWER TRIANGULAR MATRIX SUCH THAT U*U' = A.
!                A & U MAY OCCUPY THE SAME LOCATIONS.
!     NULLTY = OUTPUT, THE RANK DEFICIENCY OF A.
!     IFAULT = OUTPUT, ERROR INDICATOR
!                 = 1 IF N < 1
!                 = 2 IF A IS NOT +VE SEMI-DEFINITE
!                 = 0 OTHERWISE
!     RMAX   = OUTPUT, AN ESTIMATE OF THE RELATIVE ACCURACY OF THE
!                DIAGONAL ELEMENTS OF U.
!     R()    = OUTPUT, ARRAY CONTAINING BOUNDS ON THE RELATIVE ACCURACY
!                OF EACH DIAGONAL ELEMENT OF U.

!     LATEST REVISION - 1 April 1985

!***************************************************************************

REAL (dp), INTENT(IN)   :: a(:)
INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(OUT)  :: u(:)
INTEGER, INTENT(OUT)    :: nullty, ifault
REAL (dp), INTENT(OUT)  :: rmax, r(:)

!     ETA SHOULD BE SET EQUAL TO THE SMALLEST +VE VALUE SUCH THAT
!     1._dp + ETA IS CALCULATED AS BEING GREATER THAN 1._dp IN THE ACCURACY
!     BEING USED.

REAL (dp), PARAMETER :: eta = EPSILON(1.0_dp), zero = 0._dp
INTEGER              :: i, icol, irow, j, k, l, m
REAL (dp)            :: rsq, w

ifault = 1
IF (n > 0) THEN
  ifault = 2
  nullty = 0
  rmax = eta
  r(1) = eta
  j = 1
  k = 0

!     FACTORIZE COLUMN BY COLUMN, ICOL = COLUMN NO.

  DO  icol = 1, n
    l = 0

!     IROW = ROW NUMBER WITHIN COLUMN ICOL

    DO  irow = 1, icol
      k = k + 1
      w = a(k)
      IF (irow == icol) rsq = (w*eta) ** 2
      m = j
      DO  i = 1, irow
        l = l + 1
        IF (i == irow) EXIT
        w = w - u(l) * u(m)
        IF (irow == icol) rsq = rsq + (u(l)**2*r(i)) ** 2
        m = m + 1
      END DO

      IF (irow == icol) EXIT
      IF (u(l) /= zero) THEN
        u(k) = w / u(l)
      ELSE
        u(k) = zero
        IF (ABS(w) > ABS(rmax*a(k))) GO TO 60
      END IF
    END DO

!     END OF ROW, ESTIMATE RELATIVE ACCURACY OF DIAGONAL ELEMENT.

    rsq = SQRT(rsq)
    IF (ABS(w) > 5.*rsq) THEN
      IF (w < zero) RETURN
      u(k) = SQRT(w)
      r(i) = rsq / w
      IF (r(i) > rmax) rmax = r(i)
    ELSE
      u(k) = zero
      nullty = nullty + 1
    END IF

    j = j + icol
  END DO
  ifault = zero
END IF
60 RETURN

END SUBROUTINE chola

SUBROUTINE print_tri_matrix(a, n, lout)

INTEGER, INTENT(IN)    :: n, lout
REAL (dp), INTENT(IN)  :: a(:)

!     Local variables
INTEGER  :: i, ii, i1, i2, l

l = 1
DO l = 1, n, 6
  ii = l * (l-1) / 2
  DO i = l, n
    i1 = ii + l
    ii = ii + i
    i2 = MIN(ii,i1+5)
    WRITE (lout,'(1X,6G13.5)') a(i1:i2)
  END DO
END DO
RETURN

END SUBROUTINE print_tri_matrix

END MODULE Nelder_Mead

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
program build_cubist_model

integer, parameter :: dx=3750, dy=3750
real :: latn(dx,dy), dtrad(dx,dy), day(dx,dy), night(dx,dy)
real :: dtime(dx,dy), lai(dx,dy), terrain(dx,dy)
character(len=4) :: arg1,arg3
character(len=3) :: arg2
character(len=255) :: arg4, basedir
!character (len=:), allocatable :: basedir
integer :: year, doy, yyyyddd, riseddd
character(len=255) :: ifile1, ifile2, ifile3, ifile4, ifile5, ifile6, ofile1
real :: predict(dx,dy)
real :: precip(dx,dy), fmax(dx,dy)
real :: fc(dx,dy)
integer :: laiddd, nweek1, cday, offset, rday

call getarg(1,arg1)
call getarg(2,arg2)
call getarg(3,arg3)
call getarg(4,arg4)
read(arg1,'(i4)') year
read(arg2,'(i3)') doy
yyyyddd=year*1000+doy
basedir=trim(arg4)

nweek1=(doy-1)/7.
cday=nweek1*7.
offset=int((doy-cday)/7.)
rday=int(((offset+nweek1)*7)+1)
riseddd=2014*1000+rday

write(ifile1,'(a,a,a,a)') basedir,'/STATIC/PRECIP/tiles/PRECIP_',arg3,'.dat'
write(ifile2,'(a,a,a,a)') basedir,'/STATIC/FMAX/tiles/FMAX_',arg3,'.dat'
open(10,file=ifile1,form='unformatted',access='direct',recl=dx*dy*4)
open(11,file=ifile2,form='unformatted',access='direct',recl=dx*dy*4)
read(10,rec=1) precip
read(11,rec=1) fmax
close(10)
close(11)

nweek1=(doy-1)/7.
cday=nweek1*8.
offset=int((doy-cday)/7.)
rday=int(((offset+nweek1)*8)+1)
laiddd=year*1000+rday

write(ifile1,'(a,a,a,a)') basedir,'/STATIC/TERRAIN_SD/tiles/TERRAIN_',arg3,'.dat'
write(ifile2,'(a,a,a,a,I7,a,a,a)') basedir,'/TILES/',arg3,'/FINAL_DAY_LST_',yyyyddd,'_',arg3,'.dat'
write(ifile3,'(a,a,a,a,I7,a,a,a)') basedir,'/TILES/',arg3,'/FINAL_NIGHT_LST_',yyyyddd,'_',arg3,'.dat'
write(ifile4,'(a,a,I7,a,a,a)') basedir,'/STATIC/LAI/tiles/MLAI_',laiddd,'_',arg3,'.dat'
write(ifile5,'(a,a,I7,a,a,a)') basedir,'/STATIC/DTIME/tiles/DTIME_',riseddd,'_',arg3,'.dat'

open(10,file=ifile1,form='unformatted',access='direct',recl=dx*dy*4)
open(11,file=ifile2,form='unformatted',access='direct',recl=dx*dy*4)
open(12,file=ifile3,form='unformatted',access='direct',recl=dx*dy*4)
open(13,file=ifile4,form='unformatted',access='direct',recl=dx*dy*4)
open(14,file=ifile5,form='unformatted',access='direct',recl=dx*dy*4)
read(10,rec=1) terrain
read(11,rec=1) day
read(12,rec=1) night
read(13,rec=1) lai
read(14,rec=1) dtime
close(10)
close(11)
close(12)
close(13)
close(14)

write(6,*) day(3001,2501), night(3001,2501)
latn(:,:) = -9999.
do j = 1, dy
do i = 1, dx

 if (precip(i,j).gt.250.and.precip(i,j).lt.500.and.fmax(i,j).ge.0.0.and.fmax(i,j).lt.1.0) then
 if (day(i,j).ne.-9999.and.night(i,j).ne.-9999.) then
 if (terrain(i,j).ne.-9999.and.lai(i,j).ne.-9999.) then

  dtrad(i,j) = day(i,j)-night(i,j)
  fc(i,j) = 1.0-exp(-0.5*lai(i,j))

if ( dtime(i,j).le.1.5666664.and.dtrad(i,j).le.14.333391) then
 latn(i,j) = (-10.1172069+(0.099*dtrad(i,j))+(0.05*day(i,j))+(-0.005*terrain(i,j))+(0.5*fc(i,j)))
endif
if ( dtime(i,j).le.1.6458384.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (-33.7260704+(0.138*dtrad(i,j))+(0.132*day(i,j)))
endif
if ( dtime(i,j).le.1.7687515.and.dtime(i,j).gt.1.5666664) then
 latn(i,j) = (-8.9760558+(0.158*dtrad(i,j))+(-0.006*terrain(i,j))+(7.7641*dtime(i,j)))
endif
if ( dtrad(i,j).le.14.333391.and.dtime(i,j).gt.1.7687515.and.dtime(i,j).le.2.5999999) then
 latn(i,j) = (-6.0488545+(0.286*dtrad(i,j))+(0.006*day(i,j))+(-2.5*fc(i,j))+(4.3457*dtime(i,j)))
endif
if ( dtime(i,j).le.2.3291681.and.dtrad(i,j).gt.14.333391.and.dtime(i,j).gt.1.6458384) then
 latn(i,j) = (-20.8916886+(0.234*dtrad(i,j))+(0.061*day(i,j))+(-0.004*terrain(i,j))+(-2.6*fc(i,j))+(4.3235*dtime(i,j)))
endif
if ( dtrad(i,j).le.14.333391.and.dtime(i,j).gt.2.5999999) then
 latn(i,j) = (-3.2000537+(0.576*dtrad(i,j))+(0.008*day(i,j))+(0.009*terrain(i,j))+(-1.2*fc(i,j))+(1.2707*dtime(i,j)))
endif
if ( day(i,j).le.296.32886.and.dtime(i,j).gt.2.3291681.and.dtime(i,j).le.3.2041662.and.fc(i,j).gt.0.099774703) then
 latn(i,j) = (-27.6182198+(0.213*dtrad(i,j))+(0.112*day(i,j))+(-0.006*terrain(i,j))+(-4.1*fc(i,j))+(1.5605*dtime(i,j)))
endif
if ( dtime(i,j).le.2.8083334.and.day(i,j).gt.296.32886.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (-19.4747+(0.348*dtrad(i,j))+(0.033*day(i,j))+(-2.5*fc(i,j))+(6.113*dtime(i,j)))
endif
if ( fc(i,j).le.0.099774703.and.dtime(i,j).gt.2.3291681.and.dtime(i,j).le.2.8916676.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (-21.7913334+(0.255*dtrad(i,j))+(0.052*day(i,j))+(-0.007*terrain(i,j))+(-0.6*fc(i,j))+(5.6417*dtime(i,j)))
endif
if ( dtime(i,j).le.2.9916666.and.dtime(i,j).gt.2.8916674.and.fc(i,j).gt.0.099774703) then
 latn(i,j) = (119.1887771+(0.408*dtrad(i,j))+(0.041*day(i,j))+(-0.007*terrain(i,j))+(-3*fc(i,j))+(-42.1144*dtime(i,j)))
endif
if ( dtime(i,j).le.2.8916674.and.dtime(i,j).gt.2.8083334.and.fc(i,j).gt.0.099774703.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (-49.9743535+(0.289*dtrad(i,j))+(0.032*day(i,j))+(-0.006*terrain(i,j))+(-6.3*fc(i,j))+(17.4634*dtime(i,j)))
endif
if ( dtime(i,j).le.3.2041662.and.dtime(i,j).gt.2.9916666.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (37.3755468+(0.355*dtrad(i,j))+(0.011*day(i,j))+(-0.003*terrain(i,j))+(-4.9*fc(i,j))+(-10.4766*dtime(i,j)))
endif
if ( dtime(i,j).gt.4.2833338.and.fc(i,j).gt.0.2488669) then
 latn(i,j) = (45.1926396+(0.439*dtrad(i,j))+(-0.081*day(i,j))+(0.006*terrain(i,j))+(-6.5*fc(i,j))+(-2.3874*dtime(i,j)))
endif
if ( terrain(i,j).le.2.4648063.and.dtime(i,j).gt.3.2041662.and.fc(i,j).gt.0.099774703.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (7.8977438+(0.45*dtrad(i,j))+(-0.027*day(i,j))+(1.197*terrain(i,j))+(-5.3*fc(i,j))+(1.5055*dtime(i,j)))
endif
if ( dtime(i,j).gt.3.3322911.and.dtime(i,j).le.3.8177083.and.day(i,j).le.319.33286.and.fc(i,j).gt.0.099774703.and.terrain(i,j).gt.2.4648063) then
 latn(i,j) = (-15.6814407+(0.328*dtrad(i,j))+(0.047*day(i,j))+(-0.007*terrain(i,j))+(-5*fc(i,j))+(3.4801*dtime(i,j)))
endif
if ( dtime(i,j).le.3.3322911.and.dtime(i,j).gt.3.2041662.and.terrain(i,j).gt.2.4648063) then
 latn(i,j) = (-34.6459961+(0.172*dtrad(i,j))+(0.064*day(i,j))+(-0.017*terrain(i,j))+(-7.5*fc(i,j))+(9.5506*dtime(i,j)))
endif
if ( dtime(i,j).gt.3.8177083.and.dtime(i,j).le.4.2833338.and.day(i,j).le.319.33286.and.fc(i,j).gt.0.099774703.and.terrain(i,j).gt.2.4648063) then
 latn(i,j) = (-12.8867856+(0.326*dtrad(i,j))+(0.028*day(i,j))+(-0.004*terrain(i,j))+(-6.7*fc(i,j))+(4.1068*dtime(i,j)))
endif
if ( fc(i,j).le.0.099774703.and.dtime(i,j).gt.2.8916676.and.dtime(i,j).le.4.2833338) then
 latn(i,j) = (-16.5686158+(0.399*dtrad(i,j))+(0.03*day(i,j))+(-0.003*terrain(i,j))+(9*fc(i,j))+(4.1221*dtime(i,j)))
endif
if ( day(i,j).gt.319.33286.and.fc(i,j).gt.0.099774703.and.dtime(i,j).le.4.2833338.and.dtime(i,j).gt.3.2041662) then
 latn(i,j) = (19.6101441+(0.384*dtrad(i,j))+(-0.059*day(i,j))+(-0.011*terrain(i,j))+(-9.4*fc(i,j))+(2.7591*dtime(i,j)))
endif
if ( dtime(i,j).gt.4.2833338.and.fc(i,j).le.0.2488669.and.dtrad(i,j).gt.14.333391) then
 latn(i,j) = (34.2715956+(0.468*dtrad(i,j))+(-0.045*day(i,j))+(0.009*terrain(i,j))+(-15.7*fc(i,j))+(-2.134*dtime(i,j)))
endif
 endif
 endif
 endif

enddo
enddo

write(ifile1,'(a,a,a,a)') basedir,'/PROCESSING/DTRAD_PREDICTION/comp1_',arg3,'.dat'
open(10,file=ifile1,form='unformatted',access='direct',recl=dx*dy*4)
write(10,rec=1) latn
close(10)

end program

max_stars_repo_name: bucricket/projectMASviirsMAC

max_issues_repo_name: bucricket/projectMASviirsMAC

max_forks_repo_name: bucricket/projectMASviirsMAC


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Find element position on a vector.
!! @author Ilario Mazzieri
!> @date September, 2013 - Creation
!> @version 1.0
!> @param[in] n
!> @param[in] vec vec vector where seraching tar
!> @param[in] tar target
!> @param[ount] ind position of tar in vec, 0 otherwise

      subroutine FIND_POS(n,vec,tar,ind)

      implicit none

      integer*4 :: n, i, tar, ind
      integer*4, dimension(n) :: vec

      ind = 0
      do i = 1, n
        if(vec(i) .eq. tar) then
          ind = i
          return
        endif
      enddo

      return

      end subroutine FIND_POS


MODULE band_unfold_module

  use aa_module, only: aa
  use bb_module, only: bb, calc_bb
  use rgrid_module, only: Ngrid, Igrid
  use band_variables, only: nbk, ak, nfki
  use parallel_module
  use wf_module, only: unk, esp
  use fft_module
  use rsdft_mpi_module
  use bz_module, only: kbb
  use watch_module

  implicit none

  PRIVATE
  PUBLIC :: read_band_unfold, init_band_unfold, band_unfold &
           ,finalize_band_unfold, iswitch_banduf

  logical :: iswitch_banduf=.false.
  logical :: iswitch_ufld_l=.false.

  integer :: unit_uf
  integer :: unit_uf_k
  integer :: unit_uf_z
  integer :: nktrj,nktrj_0
  integer :: nbk_pc
  integer :: mg_pc, mg_sc, mg_pc_0, mg_pc_1
  integer,allocatable :: LG_pc(:,:)
  integer,allocatable :: map_ktrj(:,:)

  real(8) :: ax_pc, aa_pc(3,3), bb_pc(3,3)
  real(8),allocatable :: kbb_pc(:,:),kxyz_pc(:,:)
  real(8),allocatable :: kbb_0(:,:),kxyz_0(:,:),kbb_1(:,:)
  real(8),allocatable :: weight_uf(:,:,:)
  real(8),allocatable :: weight_uf_k(:,:,:)
  real(8),allocatable :: weight_uf_z(:,:,:,:)

#ifdef _DRSDFT_
  integer,parameter :: TYPE_MAIN=MPI_REAL8
  real(8),parameter :: zero=0.0d0
#else
  integer,parameter :: TYPE_MAIN=RSDFT_MPI_COMPLEX16
  complex(8),parameter :: zero=(0.0d0,0.0d0)
#endif

  integer :: job_ctrl=0

CONTAINS

  SUBROUTINE read_band_unfold( rank, unit )
    implicit none
    integer,intent(IN) :: rank,unit
    character(6) :: cbuf,ckey
    integer :: i
    call write_border( 0, " read_band_unfold(start)" )
    if ( rank == 0 ) then
       rewind unit
       do i=1,10000
          read(unit,*,END=999) cbuf
          call convert_capital(cbuf,ckey)
          if ( ckey(1:6) == "BANDUF" ) then
             backspace(unit)
             read(unit,*) cbuf
             read(unit,*) ax_pc
             read(unit,*) aa_pc(1:3,1)
             read(unit,*) aa_pc(1:3,2)
             read(unit,*) aa_pc(1:3,3)
             iswitch_banduf=.true.
          else if ( ckey(1:6) == "UFLD_L" ) then
             iswitch_ufld_l=.true.
          end if
       end do
999    continue
       write(*,*) "iswitch_banduf=",iswitch_banduf
    end if
    call mpi_bcast(iswitch_banduf,1,MPI_LOGICAL,0,MPI_COMM_WORLD,i)
    call mpi_bcast(iswitch_ufld_l,1,MPI_LOGICAL,0,MPI_COMM_WORLD,i)
    call mpi_bcast(ax_pc,1,MPI_REAL8,0,MPI_COMM_WORLD,i)
    call mpi_bcast(aa_pc,9,MPI_REAL8,0,MPI_COMM_WORLD,i)
    call write_border( 0, " read_band_unfold(end)" )
  END SUBROUTINE read_band_unfold

  SUBROUTINE init_band_unfold &
       ( nktrj_io, ktrj_out, unit_in, disp_switch, job_ctrl_in )

    implicit none
    integer,intent(IN)  :: unit_in
    logical,intent(IN)  :: disp_switch
    integer,intent(INOUT) :: nktrj_io
    real(8),intent(OUT) :: ktrj_out(6,nktrj_io)
    integer,optional,intent(IN) :: job_ctrl_in
    integer :: iktrj,MB,MS,MZ

    call write_border( 0, " init_band_unfold(start)" )

    unit_uf   = unit_in
    unit_uf_k = unit_in + 1
    unit_uf_z = unit_in + 2
    if ( myrank == 0 ) then
       open(unit_uf  ,file="band_ufld")
       open(unit_uf_k,file="band_ufld_k")
       if ( iswitch_ufld_l ) open(unit_uf_z,file="band_ufld_z")
    end if

    if ( present(job_ctrl_in) ) job_ctrl=job_ctrl_in

    aa_pc(:,:) = ax_pc*aa_pc(:,:)

    call calc_bb( aa_pc, bb_pc )

    nktrj_0 = nktrj_io

    call init0_band_unfold

!    call possible_kpoints_in_PC( kbb_1 )

    ktrj_out=0.0d0
    ktrj_out(1:3,1:nktrj) = kbb_1(1:3,1:nktrj)

    nktrj_io = nktrj

    call init1_band_unfold ! generate G-vectors of PC

    MB=sum(ir_band)
    MS=sum(ir_spin)
    MZ=Ngrid(3)
    allocate( weight_uf(nktrj_0,MB,MS) ) ; weight_uf=0.0d0
    allocate( weight_uf_k(nktrj_0,MB,MS) ) ; weight_uf_k=0.0d0

    if ( iswitch_ufld_l ) then
       allocate( weight_uf_z(nktrj_0,MB,MS,0:MZ-1) ) ; weight_uf_z=0.0d0
       if ( disp_switch ) write(*,*) &
            "size[weight_uf_z](MB)=",size(weight_uf_z)*8.0d0/1024.d0**2
    end if

    if ( disp_switch ) then
       write(*,*) "size[weight_uf](MB)=",size(weight_uf)*8.0d0/1024.d0**2
       write(*,*) "mg_pc=",mg_pc
       write(*,*) "(aa_pc)"
       write(*,'(1x, f20.15)') ax_pc
       write(*,'(1x,3f20.15)') aa_pc(1:3,1)
       write(*,'(1x,3f20.15)') aa_pc(1:3,2)
       write(*,'(1x,3f20.15)') aa_pc(1:3,3)
       write(*,*) "(bb_pc)"
       write(*,'(1x,3f20.15)') bb_pc(1:3,1)
       write(*,'(1x,3f20.15)') bb_pc(1:3,2)
       write(*,'(1x,3f20.15)') bb_pc(1:3,3)
       write(*,'(1x,a4,1x,2(10x,a10,10x))') "iktrj","k(SC)","k(PC)"
       do iktrj=1,nktrj_0
          write(*,'(1x,i4,1x,2(3f10.6,1x))') &
               iktrj,kbb_0(1:3,iktrj),kbb_pc(1:3,iktrj)
       end do
       write(*,'(1x,a4,1x,10x,a10,10x,a4)') "iktrj","k'(SC)"," "
       do iktrj=1,nktrj
          write(*,'(1x,i4,1x,3f10.6,1x,i4)') &
               iktrj,kbb_1(1:3,iktrj),count(map_ktrj(:,1)==iktrj)
       end do
    end if

    call write_border( 0, " init_band_unfold(end)" )

  END SUBROUTINE init_band_unfold

  SUBROUTINE possible_kpoints_in_PC( ksc )
    implicit none
    real(8),intent(IN) :: ksc(:,:)
    integer :: nksc,i1,i2,i3,k
    real(8) :: q1,q2,q3,qxyz_sc(3),pi2,kpc(3)
    real(8),allocatable :: kpc_0(:,:)
    pi2=2.0d0*acos(-1.0d0)
    nksc=size(ksc,2)
    do k=1,nksc
       do i3=-(Ngrid(3)-1)/2,Ngrid(3)/2
       do i2=-(Ngrid(2)-1)/2,Ngrid(2)/2
       do i1=-(Ngrid(1)-1)/2,Ngrid(1)/2
          q1 = ksc(1,k) + i1
          q2 = ksc(2,k) + i2
          q3 = ksc(3,k) + i3
          qxyz_sc(1) = bb(1,1)*q1 + bb(1,2)*q2 + bb(1,3)*q3
          qxyz_sc(2) = bb(2,1)*q1 + bb(2,2)*q2 + bb(2,3)*q3
          qxyz_sc(3) = bb(3,1)*q1 + bb(3,2)*q2 + bb(3,3)*q3
          kpc(1) = sum( aa_pc(:,1)*qxyz_sc(:) )/pi2
          kpc(2) = sum( aa_pc(:,2)*qxyz_sc(:) )/pi2
          kpc(3) = sum( aa_pc(:,3)*qxyz_sc(:) )/pi2
          call fold_into_primitive( kpc )
!          if ( disp_switch_parallel ) write(*,'(1x,3f10.5)') kpc(:)
          call store_independent_vectors( kpc, kpc_0 )
       end do
       end do
       end do
    end do ! k
    do k=1,size(kpc_0,2)
       if ( disp_switch_parallel ) write(*,'(1x,i8,3f15.5)') k,kpc_0(:,k)
    end do
    deallocate( kpc_0 )
    call end_mpi_parallel
    stop
  END SUBROUTINE possible_kpoints_in_PC

  SUBROUTINE store_independent_vectors( v, v0 )
    implicit none
    real(8),intent(IN) :: v(:)
    real(8),allocatable,intent(INOUT) :: v0(:,:)
    real(8),allocatable :: vtmp(:,:)
    integer :: m,n,i
    m=size( v )
    if ( .not.allocated(v0) ) then
       allocate( v0(m,1) ) ; v0=0.0d0
       v0(:,1)=v(:)
    else
       n=size( v0, 2 )
       do i=1,n
          if ( all( abs(v-v0(:,i)) < 1.d-6 ) ) return
       end do
       allocate( vtmp(m,n) ) ; vtmp=0.0d0
       vtmp=v0
       deallocate( v0 )
       allocate( v0(m,n+1) ) ; v0=0.0d0
       v0(:,1:n)=vtmp
       v0(:,n+1)=v
       if ( disp_switch_parallel ) write(*,'(1x,i4,3f15.5)') n+1,v
       deallocate( vtmp )
    end if
  END SUBROUTINE store_independent_vectors

  SUBROUTINE fold_into_primitive( v )
    implicit none
    real(8),intent(INOUT) :: v(3)
    integer,parameter :: maxloop=100000
    integer :: i,loop
    do i=1,3
       do loop=1,maxloop
          if ( v(i) > 0.5001d0 ) then
             v(i) = v(i) - 1.0d0
          else if ( v(i) <= -0.4999d0 ) then
             v(i) = v(i) + 1.0d0
          else
             exit
          end if
       end do ! loop
       if ( loop > maxloop ) then
          write(*,*) i,v(i)
          call end_mpi_parallel
          stop "error!"
       end if
    end do
  END SUBROUTINE fold_into_primitive

  SUBROUTINE init0_band_unfold

    implicit none
    real(8) :: pi2,x
    integer :: ibk,j,i,iktrj,loop

! --- k points in PC ---

    allocate( kbb_pc(3,nktrj_0)  ) ; kbb_pc=0.0d0
    allocate( kxyz_pc(3,nktrj_0) ) ; kxyz_pc=0.0d0

    if ( job_ctrl == 2 ) then

       kbb_pc(:,:) = kbb(:,:)

    else

       j=0
       do ibk=1,nbk
          do i=0,nfki(ibk)-1
             j=j+1
             kbb_pc(1,j) = ak(1,ibk) + i*( ak(1,ibk+1) - ak(1,ibk) )/nfki(ibk)
             kbb_pc(2,j) = ak(2,ibk) + i*( ak(2,ibk+1) - ak(2,ibk) )/nfki(ibk)
             kbb_pc(3,j) = ak(3,ibk) + i*( ak(3,ibk+1) - ak(3,ibk) )/nfki(ibk)
          end do ! i
       end do ! ibk
       kbb_pc(1,nktrj_0) = ak(1,nbk+1)
       kbb_pc(2,nktrj_0) = ak(2,nbk+1)
       kbb_pc(3,nktrj_0) = ak(3,nbk+1)

    end if

    do iktrj=1,nktrj_0
       kxyz_pc(1,iktrj) = sum( bb_pc(1,:)*kbb_pc(:,iktrj) )
       kxyz_pc(2,iktrj) = sum( bb_pc(2,:)*kbb_pc(:,iktrj) )
       kxyz_pc(3,iktrj) = sum( bb_pc(3,:)*kbb_pc(:,iktrj) )
    end do

    if ( disp_switch_parallel ) then
       x=0.0d0
       write(*,'(1x,i5,4f15.10)') 1,kbb_pc(:,1),x
       do iktrj=2,nktrj_0
          x=x+sqrt( sum((kxyz_pc(:,iktrj)-kxyz_pc(:,iktrj-1))**2) )
          write(*,'(1x,i5,4f15.10)') iktrj, kbb_pc(:,iktrj), x
       end do
    end if

! --- corresponding k-points in SC ---

    pi2 = 2.0d0*acos(-1.0d0)

    allocate( kbb_0(3,nktrj_0)   ) ; kbb_0=0.0d0
    allocate( kxyz_0(3,nktrj_0)  ) ; kxyz_0=0.0d0

    do iktrj=1,nktrj_0
       kbb_0(1,iktrj) = sum( aa(:,1)*kxyz_pc(:,iktrj) )/pi2
       kbb_0(2,iktrj) = sum( aa(:,2)*kxyz_pc(:,iktrj) )/pi2
       kbb_0(3,iktrj) = sum( aa(:,3)*kxyz_pc(:,iktrj) )/pi2
    end do

    do iktrj=1,nktrj_0
       do i=1,3
          do loop=1,1000
             if ( kbb_0(i,iktrj) > 0.5d0 ) then
                kbb_0(i,iktrj) = kbb_0(i,iktrj) - 1.0d0
             else if ( kbb_0(i,iktrj) <= -0.50001d0 ) then
                kbb_0(i,iktrj) = kbb_0(i,iktrj) + 1.0d0
             else
                exit
             end if
          end do ! loop
       end do ! i
       kxyz_0(1:3,iktrj) = matmul( bb(1:3,1:3), kbb_0(1:3,iktrj) )
    end do ! iktrj

! ---

    allocate( kbb_1(3,nktrj_0)    ) ; kbb_1=0.0d0
    allocate( map_ktrj(nktrj_0,2) ) ; map_ktrj=0

    nktrj=1
    kbb_1(1:3,1)=kbb_0(1:3,1)
    map_ktrj(1,1)=1
    map_ktrj(1,2)=0

    do i=2,nktrj_0
       do j=1,nktrj
          if ( all( abs(kbb_0(:,i)-kbb_1(:,j)) < 1.d-12 ) ) then
             map_ktrj(i,1)=j
             map_ktrj(i,2)=0
             exit
          else if ( all( abs(kbb_0(:,i)+kbb_1(:,j)) < 1.d-12 ) ) then
             map_ktrj(i,1)=j
             map_ktrj(i,2)=1
             exit
          end if
       end do
       if ( j > nktrj ) then
          nktrj=nktrj+1
          kbb_1(:,nktrj)=kbb_0(:,i)
          map_ktrj(i,1)=j
          map_ktrj(i,2)=0
       end if
    end do

    do i=1,nktrj
       write(*,'(1x,i4,3f10.5)') i,kbb_1(:,i)
    end do

  END SUBROUTINE init0_band_unfold

  SUBROUTINE init1_band_unfold

    implicit none
    integer :: i,j,i1,i2,i3
    integer,allocatable :: LG_tmp(:,:)
    real(8) :: pi2
    real(8) :: matab(3,3),vtmp(3),utmp(3)
    integer,allocatable :: ir(:),id(:)

    pi2 = 2.0d0*acos(-1.0d0)

    do j=1,3
    do i=1,3
       matab(i,j)=sum( aa_pc(1:3,i)*bb(1:3,j) )/pi2
    end do
    end do
!    do i=1,3
!       write(*,'(1x,3f15.10)') matab(i,1:3)
!    end do

    mg_sc = Ngrid(1)*Ngrid(2)*Ngrid(3)

    allocate( LG_tmp(3,mg_sc) ) ; LG_tmp=0

    mg_pc=0
    do i3=-(Ngrid(3)-1)/2,Ngrid(3)/2
    do i2=-(Ngrid(2)-1)/2,Ngrid(2)/2
    do i1=-(Ngrid(1)-1)/2,Ngrid(1)/2

       vtmp(1) = i1
       vtmp(2) = i2
       vtmp(3) = i3

       utmp(:) = matmul( matab(:,:), vtmp(:) )

       vtmp(:) = nint( utmp(:) )

       if ( all( abs(vtmp-utmp) < 1.d-13 ) ) then
          mg_pc=mg_pc+1
          LG_tmp(1,mg_pc)=i1
          LG_tmp(2,mg_pc)=i2
          LG_tmp(3,mg_pc)=i3
       end if

    end do
    end do
    end do

    allocate( LG_pc(3,mg_pc) ) ; LG_pc=0
    LG_pc(:,:) = LG_tmp(:,1:mg_pc)

    deallocate( LG_tmp )

! ---

    allocate( ir(0:np_grid-1) ) ; ir=0
    allocate( id(0:np_grid-1) ) ; id=0
    do i=1,mg_pc
       j=mod(i-1,np_grid)
       ir(j)=ir(j)+1
    end do
    do i=0,np_grid-1
       id(i)=sum(ir(0:i))-ir(i)
    end do
    mg_pc_0 = id(myrank_g) + 1
    mg_pc_1 = id(myrank_g) + ir(myrank_g)
    deallocate( id, ir )

  END SUBROUTINE init1_band_unfold

  SUBROUTINE band_unfold( jktrj, disp_switch )

    implicit none
    integer,intent(IN) :: jktrj
    logical,intent(IN) :: disp_switch
    complex(8),allocatable :: zwork0(:,:,:),zwork1(:,:,:)
    complex(8),allocatable :: zwork2(:,:),zwork3(:,:)
    integer :: ML,ML1,ML2,ML3,MSP_0,MSP_1,MBZ_0,MBZ_1,MB_0,MB_1,ML_0,ML_1
    integer :: MB,MS,s,k,n,i,i1,i2,i3,j3,ierr,LG_sc(3),iktrj
    integer,allocatable :: iktrj_2_k(:)
    real(8) :: vtmp(3),utmp(3),pi2,sum0
    real(8) :: t1(2),t2(2),tt(2),ttt(2,5)
    real(8) :: dz,dG,z,G
    complex(8) :: ztmp, phase

    if ( .not.iswitch_banduf ) return

    call write_border( 0, " band_unfold(start)" )

    ML  = Ngrid(0)
    ML1 = Ngrid(1)
    ML2 = Ngrid(2)
    ML3 = Ngrid(3)

    ML_0  = id_grid(myrank_g) + 1
    ML_1  = id_grid(myrank_g) + ir_grid(myrank_g)
    MB_0  = id_band(myrank_b) + 1
    MB_1  = id_band(myrank_b) + ir_band(myrank_b)
    MB    = sum(ir_band)
    MBZ_0 = id_bzsm(myrank_k) + 1
    MBZ_1 = id_bzsm(myrank_k) + ir_bzsm(myrank_k)
    MSP_0 = id_spin(myrank_s) + 1
    MSP_1 = id_spin(myrank_s) + ir_spin(myrank_s)
    MS    = sum(ir_spin)

    call init_fft

! ---

    pi2 = 2.0d0*acos(-1.0d0)

    dG = pi2/aa(3,3)

    dz = aa(3,3)/Ngrid(3)

! ---

    allocate( zwork0(0:ML1-1,0:ML2-1,0:ML3-1) ) ; zwork0=(0.0d0,0.0d0)
    allocate( zwork1(0:ML1-1,0:ML2-1,0:ML3-1) ) ; zwork1=(0.0d0,0.0d0)
    allocate( iktrj_2_k(nktrj_0)              ) ; iktrj_2_k=0

    if ( iswitch_ufld_l ) then
       allocate( zwork2(0:ML1-1,0:ML2-1) ) ; zwork2=(0.0d0,0.0d0)
       allocate( zwork3(0:ML1-1,0:ML2-1) ) ; zwork3=(0.0d0,0.0d0)
    end if

    weight_uf(:,:,:)=0.0d0
    weight_uf_k(:,:,:)=0.0d0

    tt=0.0d0

    do iktrj=1,nktrj_0

       if ( map_ktrj(iktrj,1) /= jktrj ) cycle

       call watchb( t1 )

       vtmp(:) = kxyz_0(:,iktrj) - kxyz_pc(:,iktrj)
       utmp(:) = matmul( vtmp, aa(:,:) )/pi2
       vtmp(:) = nint( utmp(:) )

       if ( any( abs(vtmp-utmp) > 1.d-13 ) ) stop "stop@band_unfold"

       LG_sc(:) = nint( utmp(:) )

! ---

       iktrj_2_k(iktrj)=MBZ_0

       do s=MSP_0,MSP_1
       do k=MBZ_0,MBZ_0
       do n=MB_0 ,MB_1

          call watchb( tt )

          zwork0(:,:,:)=zero
          i=ML_0-1
          do i3=Igrid(1,3),Igrid(2,3)
          do i2=Igrid(1,2),Igrid(2,2)
          do i1=Igrid(1,1),Igrid(2,1)
             i=i+1
             zwork0(i1,i2,i3) = unk(i,n,k,s)
          end do
          end do
          end do
          if ( map_ktrj(iktrj,2) == 1 ) zwork0=conjg(zwork0)

          call watchb( tt, ttt(:,1) )

          call rsdft_allreduce_sum( zwork0, comm_grid )

          call watchb( tt, ttt(:,2) )

          if ( iswitch_ufld_l ) then
#ifdef TEST
             do i2=0,ML1-1
             do i1=0,ML1-1
                do j3=0,ML3-1 !(R)
                   z = j3*dz
                   ztmp=zero
                   do i3=0,ML3-1 !(G)
                      G = i3*dG
                      phase = dcmplx( cos(G*z), sin(G*z) )
                      ztmp = ztmp + zwork0(i1,i2,i3)*phase
                   end do !(G)
                   zwork1(i1,i2,j3) = ztmp
                end do !(R)
             end do ! i1
             end do ! i2
             sum0=0.0d0
             do j3=0,ML3-1
                do i=1,mg_pc
                   i1 = mod( LG_pc(1,i)-LG_sc(1)+ML1, ML1 )
                   i2 = mod( LG_pc(2,i)-LG_sc(2)+ML2, ML2 )
                   sum0=sum0+abs( zwork1(i1,i2,j3) )**2
                end do
                weight_uf_z(iktrj,n,s,j3) = sum0
             end do
#endif
             do i3=0,ML3-1
                zwork2(:,:) = zwork0(:,:,i3)
                call forward_2d_fft( zwork2, zwork3 )
                sum0=0.0d0
                do i=1,mg_pc
                   i1 = mod( LG_pc(1,i)-LG_sc(1)+ML1, ML1 )
                   i2 = mod( LG_pc(2,i)-LG_sc(2)+ML2, ML2 )
                   sum0=sum0+abs( zwork2(i1,i2) )**2
                end do
                weight_uf_z(iktrj,n,s,i3) = sum0
             end do
          end if ! iswitch_ufld_l

          call watchb( tt )

          call forward_fft( zwork0, zwork1 )

          call watchb( tt, ttt(:,3) )

          sum0=0.0d0
          do i=mg_pc_0,mg_pc_1
             i1 = mod( LG_pc(1,i)-LG_sc(1)+ML1, ML1 )
             i2 = mod( LG_pc(2,i)-LG_sc(2)+ML2, ML2 )
             i3 = mod( LG_pc(3,i)-LG_sc(3)+ML3, ML3 )
             sum0=sum0+abs( zwork0(i1,i2,i3) )**2
          end do
          call rsdft_allreduce_sum( sum0, comm_grid )
          weight_uf(iktrj,n,s)=sum0

          call watchb( tt, ttt(:,4) )

          i1 = mod( -LG_sc(1)+ML1, ML1 )
          i2 = mod( -LG_sc(2)+ML2, ML2 )
          i3 = mod( -LG_sc(3)+ML3, ML3 )
          weight_uf_k(iktrj,n,s) = abs( zwork0(i1,i2,i3) )**2

          call watchb( tt, ttt(:,5) )

       end do ! n
       end do ! k
       end do ! s

       call watchb( t2 )
       if ( disp_switch_parallel ) then
          write(*,'(1x,"iktrj/nktrj_0=",i4," /",i4,4x,"time=",2f10.3)') &
          iktrj,nktrj_0,t2-t1
          do i=1,5
             write(*,'(1x,"(time",i1,")=",2f10.5)') i,ttt(1,i),ttt(2,i)
          end do
       end if

    end do ! iktrj

! ---

    call rsdft_allreduce_sum( weight_uf, comm_spin )
    call rsdft_allreduce_sum( weight_uf, comm_bzsm )
    call rsdft_allreduce_sum( weight_uf, comm_band )
    call rsdft_allreduce_sum( iktrj_2_k, comm_bzsm )

    call rsdft_allreduce_sum( weight_uf_k, comm_spin )
    call rsdft_allreduce_sum( weight_uf_k, comm_bzsm )
    call rsdft_allreduce_sum( weight_uf_k, comm_band )

    if ( iswitch_ufld_l ) then
       do j3=0,ML3-1
          call rsdft_allreduce_sum( weight_uf_z(:,:,:,j3), comm_spin )
          call rsdft_allreduce_sum( weight_uf_z(:,:,:,j3), comm_bzsm )
          call rsdft_allreduce_sum( weight_uf_z(:,:,:,j3), comm_band )
       end do
    end if

! --- write unfolding data ---

    do iktrj=1,nktrj_0

       k = iktrj_2_k(iktrj)
       if ( k == 0 ) cycle

       if ( myrank == 0 ) then
          write(unit_uf,'(1x,"iktrj=",i6,2x,2(3f10.6,2x))') &
               iktrj,kxyz_pc(1:3,iktrj),kxyz_0(1:3,iktrj)
          do n=1,MB
             write(unit_uf,'(1x,2i6,2(f20.15,2x,g20.10,2x))') &
                  k,n,( esp(n,k,s),weight_uf(iktrj,n,s),s=1,MS )
          end do
          write(unit_uf_k,'(1x,"iktrj=",i6,2x,2(3f10.6,2x))') &
               iktrj,kxyz_pc(1:3,iktrj),kxyz_0(1:3,iktrj)
          do n=1,MB
             write(unit_uf_k,'(1x,2i6,2(f20.15,2x,g20.10,2x))') &
                  k,n,( esp(n,k,s),weight_uf_k(iktrj,n,s),s=1,MS )
          end do
       end if

       if ( iswitch_ufld_l .and. myrank == 0 ) then
          write(unit_uf_z,'(1x,"iktrj=",i6,2x,2(3f10.6,2x))') &
               iktrj,kxyz_pc(1:3,iktrj),kxyz_0(1:3,iktrj)
          do i3=0,Ngrid(3)-1
          do n=1,MB
             write(unit_uf_z,'(1x,3i6,2(f20.15,2x,g20.10,2x))') &
                  k,n,i3,( esp(n,k,s),weight_uf_z(iktrj,n,s,i3),s=1,MS )
          end do
          end do
       end if

    end do ! iktrj

! ---

    if ( iswitch_ufld_l ) then
       deallocate( zwork3 )
       deallocate( zwork2 )
    end if
    deallocate( iktrj_2_k )
    deallocate( zwork1    )
    deallocate( zwork0    )

    call finalize_fft

    call write_border( 0, " band_unfold(end)" )

  END SUBROUTINE band_unfold

  SUBROUTINE finalize_band_unfold
    if ( iswitch_banduf .and. myrank == 0 ) close(unit_uf)
    if ( iswitch_banduf .and. myrank == 0 ) close(unit_uf_k)
    if ( iswitch_ufld_l .and. myrank == 0 ) close(unit_uf_z)
  END SUBROUTINE finalize_band_unfold

END MODULE band_unfold_module


text:
      SUBROUTINE slCMBN ( NSEL, NCAND, LIST, J )
*+
*     - - - - - -
*      C O M B N
*     - - - - - -
*
*  Generate the next combination, a subset of a specified size chosen
*  from a specified number of items.
*
*  Given:
*     NSEL     i        number of items (subset size)
*     NCAND    i        number of candidates (set size)
*
*  Given and returned:
*     LIST     i(NSEL)  latest combination, LIST(1)=0 to initialize
*
*  Returned:
*     J        i        status: -1 = illegal NSEL or NCAND
*                                0 = OK
*                               +1 = no more combinations available
*
*  Notes:
*
*  1) NSEL and NCAND must both be at least 1, and NSEL must be less
*     than or equal to NCAND.
*
*  2) This routine returns, in the LIST array, a subset of NSEL integers
*     chosen from the range 1 to NCAND inclusive, in ascending order.
*     Before calling the routine for the first time, the caller must set
*     the first element of the LIST array to zero (any value less than 1
*     will do) to cause initialization.
*
*  2) The first combination to be generated is:
*
*        LIST(1)=1, LIST(2)=2, ..., LIST(NSEL)=NSEL
*
*     This is also the combination returned for the "finished" (J=1)
*     case.
*
*     The final permutation to be generated is:
*
*        LIST(1)=NCAND, LIST(2)=NCAND-1, ..., LIST(NSEL)=NCAND-NSEL+1
*
*  3) If the "finished" (J=1) status is ignored, the routine
*     continues to deliver combinations, the pattern repeating
*     every NCAND!/(NSEL!*(NCAND-NSEL)!) calls.
*
*  4) The algorithm is by R.F.Warren-Smith (private communication).
*
*  P.T.Wallace   Starlink   25 August 1999
*
*  Copyright (C) 1999 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      INTEGER NSEL,NCAND,LIST(NSEL),J

      INTEGER I,LISTI,NMAX,M
      LOGICAL MORE

*  Validate, and set status.
      IF (NSEL.LT.1.OR.NCAND.LT.1.OR.NSEL.GT.NCAND) THEN
         J = -1
         GO TO 9999
      ELSE
         J = 0
      END IF

*  Just starting?
      IF (LIST(1).LT.1) THEN

*     Yes: return 1,2,3...
         DO I=1,NSEL
            LIST(I) = I
         END DO

      ELSE

*     No: find the first selection that we can increment.

*     Start with the first list item.
         I = 1

*     Loop.
         MORE = .TRUE.
         DO WHILE (MORE)

*        Current list item.
            LISTI = LIST(I)

*        Is this the final list item?
            IF (I.GE.NSEL) THEN

*           Yes:  comparison value is number of candidates plus one.
               NMAX = NCAND+1
            ELSE

*           No:  comparison value is next list item.
               NMAX = LIST(I+1)
            END IF

*        Can the current item be incremented?
            IF (NMAX-LISTI.GT.1) THEN

*           Yes:  increment it.
               LIST(I) = LISTI+1

*           Reinitialize the preceding items.
               DO M=1,I-1
                  LIST(M) = M
               END DO

*           Break.
               MORE = .FALSE.
            ELSE

*           Can't increment the current item:  is it the final one?
               IF (I.GE.NSEL) THEN

*              Yes:  set the status.
                  J = 1

*              Restart the sequence.
                  DO I=1,NSEL
                     LIST(I) = I
                  END DO

*              Break.
                  MORE = .FALSE.
               ELSE

*              No:  next list item.
                  I = I+1
               END IF
            END IF
         END DO
      END IF
 9999 CONTINUE

      END

ext: f

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
subroutine sup_reanut(a,itask)
!-----------------------------------------------------------------------
!****f* Nstinc/nsi_reanut
! NAME
!    nsi_reanut
! DESCRIPTION
!    This routine reads the numerical treatment for the incompressible NS
!    equations.
! USES
!    a%Listener%listen
! USED BY
!    nsi_turnon
!***
!-----------------------------------------------------------------------
   use typre
   use Mod_Listen
   use Mod_ThreeField
   implicit none

   class (ThreeFieldNSProblem) :: a
   integer(ip) :: itask

   integer(ip) :: istab,ipart,rpalo,auxit

   real(rp)  :: auxre,auxre2
   !For Listener
   real(rp), pointer     :: param(:)
   character(5), pointer :: words(:)
   integer(ip), pointer  :: nnpar,nnwor

   call a%Listener%getarrs(words,param,nnpar,nnwor)

   if (itask == 0) then                                  ! Initializations (defaults)

      a%staco(1)  = 4.0_rp                               ! Viscous part of tau1
      a%staco(2)  = 2.0_rp                               ! Convective part of tau1
      a%staco(3)  = 0.25_rp                              ! Rotation part of tau1
      a%staco(4)  = 0.25_rp                              ! Coefficient of tau2
      a%kfl_repro = 0                                    ! Res. projection not used
      a%kfl_reproBoundzero=1
      a%kfl_penal = 0                                    ! Penalty form
      a%kfl_repro_SkipFE = 0                             ! Default is do not skip the Finite Element Part
      a%kfl_shock = 0                                    ! Shock capturing off
      a%shock     = 0.0_rp                               ! SC parameter
      a%kfl_wtemp = 1                                    ! dT/dt is in the residual
      a%kfl_wlapl = 1                                    ! nu Lapl(u) is in the residual
      a%kfl_tiacc = 1                                    ! First order time integ.

      a%kfl_hdifumin=0
      a%kfl_tausm=0                                      ! Smooth taus
      a%kfl_ntausmooth=1                                 ! Number of consecutive taus smooths

      a%kfl_trasg = 0                                    ! Tracking of subscales
      a%mtrit     = 1                                    ! Tracking iterations
      a%kfl_tacsg = 0                                    ! Time accuracy of subscales
      a%kfl_nolsg = 0                                    ! Non-linear subscales
      a%tosgs = 1.0_rp                                   ! Tracking tolerance
      a%relsg = 1.0_rp                                   ! Tracking relaxation
      a%kfl_nolsgNewtonRaphson = 0                       ! Non-linear subscales, newton-Raphson scheme
      a%kfl_stabm = 1                                    ! Default is Variational Multiscale
      a%subrelax  = 1.0_rp                               ! Relaxation Parameter
      a%incremental = 60_ip                              ! Incremental parameter
      a%penal = 6
      a%kfl_splitOSSMomentum = 1                         ! Split-OSS in momentum equation. 1=Off, 0=On
      a%kfl_splitOSSConstitutive = 1                     ! Split-OSS in constitutive equation. 1=Off, 0=On
      a%kfl_linearConstitutiveTerms = 1                  ! Linearization terms in constitutive equation. 0=Picard, 1=Newton
      a%kfl_linearConvectiveTerm = 1                     ! Linearization terms of convective term.  0=Picard, 1=Newton
      a%kfl_reproTemporalTermZero=1                      ! Temporal term in residual. 1=Off, 0=On

   elseif (itask == 1) then     !Inside the Numerical Treatment Block
      if(words(1)=='STABI') then
        do istab = 1,4
           a%staco(istab) = param(istab)
        end do

      else if(words(1).eq.'HDIFU') then
         if(a%Listener%exists('MINIM')) then
            a%kfl_hdifumin = 1
         endif

      else if(words(1).eq.'TAUSM') then
         if(a%Listener%exists('ON   ')) then
             a%kfl_tausm = 1
             a%kfl_ntausmooth=param(2)
         end if

      else if(words(1) == 'PENAL') then
         if(words(2)=='CLASS') a%kfl_penal = 1
            a%penal=param(2)

      else if(words(1)=='TYPEO') then
         if(words(2)=='RESID') then
            a%kfl_repro = 1 !OSS
            if(words(3) == 'SKIPF') then
               a%kfl_repro_SkipFE = 1
            elseif(words(3) == 'NOSKI') then
               a%kfl_repro_SkipFE = 1
            endif
         endif

         if(words(2)=='SPLIT'.and. words(3) == 'OSS')then
            a%kfl_repro = 2 !Split OSS
            a%kfl_splitOSSMomentum = 0
            !a%kfl_splitOSSConstitutive=0

            if (a%MatProp(1)%lawvi>=0) then
               a%kfl_repro=4
               a%kfl_splitOSSConstitutive=0
            end if
         end if

         if(words(2)=='SPLIT' .and. words(3) == 'TOTAL')then
            a%kfl_repro = 3 !Split ASGS
            a%kfl_splitOSSMomentum = 0
         end if

      else if(words(1)=='BOUND' .and. words(2)=='OFF') then
         a%kfl_reproBoundzero=0
          write(*,*) 'Projection on the boundaries is available now.'

      else if(words(1)=='TESTF') then
         if(words(2)=='VMS  ') then
            a%kfl_stabm = 1
         elseif(words(2)=='SUPG ') then
            a%kfl_stabm = 0
         elseif(words(2) == 'GLS  ') then
            a%kfl_stabm = -1
         endif

      else if(words(1)=='SHOCK') then
         if(words(2) == 'ON ') then
            a%kfl_shock = 1
            a%shock(1:2) = param(2:3)
         else if(words(2) == 'GRADU  ') then
            a%kfl_shock = 2
            a%shock(1:2) = param(2:3)
         else if(words(2) == 'DEVSS  ') then
            a%kfl_shock = 3
            a%shock(1:2) = param(2:3)
         else if(words(2) == 'ALL') then
            a%kfl_shock = 4
            a%shock(1:2) = param(2:3)
         end if

      else if(words(1) == 'SUBRE')then
         if(words(2) == 'ON ') a%subrelax = param(2)

      else if(words(1) == 'INCRE')then
         if(words(2) == 'ON ') a%incremental = param(2)

      else if(words(1)=='TEMPO') then
         if(a%Listener%exists('OFF  ')) a%kfl_wtemp = 0

      else if(words(1)=='VISCO') then
         if(a%Listener%exists('OFF  ')) a%kfl_wlapl = 0

      else if(words(1)=='TRACK') then
         if (words(2) == 'ON   ') then
            a%kfl_trasg = 1
         else
            a%kfl_trasg = 0
         end if
      else if(words(1)=='NONLI') then
         if(a%Listener%exists('LINEA')) then
           a%kfl_nolsg=0
         else
           a%kfl_nolsg=1
           a%mtrit = a%Listener%getint('ITERA',20,'#Tracking iterations')
           a%tosgs = a%Listener%getrea('TOLER',1.0e-8_rp,'#Tracking tolerance')
           if (a%Listener%exists('NEWTO')) a%kfl_nolsgNewtonRaphson = 1
         end if
      else if(words(1)=='RELAX') then
        a%relsg = a%Listener%getrea('RELAX',1.0_rp,'#Tracking relaxation')

      else if(words(1)=='ACCUR') then
         if(a%Listener%exists('DYNAM')) then
            a%kfl_tacsg=1
            if (a%kfl_repro==1 .or. a%kfl_repro==2) a%kfl_reproTemporalTermZero=0 !Switch off temporary term at the residual.
         end if

         if(a%Listener%exists('QUASI')) a%kfl_tacsg=0

      else if(words(1)=='LINEA') then
         if(words(2)=='NEWTO') then
            a%kfl_linearConvectiveTerm=1
            a%kfl_linearConstitutiveTerms=1
         else
            a%kfl_linearConvectiveTerm=1
            a%kfl_linearConstitutiveTerms=0
         end if

       else if(words(1)=='LINCO') then
         if(words(2)=='NEWTO') then
            a%kfl_linearConstitutiveTerms=1
         else
            a%kfl_linearConstitutiveTerms=0
         end if
      endif

   elseif (itask == 100) then   !Finalize reading operations
      if (a%kfl_tacsg > 0) a%kfl_trasg = 1
      if (a%kfl_nolsg > 0) a%kfl_trasg = 1

      if (a%kfl_advec==0) a%kfl_linearConvectiveTerm=0

   endif

end subroutine sup_reanut

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
C
       real*4 function ct_enh_inv_grey( rgb, ci, cimin, cimax )
C      --------------------------------------------------------
C
C Function defining an enhanced grey-scale.
C
C defining:
C             x = (cimax-ci)/(cimax-cimin)
C then:
C   ct_enh_inv_grey = ((x-lower_bound)/(upper_bound-lower_bound))**power
C
C for:
C  lower_bound  <=  x  <=  upper_bound
C
C otherwise:
C   ct_enh_inv_grey = 0.0
C
C Where the parameters upper_bound, lower_bound and power_law are passed
C in common.
*-

       include '../include/tv_modify.inc'

       integer        rgb, ci, cimin, cimax
       real*4         x
       x  = (float(ci-cimin)/float(cimax-cimin))
       if (x.ge.lower_bound .and. x.le.upper_bound) then
         ct_enh_inv_grey = ((upper_bound-x)/(upper_bound-lower_bound))**
     *                 power_law
       else if (x.lt.lower_bound) then
         ct_enh_inv_grey = 1.0
       else
         ct_enh_inv_grey = 0.0
       end if
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
      SUBROUTINE DROUT(A,B,AD,BD,NDAMPR)
C
C        'DROUT' CHECKS FOR MAXIMUM MODE SHAPE AND SAVES FOR UPPER
C        LIMIT TO OUTPUT DEFLECTIONS
C
      IMPLICIT REAL * 8 (A-H,O-Z)
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      DIMENSION A(10,3),B(10,3),AD(10,3),BD(10,3)
C
C
      N=NDAMPR
      M=0
      DO 10 I=1,N
   10 M=MAX0(NK(I),M)
C
      CALL AVAL('DIN     ',3,A,10,N,M,2)
C
      CALL AVAL('DINDOT  ',6,AD,10,N,M,2)
C
      CALL AVAL('DOUT    ',4,B,10,N,M,2)
C
      CALL AVAL('DOUTDT  ',6,BD,10,N,M,2)
C
C
      RETURN
      END

ext: for