text:
Spencer Higgins was the unsuccessful Student Focus Candidate for ASUCD President in the Winter 2007 ASUCD Election, losing by nearly 2,000 votes (but getting more than double Rob Dont Vote for Me Roys total). He will be interning for the Bush Administration during the Fall quarter of 2007.

Spencer was an Academic Affairs Commission member, Theta Xi Philanthropy Chair, and an ASUCD intern before being elected to ASUCD Senate in the Fall 2005 ASUCD Election from the Student Focus slate. Higgins was also an intern for Governor Schwarzenegger in 2006.  Higgins lost in his Independent Senate bid in the Fall of 2004 in the 6th round.

Spencer http://www.californiaaggie.com/media/storage/paper981/news/2006/08/07/CampusNews/Yearbook.Support.Group.Takes.Shape.On.Facebook2143670.shtml?norewrite200612130940&sourcedomainwww.californiaaggie.com advised Mitty Chang to form a Yearbook http://ucdavis.facebook.com/group.php?gid2204741944 Facebook to gain support for the project.

 June 2008 City Council

A group of students are trying to draft Spencer to run for City Council in 2008. It is currently a closed group. Image(spencergroup.jpg, thumbnail, right, 149) The draft is chaired by Users/WencheMolenaar Wenche Molenaar for the month of July. The Facebook can be found http://ucdavis.facebook.com/group.php?gid2426414607 here.

The group was unsuccessful in its attempts to persuade Spencer Higgins to run for City Council, he told them no.

A Haiku About Spencer

The Spencer Higgins
Warrior of light and wind
Show us happiness

J. Bleckman

Coffee House Controversy

    http://www.youtube.com/watch?v3juDvo2HbHk A video where he and Tiny Sanders argue that increases to both coffee house wages and to ASUCD stipends will only cost ASUCD $10,000. While former ASUCD Controller and Student Focus Campaign manager Kai SavareeRuess stood by these figures, Paul Harms, later appointed ASUCD Controller, and CoHo management proved these figures to be tens of thousands of dollars off; please see here for more of Paul Harms info: File(HarmsAnalysis.xls). Notes on this analysis by Paul Harms are also provided below in the comments section.

  The employee discounts proposed could reduce income or raise costs (depending on how they would be budgeted) by at least $5,000, assuming that 1,000 ASUCD employees would save $5 per year through employee discounts.

ASUCD Presidential Forum Videos

   Video of/about ?
   Video of/about ?
   Spencer commits all Focus candidates elected (Molly and Cem) to personally outreach to different club every week
  Video of/about ?
 (videos deleted by Steve Ostrowski)

 Winter 2007 Candidate Statement

  My Fellow Students,

  ASUCD needs an overhaul. For too long the students have been neglected, the general population underrepresented and too long has our money been misspent and perverted. Our budget is sporadic and inefficient; it needs to be taken in an entirely new direction. More money needs to be distributed to campus clubs, organizations and students’ hands.

  As ASUCD President, I will give more money back to students through services, entertainment and cash. I propose wage increases, discounts for student employees and more money for Club Finance Council. I will work to reinvigorate on campus entertainment and will work for more student outreach and advocacy.

  I want to do these things because we deserve it. We, as students of UCDavis, deserve a student government that gives back. We deserve to leave this school knowing that these were the best years of our lives. We are proud, We are Aggies!

  On February 21, vote for the party that will make you be proud to be an Aggie, Vote Student Focus.

    It’s Time for Change,
      Spencer Higgins
      ASUCD President

 Fall 2005 Candidate Statement

  Whats up everybody?!!! My name is SPENCER HIGGINS and Im running for Senate on the STUDENT FOCUS slate. Im a second year double major in International Relations and Economics. My on campus experience consists of:

  ASUCD Academic Affairs Commission
  Theta Xi Philanthropy Chair
  ASUCD Senate Intern

  Those are the reasons why I am qualified to run for ASUCD Senate.  These are some of the things I want to do:

  More fun on campus!
    Concerts and Movie premieres!
    Keeping Aggie Pack Alive
  Better Savings
    Roll over Pages  where your left over campus printing pages roll over until the end of each quarter!
  Expanding UCD
    Creation of a Publicity Committee
    Creation of a Business Major

  TOGETHER WE CAN GET THINGS DONE!  So dont forget to vote November 16th and 17th at http://elections.ucdavis.edu elections.ucdavis.edu for SPENCER HIGGINS and the rest of the wonderful STUDENT FOCUS slate.

Campaign Promises and Record as Senator

    Concerts and Movie premieres: Not achieved (campus cinema still defunct)
    Keeping Aggie Pack Alive: Achieved status quo funding as opposed to cuts, by managing to save them from being purged during budget hearings
    Roll over Pages: Not achieved
    Creation of a Publicity Committee: Not achieved
    Creation of a Business Major: Not achieved
    The bill, 20062007 ASUCD Senate Bill #9, to give ASUCD funding to host a paint ball competition as a Theta Xi philanthropy event, this bill has vanished from all records, it was  withdrawn before consideration.
    The only other bill he wrote as a Senator was to fund http://asucd.ucdavis.edu/government/govdocuments/senatebills/20052006/sb46winter06.doc portapotties by the Marya Welch Tennis Center.

I didnt advise mitty to create a yearbook, i just helped him out with information regarding our previous yearbook.  His goal was to bring back a yearbook. SpencerHiggins

Former intern/assistant to two ASUCD Controllers and fourtime LEAD campaign treasurer Paul Harms calculates the cost to be $59,400 ($59,433.06 to be more exact, calculated with a $0.25/hr raise for all CoHo employees and a raise of one pay scale for all 100 ASUCD nonaggie employees on stipend pay, excluding the external secretary, a position likely to be cut. This includes employee benefits of 2.79% and assumes that working hours budgeted at the CoHo remain constant at 124,349 per year from FY20062007.)  The Coffee House management also independently verified Pauls calculations for the proposed CoHo wage increase, which alone would cost $31,954.58.  The detailed calculations, in MS Excel Spreadsheets, are available on request; please contact Paul at prharms@ucdavis.edu.

20070219 10:19:46 nbsp There are, of course, two ways they could make their wage/stipend increase proposal cost only $10,000.  Unfortunately, both ways are ridiculous.  First, the raises could be so minimal and insignificant (to the tune of an extra dollar per week) that, added together, they only cost $10,000.  The second and really scary possibility is that they plan to invent $50,000 in income to offset the $60,000 cost of wage increases.  Well, actually, theyd have to invent $80,000 in income to do this, since we are guaranteed to lose $30,000 in income in the next budget.  Such an invention (what I like to call BUDGET LIES, used extensively by the last Student Focus administration...just look at the FY20042005 ASUCD Creative Media budget) would jeopardize the financial health of ASUCD by drawing money out of our reserves to pay for operational expenses and give a minimal benefit to only a few hundred students, certainly not benefiting the 24,000 students ASUCD represents this year. Users/PaulHarms

Steve, no one was opposing Spencer in trying to create a Business Major. No one was opposing Spencer in increasing roll over pages. No one was opposing Spencer from bringing more concerts. Commissioners, Senators, and Unit members would have welcomed him and the legislation. I would agree that the publicity committee idea, if brought to a commission or Senate would have failed.

i would like to point out that.... I fought for all of those changes and was quickly told no chance by the administration on the Business Major and Roll Over Pages then was in fact COMPLETELY OPPOSED from working with Entertainment Council We only work with LEAD was the response. completely ridiculous.  And a publicity unit is still a great idea, but good luck passing things in a Senate where every idea i brought to the table was instantly opposed by james schwab and therefore the majority of the LEAD Slate.

and Student Focus and I are the ONLY reason aggie pack stayed at the minimal funding it had.... as Darnell proposed it be almost entirely cut... we chose our battles in Budget Hearings because we knew that the LEAD majority would only allow us to argue sooo long until they would just close the budget.  We were such a minority that LEAD had unlimited opportuity to make GREAT changes for the Association... and failed miserably.... take a look at COOO (LEAD Baby), it failed miserably, TWICE... and LEAD attempts to cover it with a larger institution (outreach assembly)... smart move.

look, if you want real change...and i mean CHANGE, then vote for Tiny and I.  We will actually end this ignorant bickering and make changes that will be seen by the students.  Users/SpencerHiggins
  As the primary author on the Outreach Assembly thingy, Id like to say that it will probably save us about $7,000.  It should have a small programmatic budget, and pay the chair $49 a week.  Thats it.  It still saves us a lot over the COOO (which everyone agrees was a LEAD failure).  I dont see how everyone keeps saying its going to cost a lot.   Users/BrentLaabs
  You should have known better that ASUCD has nothing to do with IET. Also, you should have known that the academic senate runs the show in regards to what majors are offered at UCD. Your campaign promises are a lot like those made by the 5th grade student council candidate who promises that hell abolish homework and have the cafeteria serve fast food every day. The office you ran for NEVER would put you in a position to do those things and you should have known that.  Users/WilliamLewis
  You are proving the very valuable point that ASUCD resolutions amount to squat when it comes to Academic policy.  The fact is that being an ASUCD Senators grants some attention ot a particular problem more than any one student could.  Users/SteveOstrowski
  Yep the COOO structure failed. Thats why were trying a new one for cheaper! LOL, dude you never brought ideas to the table, you only authored one piece of legislation. Furthermore http://asucd.ucdavis.edu/government/govdocuments/budget/fy200607/generalprograms.pdf Check this out: The budget shows that Darnell proposed to maintain the same funding for Aggie pack as it had been the year before, total of $6,000. However the approved was $14,000. I would say that is a huge increase, who was responsibile for the increase, Spencer you most likely pushed for it, and then it was passed by a LEAD majority. Entertainment is not opposed to working with Student Focus, the Director and one of your current Senators (Tung and Frick) are really good friends. Entertainment Council, like most people is opposed to working with creepers.The administration always says no at first. They told Rob no, when requested staplers. They told Darnell no when he requested a fix route tipsy. But guess what they worked hard, and were determined to succeed and they did. I dont think people want a leader who quickly gives up.So in summary, you were wrong about budget hearings, wrong about Entertainment Council, and admitted to failing to work hard for your goals. But you were right about the COOO, which everyone already agrees on. Good luck becoming Pres. Users/JamesSchwab
  The presidents proposed budget shows the amount Aggie Pack was cut to before the LEAD administration headed by Caliph wrote the Money From Nowhere Act, that took money from Capital Reserves to fund the budget, in order to increase Aggie Pack to its previous years level of $14,000. And Spencer and I had to then sit and fight at budget hearings and head up negotiations just to get Aggie Pack back to the level they have received for the last 3+ years because people couldnt understand that LEAD the previous year took money that didnt really exist and put it into the operating budget and that is why it wasnt reflected in the final budget numbers for that year. $14,000 wasnt a huge increase, it was simply the same amount of money that had received before. And furthermore, with LEAD touting their fiscal responsibility and socalled balanced budgets, this tactic used two years ago by a LEAD administration is essentially a budget deficit. You take money from the reserves to fund an operating budget, thats a deficit. Just because you want to make everyone happy at budget hearings, instead of doing the responsible thing and passing a budget with the money youre given doesnt make it right. So, in summary, you were wrong about budget hearings. Users/ChrisHerold
       Actually, that is not true. From 04/05 to 05/6 Aggie Pack was reduced from $18k to $6. After that budget was passed, Caliph and Darnell each took money of of their budgets to make up for some of the loss. Then for 06/07 Darnell and Kai proposed to keep it at $6,000. The Senate took it back to $14k with the extra money Darnell and Kai left them to work with in the budget (though Natalia Farhadmotamed wanted to also cut money from Scholarships to increase Aggie Pack funding). Furthermore, the original deduction form 04/05 to 05/06 was passed under a nonLEAD Senate, while the increase in the 06/07 was passed by a LEAD dominated Senate. James Schwab
        Admittedly, the Money from Nowhere Act was kind of shady.  Essentially, youre both right about the increase.  But really, this stuff is not very much out of the norm for ASUCD budgeting.  We have pieinthesky income line items that never materialize, and unit directors routinely run over.  We only noticed months later that the Aggie Student Store was selling items under cost.  Now this happened with a LEAD President and a Student Focus Controller, so I think that theres plenty of blame to go around.  And every time we do this kind of voodoo budgeting, its like cutting money from the Coffee House (and, incidentally, it makes Paul Harms eat a baby).  ASUCD has had unbalanced budgets since before LEAD and Focus existed, so theres plenty of blame to go around.  So instead of blaming each other for something thats essentially everybodys fault, could we all work together next year to get a decent budget?   Users/BrentLaabs

20070221 19:49:27 nbsp working with Entertainment Council We only work with LEAD was the responseSpencer Higgins

Just for the record, that is completely false. Heres what happened which may have triggered Spencers statement on EC: Spencer approached me in mid August (?) to host a benefit concert for Cal Aggie Kids camp. He expressed interest in bringing Jurassic 5 and having them at the Mondavi center. Needless to say, it was a bit much for me to take in. Not only had James JUST approach me asking me for help with the ASUCD Voter Registration event, but I had also just become acclamated with my position as the director of the EC (having no previous experience with it in the past). A idea like that was out of the questionespecially with the limited funding the EC gets. After I turned down Spencers idea (not Student Focus, but Spencer), he asked for my assistance in helping him PLAN the event and teach him how to go about putting on a show. I agreed to help him with whatever he needed, but soon after, I never heard from him again. Never once in our conversation did I say we only work with LEAD nor did I ever give off that image.

Im not a big fan of bipartisanships to begin with, but that doesnt mean I go out of my way to favor one over the other. Yes, I will admit that many LEAD candidates have approached me and been helpful with the things Ive tried to do, but the same could be said for many of the Student Focus representatives. Like James mentioned above, my adoptaunit senator is Alex Frick (Student Focus) and Eric Friedman and Mara Harris have been avid supporters of what Im trying to do. In closing, I hope any Student Focus or LEAD candidate will understand that my opinions and views do not express that of the Entertainment Council and its constituents. EC prides itself in staying neutral, as any ASUCD unit should. Fasho.
Users/EmilyTung

20070221 20:29:52 nbsp I dont think Spencer gets enough credit for voting no for the Keak da Sneak voter registration concert, which was a huge failure.   Users/GregWebb

  Wait, Spencer voted no on a bill that brought a concert to campus!  Wasnt bringing concerts to campus part of his platform when he ran for Senate (see above)?  Why did he vote no? Users/PaulHarms
  Oh yeah, and PS how was it a failure? I saw more people at that event than any other ASUCD event Ive been to before!  I was there...where were you, Greg Webb?  Users/PaulHarms
  I cant speak for Spencer but it was a large allocation of funds for what couldve been a more expanded crowd.  The purpose was voter registration and the concert failed to do that, how many people actually registered to vote? .  Users/GregWebb
  I would be careful making a statement like that and if youre going to make a statement like that, make sure to have concrete evidence as to why you think so. And the Voter Registration event was to increase statewide voter registration among college students (not to be confused with ASUCD voting because ASUCD elections dont require a separate registration). Unfortunately, our government failed at every capacity to increase the voter turnout in the last gubernatorial elections. Unfortunately, shit like that is out of ASUCDs control.   Users/EmilyTung
  And Greg, I would advise you to NOT insult the accomplishments of any ASUCD unitsespecially during YOUR elections. Criticisms about the event are one thing (I think the technical aspects of the event couldve been orchestrated much better), but saying that it was a huge failure is another.   Users/EmilyTung
  I apologize for the language I have used, it was not need.  I am optimistic about what the EC, and ASUCD can do in terms of voter registration and Ill be sure not to miss that chance to help when 2008 rolls around.  Users/GregWebb
  While shit like that is out of ASUCDs control they did/do(?) have a voter registration coordinator, and to be frank there was little to no registration happening at the event. I, however, am not blaming the coordinator for that, but rather the nature and circumstances of the event. Fortunately there was such a show of force from CalPIRG and other sources at other times, but because that event was about registering to vote, and it was about the Green/Lib debate, I personally wouldnt consider it a success for those purposes. For having a good time? Sure. For Bringing Keak Da Sneak to campus? Great! But if we consider the 3 or 5 thousand dollars it cost to put on the event, each of the few voters registered was quite expensive, and the resources could have been used better elsewhere. PS: I was at the event, tabling for my party and with some voter reg. forms, but there coulda been more just voter reg people, and the vast majority of waiting students werent waiting in easily registerable places (like a waaaay too crowded hall). Users/MaxMikalonis
  When you add up all the costs, not just the bill from Senate Reserves, the total is close to $10K. Now, Im all for bringing big names in entertainment to campus, but not to the ARC ballroom where 300 max can get in and ASUCD officials are on the list to an already overcrowded event. That means ASUCD spent over $30/person to let 300 people go see Keak and the other ~21,700 students on this campus had their money spent on something they recieved absolutely no benefit from. However, I dont blame EC for this, they were just doing their job given what they were told. Users/ChrisHerold
        Yup, the concert could have been a lot better. The concert had several purposes, register voters, educate students, hype the election, and entertain. The registration part was lacking but students got educated, entertained and the media covered it from Sacramento to the Bay. Oh ya and the ballroom holds about 800, 300 when there are seats. About 700 students came. I will completely agree that the event could have been so much better, but you learn from your mistakes and make it better. But this is what happens when you actually try and do stuff on campus, and not complaing that people told you no. James Schwab
        ~21,700 students on this campus had their money spent on something they recieved absolutely no benefit from.  Sorta like spending $59,000 to raise ASUCD Stipends and CoHo wages?  :)  Users/PaulHarms

20070724 00:01:09 nbsp Portapotties across the City! Paintball matches in Central Park!  Users/PaulHarms
  Philanthropic Paintball matches... Users/GregWebb
  You are seriously underestimating Higgins interest in City politics.   Users/SteveOstrowski
   I dont underestimate his interest.  Im sure hes quite interested.  Problem for him is, there are lots of people interested in him losing.  Big time. Again.\
     If you are uncomfortable with the situation then you should run yourself or recruit someone to run for City Council.   Users/SteveOstrowski
                  We already have :) Users/JimSchwab

20070724 00:24:04 nbsp This draft Higgins thing is an Ostrowski joke right? Users/GregWebb
                                              No this is not a joke and why is always assumed that I am behind these kinds of things.   Users/SteveOstrowski

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE GDPWND  ( wind, wintyp, iwnclr, iret )
C************************************************************************
C* GDPWND								*
C*									*
C* This subroutine decodes the input for WIND.  The variable has	*
C* several parts separated by a slashes.  The first part contains the	*
C* wind color number.  							*
C* The second part contains the size, width, type of the arrow or barb,	*
C* and the arrowhead size separated by slashes.  The arrow/barb size	*
C* is a multiple of the base size.  Type 1 plots a circle or an		*
C* arrowhead for calm winds.  Type 2 does not plot anything for calm	*
C* winds.  The arrowhead size is a multiple of the base arrowhead size. *
C*									*
C* An example of the wind string is: 6/1.0/5/2				*
C*									*
C* GDPWND  ( WIND, WINTYP, IWNCLR, IRET )				*
C*									*
C* Input parameters:							*
C*	WIND		CHAR*		Wind input			*
C*	WINTYP		CHAR*1		Wind type			*
C*					  B = wind barb 		*
C*					  A = wind arrow		*
C* Output parameters:							*
C*	IWNCLR		INTEGER		Wind color 			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* I. Graffman/RDS	12/84	Original source IP_WIND			*
C* M. desJardins/GSFC	 4/85	Changed default to MKS; changed code.	*
C* I. Graffman/RDS	 5/86	Renamed IN_WIND				*
C* M. desJardins/GSFC	 6/88	Added size				*
C* S. Schotz/GSC	 1/90	Added width				*
C* M. desJardins/GSFC	 5/90	Fixed problem with width but no size	*
C* S. Schotz/GSC	 5/90	Added type				*
C* S. Schotz/GSC	 8/90	Added arrow head size			*
C* K. Brill/NMC		01/93	Return N as output in WINUNI		*
C* D.W.Plummer/NCEP	11/96	From IN_WIND, wintyp input, not output	*
C* D.W.Plummer/NCEP	 7/97	Bug fix stripping color from string	*
C* D.W.Plummer/NCEP	 5/98	Added call to set directional arrows	*
C* D.W.Plummer/NCEP	 2/99	Bug fix when input wind vrbl is blank	*
C************************************************************************
       INCLUDE         'GEMPRM.PRM'
C
	CHARACTER* (*)	wind, wintyp
C*
	CHARACTER	win*24, win2*24, w*1
	REAL		rwind (4)
C------------------------------------------------------------------------
	iret  = 0
C
C*	Check for size.
C
	ibreak = INDEX  ( wind, '/' )
	IF  ( ibreak .eq. 0 )  THEN
	    CALL ST_LSTR ( wind, lens, iret )
	    i = 1
	    w = wind (i:i)
	    DO WHILE ( i .le. lens .and.
     +			.not. (( w .ge. '0' ) .and. ( w .le. '9' )) )
		i = i + 1
	        w = wind (i:i)
	    END DO
	    win = wind (i:)
	    CALL ST_NUMB ( win, iwnclr, ier )
	    IF ( iwnclr .eq. IMISSD )  iwnclr = 1
	    size   = 0.
	    iwidth = 0
   	    itype = 0
            sizehd = 0.
	  ELSE
	    win2 = wind ( : ibreak-1 )
	    CALL ST_LSTR ( win2, lens, iret )
	    i = 1
	    w = win2 (i:i)
	    DO WHILE ( i .le. lens .and.
     +			.not. (( w .ge. '0' ) .and. ( w .le. '9' )) )
		i = i + 1
	        w = win2 (i:i)
	    END DO
	    win = win2 (i:)
	    CALL ST_NUMB ( win, iwnclr, ier )
	    IF ( iwnclr .eq. IMISSD )  iwnclr = 1
	    win2 = wind ( ibreak+1: )
	    CALL ST_RLST  ( win2, '/', 0., 4, rwind, n, ier )
	    size   = rwind (1)
	    iwidth = NINT ( rwind (2) )
	    itype  = NINT ( rwind (3) )
            sizehd = rwind (4)
	END IF
C
C*	Now, set size, width, type, and/or arrow head size
C
	IF  ( ( size .gt. 0. ) .or. ( iwidth .ne. 0 ) .or.
     +        ( itype .ne. 0 ) .or. ( sizehd. gt. 0 ) )  THEN
C	    IF  ( wintyp .eq. 'B' )  THEN
		CALL GSBARB  ( size, iwidth, itype, ier )
C	      ELSE IF  ( wintyp .eq. 'A' )  THEN
		CALL GSARRW  ( size, sizehd, iwidth, itype, ier )
C	      ELSE IF  ( wintyp .eq. 'D' )  THEN
		CALL GSDARR  ( size, sizehd, iwidth, itype, ier )
C	    END IF
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Describe Users/BobBlumenfields here.

20101021 15:36:02 nbsp Welcome to the Wiki, Bob! Users/TomGarberson


text:
	SUBROUTINE TI_DAYW  ( idtarr, idayw, iret )
C************************************************************************
C* TI_DAYW								*
C*									*
C* This subroutine returns the day of the week, IDAYW, given an integer	*
C* time.  IDAYW is set to 1 for Sunday, 2 for Monday, 3 for Tuesday,	*
C* 4 for Wednesday, 5 for Thursday, 6 for Friday, and 7 for Saturday.	*
C*									*
C* This subroutine does not check that a valid time was entered.  The 	*
C* year must be a full four-digit year.					*
C*									*
C* TI_DAYW  ( IDTARR, IDAYW, IRET )					*
C*									*
C* Input parameters:							*
C*	IDTARR (5)	INTEGER		Time array (YYYY,MM,DD,HH,MM)	*
C*									*
C* Output parameters:							*
C*	IDAYW		INTEGER		Day of week (1 - Sun,...)	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	12/87	Adapted from DAYOFWEEK in TRAFIC	*
C* M. desJardins/GSFC	 9/88	Add 1900 to YY				*
C* I. Durham/GSC	 9/98   Added March leap year correction	*
C* I. Durham/GSC	 9/98	Rewrote using new computation		*
C* S. Jacobs/NCEP	11/99	Simplified calculation			*
C************************************************************************
	INTEGER		idtarr (*)
C*
	INTEGER		jdtarr (5)
C------------------------------------------------------------------------
	iret = 0
C
C*	Get the number of days for the first of the given month.
C
	DO  i = 1, 5
	    jdtarr (i) = idtarr (i)
	END DO
	jdtarr (3) = 1
	CALL TI_ITOJ ( jdtarr, jyear, jday, ier )
C
C*	Compute the offset for this month.
C
	ioff  = MOD ( jday-1, 7 )
C
C*	Compute the day of the week.
C
	iy    = idtarr (1)
	id    = idtarr (3)
	idayw = MOD ( ( iy + (iy-1)/4 - (iy-1)/100 + (iy-1)/400 +
     +		        ioff + (id-1) ), 7 ) + 1
C
	RETURN
	END


text:
+incdir+${OC_WB_IP}/ve/wb_uart/tb
+incdir+${OC_WB_IP}/ve/wb_uart/tests

// ${OC_WB_IP}/ve/wb_uart/tb/wb_uart_env_pkg.sv

// ${OC_WB_IP}/ve/wb_uart/tests/wb_uart_tests_pkg.sv

-f ${STD_PROTOCOL_IF}/rtl/wb/wb.f
-f ${OC_WB_IP}/rtl/wb_uart/rtl.f
-f ${SV_BFMS}/wb/wb.f
-f ${SV_BFMS}/event/event.f
-f ${SV_BFMS}/uart/sv.f

${OC_WB_IP}/ve/wb_uart/tb/wb_uart_tb_hdl.sv


text:
c gmm01f
c  THIS FILE IS PART OF THE ORIGINAL GMM01F DISTRIBUTION
c
c  This subroutine is the implementation of Mackowski's three-step
c  method for decomposition of translation matrix into rotational
c  and axial translational parts.
c  The rotation-translation-rotation formulation is described
c  in Mackowski, Proc. R. Soc. Lond. A 433, pp.611-612 (1991).
c  Yu-lin Xu   12/2000
c
      subroutine rtr(anpt,nodrj,nodri,ekt,drot)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter(nmp=np*(np+2))
      double precision drot(*)
      complex*16 anpt(2,*),ant(2,2*np),amt(2,-np:np),a,b,
     +           ekt(*),ek(-np:np),atr(2,np,nmp)
      common/tran/atr
c
      ek(0)=1.d0
      nmax=max(nodrj,nodri)
      do m=1,nmax
         ek(m)=ekt(m)
         ek(-m)=dconjg(ek(m))
      enddo
      irc=0
      do n=1,nodrj
         n1=n*(n+1)
         do m=-n,n
            amt(1,m)=0.d0
            amt(2,m)=0.d0
         enddo
         do k=-n,n
            kn=n1+k
            a=ek(k)*anpt(1,kn)
            b=ek(k)*anpt(2,kn)
            do m=-n,n
               irc=irc+1
               amt(1,m)=amt(1,m)+a*drot(irc)
               amt(2,m)=amt(2,m)+b*drot(irc)
            enddo
         enddo
         do m=-n,n
            imn=n1+m
            anpt(1,imn)=amt(1,m)
            anpt(2,imn)=amt(2,m)
         enddo
      enddo
      mmax=min(nodrj,nodri)
      do m=-mmax,mmax
         n1=max(1,iabs(m))
         do n=n1,nodrj
            imn=n*(n+1)+m
            do ip=1,2
               ant(ip,n)=anpt(ip,imn)
            enddo
         enddo
         do n=n1,nodri
            imn=n*(n+1)+m
            a=0.d0
            b=0.d0
            do l=n1,nodrj
               ml=l*(l+1)+m
               a=a+atr(1,n,ml)*ant(1,l)
     1            +atr(2,n,ml)*ant(2,l)
               b=b+atr(1,n,ml)*ant(2,l)
     1            +atr(2,n,ml)*ant(1,l)
            enddo
            anpt(1,imn) = a
            anpt(2,imn) = b
         enddo
      enddo
      in=1
      irc=0
      do n=1,nodri
         in=-in
         n1=n*(n+1)
         do m=-n,n
            amt(1,m)=0.d0
            amt(2,m)=0.d0
         enddo
         sik=-in
         do k=-n,n
        sik=-sik
        kn=n1+k
        a=sik*anpt(1,kn)
            b=sik*anpt(2,kn)
            do m=-n,n
               irc=irc+1
               amt(1,m)=amt(1,m)+a*drot(irc)
               amt(2,m)=amt(2,m)+b*drot(irc)
            enddo
         enddo
         sik=-in
         do m=-n,n
            sik=-sik
            imn=n1+m
            anpt(1,imn)=amt(1,m)*ek(-m)*sik
            anpt(2,imn)=amt(2,m)*ek(-m)*sik
         enddo
      enddo
      return
      end

max_stars_repo_name: m-ringler/GMM-DIP

max_issues_repo_name: m-ringler/GMM-DIP

max_forks_repo_name: m-ringler/GMM-DIP


text:
program main  ! global and partial sum

  implicit none
  include 'mpif.h'
  integer :: size, rank, ierr
  integer :: sum, i

  call mpi_init(ierr)
  call mpi_comm_size(MPI_COMM_WORLD, size, ierr)
  call mpi_comm_rank(MPI_COMM_WORLD, rank, ierr)

  call mpi_reduce(rank, sum, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr)

  if (rank == 0) write(6,*) "global sum = ", sum

  call mpi_scan(rank, sum, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)

  if (rank /= 0) then
     call mpi_send(sum, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, ierr)
  else
     do i = 0, size - 1
        if (i > 0) then
           call mpi_recv(sum, 1, MPI_INTEGER, i, 0, &
                         MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
        endif
        write(6,*) " partial sum(", i, ") = ", sum
     enddo
  endif

  call mpi_finalize(ierr)

end

max_stars_repo_name: David-Loibl/cprog

max_issues_repo_name: David-Loibl/cprog

max_forks_repo_name: David-Loibl/cprog


MODULE localpot_module

!  use rgrid_module, only: Ngrid, Igrid
  use xc_module, only: Vxc
  use hartree_variables, only: Vh
  use ps_local_module, only: Vion
 !use localpot2_module

  implicit none

  PRIVATE
  PUBLIC :: Vloc, init_localpot, construct_matrix_localpot
  PUBLIC :: op_localpot
  PUBLIC :: update_localpot

  real(8),allocatable :: Vloc(:,:)
  integer :: ML_0,ML_1,MS_0,MS_1
  integer :: MS
  integer,allocatable :: Ngrid(:),Igrid(:,:)

CONTAINS

  SUBROUTINE init_localpot( ML_0_in,ML_1_in, MS_0_in,MS_1_in )
    implicit none
    integer,intent(IN) :: ML_0_in,ML_1_in,MS_0_in,MS_1_in
    call write_border( 0, " init_localpot(start)" )
    ML_0 = ML_0_in
    ML_1 = ML_1_in
    MS_0 = MS_0_in
    MS_1 = MS_1_in
    allocate( Vloc(ML_0:ML_1,MS_0:MS_1) )
    Vloc=0.0d0
    call write_border( 0, " init_localpot(end)" )
  END SUBROUTINE init_localpot

  SUBROUTINE op_localpot( f, Vf, s_in )
    implicit none
#ifdef _DRSDFT_
    real(8),intent(IN) :: f(:,:)
    real(8),intent(INOUT) :: Vf(:,:)
#else
    complex(8),intent(IN) :: f(:,:)
    complex(8),intent(INOUT) :: Vf(:,:)
#endif
    integer,optional,intent(IN) :: s_in
    integer :: s,m,n,i,j

    m = size( f, 1 )
    n = size( f, 2 )
    s = 1 ; if ( present(s_in) ) s=s_in

!    if ( flag_localpot2 ) then
!
!       allocate( ft(Ngrid(0)) )
!       ft(:)=zero
!
!!!$OMP PARALLEL
!       do n=1,nn
!!$OMP SINGLE
!       call mpi_allgatherv(f(1,n),mm,TYP,ft,ir_grid,id_grid,TYP,comm_grid,ierr)
!!$OMP END SINGLE
!!$OMP DO
!          do i=1,mm
!          do j=1,MLpot
!             vf(i,n)=vf(i,n)+vloc_nl(j,i+ML_0-1)*ft(Lpot(j,i+ML_0-1))
!          end do
!          end do
!!$OMP END DO
!       end do
!!!$OMP END PARALLEL
!
!       deallocate( ft )
!
!    else

!!$OMP PARALLEL
    do j=1,n
!$OMP DO
       do i=1,m
          Vf(i,j) = Vf(i,j) + Vloc(ML_0-1+i,s)*f(i,j)
       end do
!$OMP END DO
    end do
!!$OMP END PARALLEL

!    end if

  END SUBROUTINE op_localpot

  SUBROUTINE read_localpot(file_name,rank)
    implicit none
    character(*),intent(IN) :: file_name
    integer,intent(IN) :: rank
    integer,parameter :: unit=80
    integer :: ML_tmp(0:3),ierr,i1,i2,i3,i,s
    integer,allocatable :: LL_tmp(:,:)
    real(8),allocatable :: rtmp(:),rtmp3(:,:,:)
    include 'mpif.h'
    if ( rank == 0 ) then
       open(unit,file=file_name,form='unformatted')
       read(unit) ML_tmp(0:3)
       write(*,*) "ML_TMP=",ML_tmp
    end if
    call mpi_bcast(ML_tmp,4,mpi_integer,0,mpi_comm_world,ierr)
    if ( ML_tmp(0) /= Ngrid(0) ) stop "at read_localpot"
    allocate( LL_tmp(3,ML_tmp(0)) )
    if ( rank == 0 ) then
       read(unit) LL_tmp(:,:)
    end if
    call mpi_bcast(LL_tmp,3*ML_tmp(0),mpi_integer,0,mpi_comm_world,ierr)
    allocate( rtmp3(0:ML_tmp(1)-1,0:ML_tmp(2)-1,0:ML_tmp(3)-1) )
    allocate( rtmp(ML_tmp(0)) )
    do s=1,MS
       if ( rank == 0 ) then
          read(unit) rtmp(:) ! rho
          read(unit) rtmp(:) ! Vloc
       end if
       call mpi_bcast(rtmp,ML_tmp(0),mpi_real8,0,mpi_comm_world,ierr)
       do i=1,ML_tmp(0)
          rtmp3( LL_tmp(1,i),LL_tmp(2,i),LL_tmp(3,i) ) = rtmp(i)
       end do
       if ( MS_0 <= s .and. s <= MS_1 ) then
          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
             Vloc(i,s) = rtmp3(i1,i2,i3)
          end do
          end do
          end do
       end if
    end do ! s
    if ( rank == 0 ) then
       close(unit)
    end if
    deallocate( rtmp   )
    deallocate( rtmp3  )
    deallocate( LL_tmp )
  END SUBROUTINE read_localpot

  SUBROUTINE construct_matrix_localpot( s, ML, Hmat )
    implicit none
    integer,intent(IN) :: s, ML
#ifdef _DRSDFT_
    real(8),intent(INOUT) :: Hmat(ML,ML)
#else
    complex(8),intent(INOUT) :: Hmat(ML,ML)
#endif
    integer :: i

    do i=1,ML
       Hmat(i,i) = Hmat(i,i) + Vloc(i,s)
    end do

  END SUBROUTINE construct_matrix_localpot

  SUBROUTINE update_localpot
    implicit none
    integer :: s
    do s=MS_0,MS_1
       Vloc(:,s) = Vion(:) + Vh(:) + Vxc(:,s)
    end do
  END SUBROUTINE update_localpot

END MODULE localpot_module


text:
module da_qscat

   use module_domain, only : domain

   use da_control, only : obs_qc_pointer,max_ob_levels,missing_r, &
      check_max_iv_print, check_max_iv_unit, v_interp_p, v_interp_h, &
      check_max_iv, missing, max_error_uv, max_error_t, rootproc, &
      qscat, max_error_p,max_error_q, trace_use_dull, &
      max_stheight_diff,missing_data,max_error_bq,max_error_slp, &
      max_error_bt, max_error_buv, anal_type_verify, kms,kme,kts,kte,&
      ob_vars,qcstat_conv_unit, fails_error_max, &
      convert_fd2uv,convert_uv2fd,max_error_spd,max_error_dir,max_omb_spd,max_omb_dir,pi,qc_rej_both,&
      wind_sd_qscat, wind_stats_sd, write_rej_obs_conv
   use da_grid_definitions, only : da_ffdduv, da_ffdduv_model, da_ffdduv_diagnose
   use da_physics, only : da_uv_to_sd_lin, da_uv_to_sd_adj
   use da_define_structures, only : maxmin_type, iv_type, y_type, jo_type, &
      bad_data_type, x_type, number_type, bad_data_type
   use da_interpolation, only : da_to_zk, &
      da_interp_lin_3d,da_interp_lin_3d_adj
   use da_par_util, only : da_proc_stats_combine
   use da_par_util1, only : da_proc_sum_int
   use da_statistics, only : da_stats_calculate
   use da_tools, only : da_max_error_qc, da_residual, da_convert_zk
   use da_tracing, only : da_trace_entry, da_trace_exit

   ! The "stats_qscat_type" is ONLY used locally in da_qscat:

   type residual_qscat1_type
      real          :: u                        ! u-wind.
      real          :: v                        ! v-wind.
   end type residual_qscat1_type

   type maxmin_qscat_stats_type
      type (maxmin_type)         :: u, v
   end type maxmin_qscat_stats_type

   type stats_qscat_type
      type (maxmin_qscat_stats_type)  :: maximum, minimum
      type (residual_qscat1_type)     :: average, rms_err
   end type stats_qscat_type

contains

#include "da_jo_and_grady_qscat.inc"
#include "da_residual_qscat.inc"
#include "da_check_max_iv_qscat.inc"
#include "da_get_innov_vector_qscat.inc"
#include "da_ao_stats_qscat.inc"
#include "da_oi_stats_qscat.inc"
#include "da_print_stats_qscat.inc"
#include "da_transform_xtoy_qscat.inc"
#include "da_transform_xtoy_qscat_adj.inc"
#include "da_calculate_grady_qscat.inc"

end module da_qscat

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


module excel_mo

  use, intrinsic :: iso_c_binding

  implicit none

  type, bind(c) :: datetime_ty

    integer(c_int) :: year
    integer(c_int) :: month
    integer(c_int) :: day
    integer(c_int) :: hour
    integer(c_int) :: min
    real(c_double) :: sec

  end type

  type, bind(c) :: image_options_ty

    integer(c_int)    :: x_offset        = 0
    integer(c_int)    :: y_offset        = 0
    real(c_double)    :: x_scale         = 1.d0
    real(c_double)    :: y_scale         = 1.d0
    integer(c_int)    :: object_position = 0
    integer(c_int)    :: decorative      = 0
    character(c_char) :: description     = c_null_char
    character(c_char) :: url             = c_null_char
    character(c_char) :: tip             = c_null_char

  end type

  interface

    integer(c_int) function name2row ( name ) &
        bind ( c, name = 'name2row_c' )
      import c_char, c_int
      character(c_char), intent(in) :: name
    end function

    integer(c_int) function name2col ( name ) &
        bind ( c, name = 'name2col_c' )
      import c_char, c_int
      character(c_char), intent(in) :: name
    end function

    integer(c_int) function name2row2 ( name ) &
        bind ( c, name = 'name2row2_c' )
      import c_char, c_int
      character(c_char), intent(in) :: name
    end function

    integer(c_int) function name2col2 ( name ) &
        bind ( c, name = 'name2col2_c' )
      import c_char, c_int
      character(c_char), intent(in) :: name
    end function

    type(c_ptr) function workbook_new ( file ) &
        bind ( c, name = 'workbook_new_c' )
      import c_ptr, c_char
      character(c_char), intent(in) :: file
    end function

    type(c_ptr) function workbook_add_worksheet ( workbook, name ) &
        bind ( c, name = 'workbook_add_worksheet_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: workbook
      character(c_char), intent(in) :: name
    end function

    subroutine worksheet_activate ( worksheet ) &
        bind ( c, name = 'worksheet_activate_c' )
      import c_ptr
      type(c_ptr), intent(in) :: worksheet
    end subroutine

    subroutine worksheet_hide ( worksheet ) &
        bind ( c, name = 'worksheet_hide_c' )
      import c_ptr
      type(c_ptr), intent(in) :: worksheet
    end subroutine

    subroutine worksheet_set_first_sheet ( worksheet ) &
        bind ( c, name = 'worksheet_set_first_sheet_c' )
      import c_ptr
      type(c_ptr), intent(in) :: worksheet
    end subroutine

    type(c_ptr) function workbook_add_format ( workbook ) &
        bind ( c, name = 'workbook_add_format_c' )
      import c_ptr
      type(c_ptr), intent(in) :: workbook
    end function

    subroutine workbook_define_name ( workbook, name, formula ) &
        bind ( c, name = 'workbook_define_name_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: workbook
      character(c_char), intent(in) :: name, formula
    end subroutine

    subroutine worksheet_write_datetime( workbook, row, col, datetime, format ) &
        bind ( c, name = 'worksheet_write_datetime_c' )
      import c_ptr, c_int, c_char, datetime_ty
      type(c_ptr),       intent(in)        :: workbook, format
      integer(c_int),    intent(in), value :: row, col
      type(datetime_ty), intent(in)        :: datetime
    end subroutine

    subroutine format_set_shrink ( format ) &
        bind ( c, name = 'format_set_shrink_c' )
      import c_ptr
      type(c_ptr), intent(in) :: format
    end subroutine

    subroutine format_set_bold ( format ) &
        bind ( c, name = 'format_set_bold_c' )
      import c_ptr
      type(c_ptr), intent(in) :: format
    end subroutine

    subroutine format_set_italic ( format ) &
        bind ( c, name = 'format_set_italic_c' )
      import c_ptr
      type(c_ptr), intent(in) :: format
    end subroutine

    subroutine worksheet_autofilter ( worksheet, first_row, first_col, last_row, last_col ) &
        bind ( c, name = 'worksheet_autofilter_c' )
      import c_ptr, c_int
      type(c_ptr), intent(in) :: worksheet
      integer(c_int), intent(in), value :: first_row, first_col, last_row, last_col
    end subroutine

    subroutine worksheet_merge_range ( worksheet, first_row, first_col, last_row, last_col, string, format ) &
        bind ( c, name = 'worksheet_merge_range_c' )
      import c_ptr, c_int, c_char
      type(c_ptr),       intent(in)        :: worksheet, format
      integer(c_int),    intent(in), value :: first_row, first_col, last_row, last_col
      character(c_char), intent(in)        :: string
    end subroutine

    subroutine worksheet_set_row ( worksheet, format, row, height ) &
        bind ( c, name = 'worksheet_set_row_c' )
      import c_ptr, c_int
      type(c_ptr), intent(in) :: worksheet, format
      integer(c_int), intent(in), value :: row
      integer(c_int), intent(in), value :: height
    end subroutine

    subroutine worksheet_set_column ( worksheet, format, col, width ) &
        bind ( c, name = 'worksheet_set_column_c' )
      import c_ptr, c_int
      type(c_ptr), intent(in) :: worksheet, format
      integer(c_int), intent(in), value :: col
      integer(c_int), intent(in), value :: width
    end subroutine

    subroutine worksheet_write_formula ( worksheet, row, col, formula, format ) &
        bind ( c, name = 'worksheet_write_formula_c' )
      import c_ptr, c_int, c_char
      type(c_ptr), intent(in) :: worksheet, format
      integer(c_int),    intent(in), value :: row, col
      character(c_char), intent(in)        :: formula
    end subroutine

    subroutine worksheet_write_number ( worksheet, format, row, col, value ) &
        bind ( c, name = 'worksheet_write_number_c' )
      import c_ptr, c_int, c_double
      type(c_ptr), intent(in) :: worksheet, format
      integer(c_int), intent(in), value :: row, col
      real(c_double), intent(in), value :: value
    end subroutine

    subroutine worksheet_write_string ( worksheet, format, row, col, string ) &
        bind ( c, name = 'worksheet_write_string_c' )
      import c_ptr, c_int, c_char
      type(c_ptr),       intent(in)        :: worksheet, format
      integer(c_int),    intent(in), value :: row, col
      character(c_char), intent(in)        :: string
    end subroutine

    type(image_options_ty) function image_set_options &
        ( x_offset, y_offset, x_scale, y_scale, position, description, url, tip ) &
        bind ( c, name = 'image_set_options_c' )
      import c_int, c_double, c_char, image_options_ty
      integer(c_int),    intent(in), value    :: x_offset, y_offset
      real(c_double),    intent(in), value    :: x_scale, y_scale
      character(c_char), intent(in), optional :: position
      character(c_char), intent(in), optional :: description, url, tip
    end function

    subroutine worksheet_insert_image_opt ( worksheet, row, col, file, options ) &
        bind ( c, name = 'worksheet_insert_image_opt_c' )
      import c_ptr, c_int, c_char, image_options_ty
      type(c_ptr),            intent(in)        :: worksheet
      type(image_options_ty), intent(in)        :: options
      integer(c_int),         intent(in), value :: row, col
      character(c_char),      intent(in)        :: file
    end subroutine

    subroutine worksheet_insert_image ( worksheet, row, col, file ) &
        bind ( c, name = 'worksheet_insert_image_c' )
      import c_ptr, c_int, c_char
      type(c_ptr),       intent(in)        :: worksheet
      integer(c_int),    intent(in), value :: row, col
      character(c_char), intent(in)        :: file
    end subroutine

    subroutine worksheet_write_comment ( worksheet, row, col, string ) &
        bind ( c, name = 'worksheet_write_comment_c' )
      import c_ptr, c_int, c_char
      type(c_ptr),       intent(in)        :: worksheet
      integer(c_int),    intent(in), value :: row, col
      character(c_char), intent(in)        :: string
    end subroutine

    subroutine worksheet_set_header ( worksheet, string ) &
        bind ( c, name = 'worksheet_set_header_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: worksheet
      character(c_char), intent(in) :: string
    end subroutine

    subroutine worksheet_set_footer ( worksheet, string ) &
        bind ( c, name = 'worksheet_set_footer_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: worksheet
      character(c_char), intent(in) :: string
    end subroutine

    subroutine workbook_close ( workbook ) &
        bind ( c, name = 'workbook_close_c' )
      import c_ptr
      type(c_ptr), intent(in) :: workbook
    end subroutine

    ! URL
    subroutine worksheet_write_url ( worksheet, row, col, url, format ) &
        bind ( c, name = 'worksheet_write_url_c' )
      import c_ptr, c_int, c_char
      type(c_ptr), intent(in) :: worksheet, format
      integer(c_int),    intent(in), value :: row, col
      character(c_char), intent(in)        :: url
    end subroutine

    ! Font size
    subroutine format_set_font_size ( format, size ) &
        bind ( c, name = 'format_set_font_size_c' )
      import c_ptr, c_int
      type(c_ptr),    intent(in)        :: format
      integer(c_int), intent(in), value :: size
    end subroutine

    ! Font name
    subroutine format_set_font_name ( format, name ) &
        bind ( c, name = 'format_set_font_name_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: name
    end subroutine

    ! Number format
    subroutine format_set_num_format ( format, style ) &
        bind ( c, name = 'format_set_num_format_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

    ! Tab color
    subroutine worksheet_set_tab_color ( worksheet, color ) &
        bind ( c, name = 'worksheet_set_tab_color_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: worksheet
      character(c_char), intent(in) :: color
    end subroutine

    ! Font color
    subroutine format_set_font_color ( format, color ) &
        bind ( c, name = 'format_set_font_color_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: color
    end subroutine

    ! Underline
    subroutine format_set_underline ( format, style ) &
        bind ( c, name = 'format_set_underline_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

    ! Pattern
    subroutine format_set_pattern ( format, pattern ) &
        bind ( c, name = 'format_set_pattern_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: pattern
    end subroutine

    ! Background color
    subroutine format_set_bg_color ( format, color ) &
        bind ( c, name = 'format_set_bg_color_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: color
    end subroutine

    ! Alignment
    subroutine format_set_align ( format, align ) &
        bind ( c, name = 'format_set_align_c' )
      import c_ptr, c_char
      type(c_ptr), intent(in) :: format
      character(c_char), intent(in) :: align
    end subroutine

    ! Border
    subroutine format_set_border ( format, style ) &
        bind ( c, name = 'format_set_border_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

    subroutine format_set_top ( format, style ) &
        bind ( c, name = 'format_set_top_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

    subroutine format_set_bottom ( format, style ) &
        bind ( c, name = 'format_set_bottom_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

    subroutine format_set_right ( format, style ) &
        bind ( c, name = 'format_set_right_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

    subroutine format_set_left ( format, style ) &
        bind ( c, name = 'format_set_left_c' )
      import c_ptr, c_char
      type(c_ptr),       intent(in) :: format
      character(c_char), intent(in) :: style
    end subroutine

  end interface

  contains

  function cs ( fs ) result ( fsnull )

    character(*), intent(in)  :: fs
    character(:), allocatable :: fsnull

    fsnull = trim(fs)//c_null_char

  end function

end module


text:
! File VersionID:
!   $Id: MacroRate.for 179 2010-05-13 15:08:34Z kroes006 $
! ----------------------------------------------------------------------
      SUBROUTINE MACRORATE(ITask,ICpBtDm,ICpBtPerZon,ICpSatGWl,
     &            ICpSatPeGWl,ICpTpPerZon,ICpTpSatZon,ICpTpWaSrDm,
     &            ArMpSsDm,AwlCorFac,FrMpWalWet,KDCrRlRef,SorpDmCp,
     &            ThtSrpRefDmCp,TimAbsCumDmCp,VlMpDm,VlMpDmCp,WaSrMp,
     &            WaSrMpDm,ZBtDm,ZWaLevDm,         flDraTub,FlEndSrpEvt,
     &            QExcMtxDmCp,QInIntSatDmCp,QInMtxSatDmCp,QInTopLatDm,
     &            QInTopPreDm,QOutDrRapCp,QOutMtxSatDmCp,QOutMtxUnsDmCp)
! ----------------------------------------------------------------------
!     Date               : april 2008
!     Purpose            :
!     Subroutines called : ABSORPTION, RAPIDDRAIN, SATFLOW
!     Functions called   : VOLUNDR
! ----------------------------------------------------------------------
      use Variables
      implicit NONE

! --- global                                                          In
      integer ICpBtDm(MaDm), ICpTpWaSrDm(MaDm), ITask
      real*8  ArMpSsDm(MaDm), AwlCorFac(Macp), FrMpWalWet(MaDm,MaCp)
      real*8  KDCrRlRef(MaDr), SorpDmCp(MaDm,MaCp)
      real*8  ThtSrpRefDmCp(MaDm,MaCp), TimAbsCumDmCp(MaDm,MaCp)
      real*8  VlMpDm(MaDm), VlMpDmCp(MaDm,MaCp), WaSrMp, WaSrMpDm(MaDm)
      real*8  ZBtDm(MaDm), ZWaLevDm(MaDm)
      logical flDraTub(Madr)
!     -                                                              Out
      real*8  QExcMtxDmCp(MaDm,MaCp), QInIntSatDmCp(MaDm,MaCp)
      real*8  QInMtxSatDmCp(MaDm,MaCp), QInTopLatDm(MaDm)
      real*8  QInTopPreDm(MaDm), QOutDrRapCp(MaCp)
      real*8  QOutMtxSatDmCp(MaDm,MaCp), QOutMtxUnsDmCp(MaDm,MaCp)
      logical FlEndSrpEvt(MaDm,MaCp)
! ----------------------------------------------------------------------
! --- local
      integer ic, ICpBtPerZon, ICpSatGWl, ICpSatPeGWl, ICpTpPerZon
      integer ICpTpSatZon, id, jd, Sq(MaDM), SqHlp
      real*8  FrFlwIn, FrFlwOut, FlwInDmExces
      real*8  FlwInDmTot, FlwInIntSatDmCpPot(MaCp), FlwInIntSatDmPot
      real*8  FlwInMtxSatDmCpPot(MaCp), FlwInMtxSatDmPot
      real*8  FlwInTopExcesTot, FlwInTopLatDmPot(MaDm), FlwInTopPot
      real*8  FlwInTopXtr, FlwInTopPreDmPot(MaDm), FlwOutDmExces(MaDm)
      real*8  FlwOutDrRapCpPot(MaCp), FlwOutDrRapPot, FlwOutDmTot(MaDm)
      real*8  FlwOutMtxSatDmCpPot(MaDm,MaCp), FlwOutMtxSatDmPot(MaDm)
      real*8  FlwOutMtxUnsDmCpPot(MaDm,MaCp), FlwOutMtxUnsDmPot(MaDm)
      real*8  FrReduQ, FrXtr, Pp, PpTot, SatDefDm(MaDm)
      real*8  SatDefDmRel(MaDm), SatDefTot, VlMpUndrDrL, VlMpUndrGwL
      real*8  VOLUNDR, WaSrMpDmMax, WaSrMpDmMin, WaSrMpDmTmp

! ----------------------------------------------------------------------
!   for FORCHECK
      t = t
      t1900 = t1900
!
!> I << ---------- INITIALISE AND CHECK MACROPORE STATUS ---------------

! --- Initialisation Fluxes and Derivatives
      QMaPo= 0.d0
      QRapDra= 0.d0
      do 10 id= 1,NumDm
         QInTopPreDm(id)= 0.d0
         QInTopLatDm(id)= 0.d0
         do 9 ic= 1, NumNod
            QExcMtxDmCp(id,ic)   = 0.d0
            QInIntSatDmCp(id,ic) = 0.d0
            QInMtxSatDmCp(id,ic) = 0.d0
            QOutMtxSatDmCp(id,ic)= 0.d0
            QOutMtxUnsDmCp(id,ic)= 0.d0
            if (id.eq.1) then
               QOutDrRapCp(ic) = 0.d0
               QExcMpMtx(ic)   = 0.d0
               dFdhMp(ic)      = 0.0
            endif
   9     continue
  10  continue

! --- Check whether the macropore status requires the working of this subroutine
      FlwInTopPot= QMpLatSs + ArMpSs * (NRaiDt+NIrd+Melt) * DT
      if (VlMp.lt.1.d-6 .or. (WaSrMp.lt.1.d-6 .and. FlwInTopPot.lt.1.d-6
     &    .and. ICpBtDm(1).le.NodGwL .and. NPeGwL.eq.0) ) then ! not required
         do 20 id= 1, NumDm
            do 19 ic= 1, NumNod
               FlEndSrpEvt(id,ic)= .true.
  19        continue
  20     continue
         QMpLatSs= 0.d0
         return
      endif

      FrReduQ= 0.1d0**(float(IDecMpRat))
!      if (IDecMpRat.gt.0) write(94,*) t,IDecMpRat

!> II < ----- CALCULATE RATES AND DERIVATIVES PER MACROPORE DOMAIN -----
!
! --- Subroutine is required; calculations for all macropore domains
!   - General initialisation
      FlwInTopExcesTot= 0.d0
      SatDefTot= 0.d0
!
      do 100 id= 1, NumDm
         FlwInDmTot= 0.d0
         FlwOutDmTot(id)= 0.d0

!- A. CALCULATE POTENTIAL IN- AND OUTGOING WATER FLOWS

! - 1. INGOING WATER FLOWS
!   1.a. INFLoW into macropores at soil surface (FlwInTop...) by infiltration:
!        directly by precipitation (Pre) and indirectly by lateral (Lat)
!        overland flow of matrix infiltration excess (runoff)
         FlwInTopPreDmPot(id)= ArMpSsDm(id) * (NRaiDt+NIrd+Melt) * DT
         FlwInTopLatDmPot(id)= 0.d0
         if (ArMpSs.gt.0.d0) FlwInTopLatDmPot(id)= ArMpSsDm(id)/ArMpSs *
     &                                             QMpLatSs

!   1.b. INFLoW into macropores from SATurated top compartments by INTerflow:
         call SATFLOW(1,ICpBtDm(id),ICpBtPerZon,ICpSatPeGWl,
     &          ICpTpWaSrDm(id),ICpTpPerZon,id,FrMpWalWet,FrReduQ,PeGWL,
     &          QInMtxSatDmCp,QOutMtxSatDmCp,ZWaLevDm(id),
     &          FlwInIntSatDmCpPot,FlwInIntSatDmPot,FlwOutMtxSatDmCpPot,
     &          FlwOutMtxSatDmPot(id))

!   1.c. INFLoW into macropores from SATurated MaTriX compartm. (exfiltration):
         call SATFLOW(2,ICpBtDm(id),NumNod,ICpSatGWl,ICpTpWaSrDm(id),
     &          ICpTpSatZon,id,FrMpWalWet,FrReduQ,GWlFlCpZo,
     &          QInMtxSatDmCp,QOutMtxSatDmCp,ZWaLevDm(id),
     &          FlwInMtxSatDmCpPot,FlwInMtxSatDmPot,FlwOutMtxSatDmCpPot,
     &          FlwOutMtxSatDmPot(id))

! - 2. OUTGOING WATER FLOWS
!   2.a. OUTFLoW out off macrop. into SATurated MaTriX compartm. (infiltrat.):
         call SATFLOW(3,ICpBtDm(id),NumNod,ICpSatGWl,ICpTpWaSrDm(id),
     &          ICpTpPerZon,id,FrMpWalWet,FrReduQ,GWlFlCpZo,
     &          QInMtxSatDmCp,QOutMtxSatDmCp,ZWaLevDm(id),
     &          FlwInMtxSatDmCpPot,FlwInMtxSatDmPot,FlwOutMtxSatDmCpPot,
     &          FlwOutMtxSatDmPot(id))

!   2.b. OUTFLoW out off macrop. into UNSaturated MaTriX compart. (absorption):
         call ABSORPTION(1,ICpBtDm(id),ICpBtPerZon,ICpTpPerZon,
     &          ICpTpSatZon,ICpTpWaSrDm(id),id,AwlCorFac,FrMpWalWet,
     &          FrReduQ,QOutMtxUnsDmCp,SorpDmCp,ThtSrpRefDmCp,
     &          TimAbsCumDmCp,FlwOutMtxUnsDmCpPot,FlwOutMtxUnsDmPot(id),
     &          FlEndSrpEvt,ZWaLevDm)

!   2.c. OUTFLoW out off macropores by RAPid DRainage. Only Main Domain (id= 1):
         call RAPIDDRAIN(id,ICpBtDm(1),ICpTpWaSrDm(1),flDraTub,
     &          FrMpWalWet,FrReduQ,KDCrRlRef,VlMpDmCp,WaSrMpDm(1),
     &          ZBtDm(1),ZWaLevDm(1),FlwOutDrRapCpPot,FlwOutDrRapPot,
     &          VlMpUndrDrL)

! For debugging
!
!         if(itel.eq.7 )then
!         write(102,'(/,i7,a1,a18,200(a1,1pe12.4:))') itel,',',
!     & 'FlwInMtxSatDmCpPot ',(',',FlwInMtxSatDmCpPot(ic) ,ic=1,numnod)
!         write(102,'(i7,a1,a18,200(a1,1pe12.4:))') itel,',',
!     & 'FlwOutMtxSatDmCpPot',(',',FlwOutMtxSatDmCpPot(1,ic),ic=1,numnod)
!         write(102,'(i7,a1,a18,200(a1,1pe12.4:))') itel,',',
!     & 'FlwOutMtxSatDmCpPot',(',',FlwOutMtxSatDmCpPot(2,ic),ic=1,numnod)
!         write(102,'(i7,a1,a18,200(a1,1pe12.4:))') itel,',',
!     & 'FlwOutMtxUnsDmCpPot',(',',FlwOutMtxUnsDmCpPot(1,ic),ic=1,numnod)
!         write(102,'(i7,a1,a18,200(a1,1pe12.4:))') itel,',',
!     & 'FlwOutMtxUnsDmCpPot',(',',FlwOutMtxUnsDmCpPot(2,ic),ic=1,numnod)
!         write(102,'(i7,a1,a18,200(a1,1pe12.4:))') itel,',',
!     & 'FlwOutDrRapCpPot',(',',FlwOutDrRapCpPot(ic),ic=1,numnod)
!         end if

!- B. CALCULATE NEW WATER STORAGE AND CHECK WITH AVAILABLE DOMAIN VOLUME

         FlwInDmTot= FlwInTopPreDmPot(id) + FlwInTopLatDmPot(id) +
     &               FlwInIntSatDmPot + FlwInMtxSatDmPot
         FlwOutDmTot(id)= FlwOutMtxSatDmPot(id) + FlwOutMtxUnsDmPot(id)
         if (id.eq.1) FlwOutDmTot(1)= FlwOutDmTot(1) + FlwOutDrRapPot
!
         WaSrMpDmTmp= WaSrMpDm(id) + FlwInDmTot - FlwOutDmTot(id)
!
         VlMpUndrGwL= 0.d0
         if (ZBtDm(id).lt.GWlFlCpZo .and. GWlFlCpZo.lt.900.d0) ! GWL = 999. indicates NO groundwater level
     &       VlMpUndrGwL= VOLUNDR(id,ICpBtDm(id),GWlFlCpZo,VlMpDmCp,
     &                            ZBtDm(id))
         VlMpUndrGwL= dmin1(VlMpUndrGwL,VlMpDm(id))

         WaSrMpDmMax= VlMpDm(id)
         if((FlwInDmTot-FlwInMtxSatDmPot)/DT.lt.1.d-7 .and.
     &       FlwInMtxSatDmPot/DT.gt.1.d-7) WaSrMpDmMax= VlMpUndrGwL
!
!   1. In case of inflow excess, decrease incoming fluxes
         FrFlwIn= 1.d0
         if (WaSrMpDmTmp.gt.WaSrMpDmMax+1.d-7) then
            FlwInDmExces= WaSrMpDmTmp - WaSrMpDmMax
            FrFlwIn= dmax1(0.d0,1.d0 - FlwInDmExces/FlwInDmTot)
            FlwInTopExcesTot= FlwInTopExcesTot + (1.d0 - FrFlwIn) *
     &                       (FlwInTopPreDmPot(id)+FlwInTopLatDmPot(id))
         endif
         QInTopPreDm(id)= FrFlwIn * FlwInTopPreDmPot(id) / DT
         QInTopLatDm(id)= FrFlwIn * FlwInTopLatDmPot(id) / DT
         do 90 ic= 1, NumNod
            QInIntSatDmCp(id,ic)= FrFlwIn * FlwInIntSatDmCpPot(ic) / DT
            QInMtxSatDmCp(id,ic)= FrFlwIn * FlwInMtxSatDmCpPot(ic) / DT
  90     continue
!
         WaSrMpDmMin= dmin1(VlMpUndrGwL,WaSrMpDm(id))
         if (id.eq.1 .and. VlMpUndrDrL.gt.1.d-6 .and.
     &       VlMpUndrDrL.lt.WaSrMpDmMin) WaSrMpDmMin= VlMpUndrDrL

         FlwOutDmExces(id)= 0.d0
         if (WaSrMpDmTmp.lt.WaSrMpDmMin) FlwOutDmExces(id)=
     &                                   WaSrMpDmMin - WaSrMpDmTmp
!      saturation deficit of domain and of total of domains
         SatDefDm(id)= dmax1(0.d0,(VlMpDm(id)-WaSrMpDmTmp))
         SatDefDmRel(id)= SatDefDm(id) / PpDmCp(id,1)
         SatDefTot= SatDefTot + SatDefDm(id)
 100  continue

!    - Redistribute excess over available volume of domains that are not filled up yet
      if (FlwInTopExcesTot.gt.1.d-6 .and. SatDefTot.gt.1.d-6) then
!      put relative saturation deficit SatDefDmRel(id) in ascending order
         do 200 id= 1, NumDm
            Sq(id)= id
 200     continue
         do 210 id= 1, NumDm-1
            do 209 jd= id+1, NumDm
               if (SatDefDmRel(Sq(id)).gt.SatDefDmRel(Sq(jd))) then
                  SqHlp= Sq(id)
                  Sq(id)= SQ(jd)
                  Sq(jd)= SqHlp
               endif
 209        continue
 210     continue

!    - Redistribute excess
         PpTot= 1.d0
         do 220 id= 1, NumDm
            jd= Sq(id)
            if (SatDefDm(jd).gt.1.d-7 .and.
     &          FlwInTopPreDmPot(jd)+FlwInTopLatDmPot(jd).gt.1.d-7) then
               Pp= PpDmCp(jd,1) / PpTot
               FlwInTopXtr= dmin1(Pp*FlwInTopExcesTot,SatDefDm(jd))
               FrXtr= FlwInTopXtr /
     &                (FlwInTopPreDmPot(jd) + FlwInTopLatDmPot(jd))
               QInTopPreDm(jd)= (1.d0+FrXtr) * FlwInTopPreDmPot(jd) / DT
               QInTopLatDm(jd)= (1.d0+FrXtr) * FlwInTopLatDmPot(jd) / DT
!
               FlwInTopExcesTot = FlwInTopExcesTot - FlwInTopXtr
               FlwOutDmExces(jd)= FlwOutDmExces(jd) - FlwInTopXtr
               FlwOutDmExces(jd)= dmax1(0.d0,FlwOutDmExces(jd))
               SatDefDm(jd)= SatDefDm(jd) - FlwInTopXtr
            else
               SatDefDm(jd)= 0.d0
            endif
            PpTot= PpTot - PpDmCp(jd,1)
 220     continue
      endif

!   2. Calculate actual outgoing fluxes on basis of (possibly) outflow excess
      do 230 id= 1, NumDm
         FrFlwOut= 1.d0
         if (FlwOutDmExces(id)/DT.gt.1.d-7) then
            FrFlwOut= 1.d0 - FlwOutDmExces(id)/FlwOutDmTot(id)
            FrFlwOut= dmax1(0.d0,FrFlwOut)
         endif
         do 228 ic= 1, NumNod
            QOutMtxSatDmCp(id,ic)= FrFlwOut *
     &                             FlwOutMtxSatDmCpPot(id,ic) / DT
            QOutMtxUnsDmCp(id,ic)= FrFlwOut *
     &                             FlwOutMtxUnsDmCpPot(id,ic) / DT
 228     continue
         if (id.eq.1) then
            do 229 ic= 1, NumNod
               QOutDrRapCp(ic)= FrFlwOut * FlwOutDrRapCpPot(ic) / DT
 229        continue
         endif
 230  continue

!   3. Remaining inflow excess is substracted from lateral inflow into macropores
!      for use in main SWAP
      QMpLatSs= QMpLatSs - FlwInTopExcesTot
      if (dabs(QMpLatSs).lt.1.0d-7) QMpLatSs= 0.0d0
!
!- C. CALCULATE FLUXES FOR USE IN SWAP
!
!   1. Aggregate fluxes for output
      do 240 id= 1, NumDm
         do 239 ic= 1, ICpBtDm(id)
            QExcMtxDmCp(id,ic)=
     &                   QOutMtxSatDmCp(id,ic) + QOutMtxUnsDmCp(id,ic) -
     &                   (QInIntSatDmCp(id,ic) + QInMtxSatDmCp(id,ic))
 239     continue
 240  continue

!   2. Aggregate fluxes for main SWAP
      do 250 ic= 1, NumNod
         do 249 id= 1, NumDm
            QExcMpMtx(ic)   = QExcMpMtx(ic) + QExcMtxDmCp(id,ic)
 249     continue
         QMapo= QMapo + QExcMpMtx(ic)
         QRapDra= QRapDra + QOutDrRapCp(ic)
 250  continue
!
      if (ITask.eq.1) return
!              ------------- END of ITask 1 -------------
!
!- D. CALCULATE DERIVATIVES FOR USE IN SWAP (Headcalc)
!
!   1. exchange with unsaturated matrix
      do 260 id= 1, NumDm
         call ABSORPTION(2,ICpBtDm(id),ICpBtPerZon,ICpTpPerZon,
     &          ICpTpSatZon,ICpTpWaSrDm(id),id,AwlCorFac,FrMpWalWet,
     &          FrReduQ,QOutMtxUnsDmCp,SorpDmCp,ThtSrpRefDmCp,
     &          TimAbsCumDmCp,FlwOutMtxUnsDmCpPot,FlwOutMtxUnsDmPot(id),
     &          FlEndSrpEvt,ZWaLevDm)
 260  continue
!
!   2. exchange with saturated matrix
      do 270 id= 1, NumDm
         call SATFLOW(4,ICpBtDm(id),NumNod,ICpSatGWl,ICpTpWaSrDm(id),
     &          ICpTpSatZon,id,FrMpWalWet,FrReduQ,GWL,QInMtxSatDmCp,
     &          QOutMtxSatDmCp,ZWaLevDm(id),FlwInMtxSatDmCpPot,
     &          FlwInMtxSatDmPot,FlwOutMtxSatDmCpPot,
     &          FlwOutMtxSatDmPot(id))

 270  continue

      return
      END

!=====================================================================
      SUBROUTINE SATFLOW(ITask,ICpBtDm,ICpBtZon,ICpSatLev,ICpTpWaSrDm,
     &              ICpTpZon,id,FrMpWalWet,FrReduQ,Lev,QInMtxSatDmCp,
     &              QOutMtxSatDmCp,ZWaLevDm,
     &              FlwInDmCpPot,FlwInDmPot,FlwOutDmCpPot,FlwOutDmPot)
! ----------------------------------------------------------------------
!     Date               : April 2008
!     Purpose            : To calculate water exchange between macropores
!                        : and the SATURATED soil matrix
!     Functions called   : -
! ----------------------------------------------------------------------
      use Variables
      implicit NONE

! --- global                                                          In
      integer id, ICpBtDm, ICpBtZon,ICpSatLev,ICpTpWaSrDm,ICpTpZon,ITask
      real*8  FrMpWalWet(MaDm,MaCp), FrReduQ, Lev, ZWaLevDm
      real*8  QInMtxSatDmCp(MaDm,MaCp), QOutMtxSatDmCp(MaDm,MaCp)
!     -                                                              Out
      real*8  FlwInDmCpPot(MaCp), FlwInDmPot
      real*8  FlwOutDmCpPot(MaDm,MaCp), FlwOutDmPot
! ----------------------------------------------------------------------
! --- local
      integer ic, Lay
      real*8  DelH, DelHDmCp(MaDm,MaCp), FlwMtxSatDmCp(MaCp), HMp, HMa
      real*8  ksatHor, Pi, RatSeepFace, RecRes, ResHor, ResRad,ResVrt

      data Pi, RatSeepFace / 3.14159d0, 10.d0 /
! ----------------------------------------------------------------------
      goto (1000,2000,3000,4000) itask

! --- Calculate potential water EXCHANGE between macropores and saturated matrix
!
 1000 continue
! - A. Calculate potential flows INCOMING into macropores
!
! --- Initialisation global Flows and contribution Derivatives
      do 10 ic= 1, NumNod
         FlwMtxSatDmCp(ic)= 0.d0
         DelHDmCp(id,ic)  = 0.d0
  10  continue

 2000 continue
! --- Initialisation local Flows
      FlwInDmPot= 0.d0
      do 15 ic= 1, NumNod
         FlwInDmCpPot(ic) = 0.d0
  15  continue

! --- Check whether this subroutine is relevant for this macropore domain
      if (ICpBtDm.ge.ICpTpZon .and. ICpBtZon.gt.0) then
         continue
      else
         return
      endif
!
!   - ICpTpZon = top compartment of relevant saturated zone
      do 100 ic= ICpTpZon, min0(ICpBtDm,ICpBtZon)
         Lay= Layer(ic)
         ksatHor=  ksatfit(Lay)
!
!   - First calculation of Heads in Macropores and Matrix, for determing flow direction
         HMp = ZWaLevDm - Z(ic)  ! head in macropore domain, assuming static equilibrium
         HMp = dmax1(0.d0,HMp)   ! no negative head in macropores
         if (HMp.lt.1.d-8) HMp = 0.d0
         HMa = H(ic)
         DelH= HMp - HMa         ! differences between heads determines flow direction
         if (HMp.lt.1.d-8 .and. DelH.gt.0.d0) DelH= 0.d0
         if (abs(DelH).lt.1.d-8) DelH = 0.d0

! - 3 flow situations:
!    1. DelH = +  -> infiltration of macropore water into matrix; Darcy flow
!    2. and 3. DelH = -  -> exfiltration of matrix water into macropores
!    2.        HMp > 0   -> Darcy flow
!    3         HMp = 0   -> seepage face; Hooghoudt above drain level
         if (Delh.gt.0.d0) then
! - 1. infiltration out of macropores into matrix: = term 2.a. (Darcy)
!           RecRes = reciprocal resistance for calculation of derivative
            RecRes = ShapeFacMp * 8.0d0 * PpDmCp(id,ic) * DZ(ic) *
     &               ksatHor / DiPoCp(ic)**2.d0
            if (ic.eq.ICpTpWaSrDm) RecRes = FrMpWalWet(id,ic) * RecRes
            FlwMtxSatDmCp(ic) = - FrReduQ * RecRes * DelH * DT
            DelHDmCp(id,ic)= DelH  ! save head difference for calculation of derivative
!
         elseif (DelH.lt.0.d0) then
!   - exfiltration out of matrix into macropores
            if (HMp.gt.0.d0) then
! - 2. flow equation according to Darcy
               RecRes = ShapeFacMp * 8.0d0 * PpDmCp(id,ic) * DZ(ic) *
     &                                       ksatHor / DiPoCp(ic)**2.d0
               if (ic.eq.ICpSatLev) RecRes = RecRes *
     &                       (GWlFlCpZo - (Z(ic)-0.5d0*DZ(ic))) / DZ(ic)
               DelHDmCp(id,ic)= DelH ! save head difference for calculation of derivative
            else
! - 3. flow equation according to seepage face: Ernst with radial resistance
!               RecRes = -4.0d0 * PpDmCp(id,ic) * ksatHor * DelH / DiPoCp(ic)**2.d0
!               DelHDmCp(id,ic)= DelH / 2.d0 ! div. by 2 because of DelH**2
               ResHor = DiPoCp(ic)**2.d0 / (8.0d0 * DZ(ic) * ksatHor)
               ResVrt = DZ(ic) / ksatHor
               ResRad = DiPoCp(ic) * DLog(RatSeepFace) / (Pi*ksatHor)
               RecRes = PpDmCp(id,ic) / (ResHor + ResVrt + ResRad)
               if (ic.eq.ICpSatLev) RecRes = RecRes *
     &                             (Lev - (Z(ic)-0.5d0*DZ(ic))) / DZ(ic)
               DelHDmCp(id,ic)= DelH ! save head difference for calculation of derivative
            endif
            FlwMtxSatDmCp(ic) = - FrReduQ * RecRes * DelH * DT
         endif
!   - save positive FlwMtxSatDmCp as INCOMING (into macropores) flux
         if (FlwMtxSatDmCp(ic).gt.0.d0) then
            FlwInDmCpPot(ic) = FlwMtxSatDmCp(ic)
            FlwInDmPot       = FlwInDmPot + FlwInDmCpPot(ic)
         endif
 100  continue

      return

 3000 continue
! - B. Calculate potential flows OUTGOING out off macropores
! --- Initialisation Flows
      FlwOutDmPot = 0.d0
      do 200 ic= 1, NumNod
         FlwOutDmCpPot(id,ic)= 0.d0
 200  continue

! --- Check whether this subroutine is relevant for this macropore domain
      if (ICpBtDm.ge.ICpTpZon) then
         continue
      else
         return
      endif

!   - save negative FlwMtxSatDmCp as OUTGOING (out off macropores) flux
!   - (calculation of FlwMtxSatDmCp(ic) done above, in ITask 1)
      do 250 ic= ICpTpZon, ICpBtDm
         if (FlwMtxSatDmCp(ic).lt.0.d0) then
            FlwOutDmCpPot(id,ic)= - FlwMtxSatDmCp(ic)
            FlwOutDmPot         = FlwOutDmPot + FlwOutDmCpPot(id,ic)
         endif
 250  continue

      return

 4000 continue
! - C. Determine derivatives

      do 300 ic= ICpTpZon, ICpBtDm
         if (abs(DelHDmCp(id,ic)).gt.1.d-14) then
            dFdhMp(ic)= dFdhMp(ic) -
     &                  (QOutMtxSatDmCp(id,ic) - QInMtxSatDmCp(id,ic)) /
     &                  DelHDmCp(id,ic)
         end if
 300  continue

      Return
      End

!=====================================================================
      SUBROUTINE ABSORPTION(ITask,ICpBtDm,ICpBtPerZon,ICpTpPerZon,
     &              ICpTpSatZon,ICpTpWaSrDm,id,AwlCorFac,FrMpWalWet,
     &              FrReduQ,QOutMtxUnsDmCp,SorpDmCp,ThtSrpRefDmCp,
     &              TimAbsCumDmCp,FlwOutMtxUnsDmCpPot,FlwOutMtxUnsDmPot,
     &              FlEndSrpEvt,ZWaLevDm)

! ----------------------------------------------------------------------
!     Date               : May 2005
!     Purpose            : To calculate absorption of macropore water into the
!                          unsaturated matrix by sorptivity under dry conditions
!                          and by pressure head gradient with Darcy
!     Functions called   : moiscap
! ----------------------------------------------------------------------
      use Variables
      implicit NONE

! --- global                                                          In
      integer id, ICpBtDm, ICpBtPerZon, ICpTpPerZon, ICpTpSatZon
      integer ICpTpWaSrDm, ITask
      real*8  AwlCorFac(Macp), FrMpWalWet(MaDm,MaCp)
      real*8  QOutMtxUnsDmCp(MaDm,MaCp), SorpDmCp(MaDm,MaCp)
      real*8  ThtSrpRefDmCp(MaDm,MaCp), TimAbsCumDmCp(MaDm,MaCp)
      real*8  ZWaLevDm(MaDm), FrReduQ
!     -                                                              Out
      real*8  FlwOutMtxUnsDmCpPot(MaDm,MaCp), FlwOutMtxUnsDmPot
      logical FlEndSrpEvt(MaDm,MaCp)
! ----------------------------------------------------------------------
! --- local
      integer ic, ICpBtUnsMtxDm, Lay
      real*8  AbsDarc, AbsSorp, CritSatDefMtx, DelH, DelHDmCp(MaDm,MaCp)
      real*8  Deriv, HMp, moiscap, QSorpMax, RecRes, SatDefMtx
      real*8  SorpAct, SorpFac, Time
      logical FlSorp(MaDm,MaCp)

!      integer tel
! ----------------------------------------------------------------------
! --- Initialization Flows and Flags
      FlwOutMtxUnsDmPot= 0.d0
      do 10 ic= 1, NumNod
         FlwOutMtxUnsDmCpPot(id,ic)= 0.d0
         FlEndSrpEvt(id,ic)        = .true.
  10  continue

! --- Initialization help variables
      ICpBtUnsMtxDm= min0(ICpBtDm,ICpTpSatZon-1)
      CritSatDefMtx= 1.0d-8
      QSorpMax     = 1.0d3
!
      goto (1000,2000) itask

 1000 continue
! --- Calculate potential Absorption
! - I By Sorptivity
      do 100 ic= ICpTpWaSrDm, ICpBtUnsMtxDm
         if (ic.gt.ICpBtPerZon .or. ic.lt.ICpTpPerZon) then
!
            Lay      = Layer(ic)
            Time     = TimAbsCumDmCp(id,ic)
            SatDefMtx= ThetaS(ic) - Theta(ic)
            if (SatDefMtx.lt.CritSatDefMtx) then
!    -  Matrix too wet: no sorptivity
               AbsSorp = 0.d0
            else
               if (Time.lt.1.d-8) then
!    -  Start a new sorptivity event; set ThtSrpRefDmCp(id,ic) and SorpDmCp
!    - ThtSrpRefDmCp = Thet_sat + Delta_thet_theor. = Thet_sat + Thet_theor. - Thet_0
!      Delta_thet_theor. = theoretical increase of theta (calculated in MacroState)
                  ThtSrpRefDmCp(id,ic)= ThetaS(ic)  ! initial, delta_theta = 0
                  SorpDmCp(id,ic)= SorpMax(Lay) *
     &                          ((ThetaS(ic)-Theta(ic)) /
     &                          (ThetaS(ic)-ThetaR(ic)))** SorpAlfa(Lay)
                  SorpAct = SorpDmCp(id,ic)
               elseif((ThtSrpRefDmCp(id,ic)-Theta(ic)).gt.CritSatDefMtx)
     &         then
!    -  Continue existing sorptivity event; set SorpAct
                  SorpAct= SorpMax(Lay) *
     &                     ((ThtSrpRefDmCp(id,ic)-Theta(ic)) /
     &                     (ThetaS(ic)-ThetaR(ic)))** SorpAlfa(Lay)
               else
!    -  End existing sorptivity event
                  SorpAct= 0.d0
               endif
               AbsSorp = SorpAct * PpDmCp(id,ic) * (4.d0 *AwlCorFac(ic)*
     &               DZ(ic) / DiPoCp(ic)) * (dsqrt(Time+DT)-dsqrt(Time))
               if (ic.eq.ICpTpWaSrDm) AbsSorp= FrMpWalWet(id,ic)*AbsSorp
! --- Limit peak amount of absorbed water at beginning of sorptiviy event
               AbsSorp = dmin1(QSorpMax*DT*DZ(ic),AbsSorp)
            endif

!   II By Darcy
            HMp = ZWaLevDm(id) - Z(ic)  ! head in macropore domain, assuming static equilibrium
            HMp = dmax1(0.d0,HMp)       ! no negative head in macropores
            if (H(ic).lt.-1.d-7 .and. HMp.gt.1.d-7) then
               DelH= dmax1(0.d0,HMp - H(ic))
            else
               DelH= 0.d0
            endif
            DelHDmCp(id,ic)= DelH     ! save head difference for calculation of derivative
                                      ! reciporocal resistance for calculation of derivative
            RecRes  = ShapeFacMp * 8.0d0 * PpDmCp(id,ic) * DZ(ic) *
     &                K(ic) / DiPoCp(ic)**2.d0
            if (ic.eq.ICpTpWaSrDm) RecRes = FrMpWalWet(id,ic) * RecRes
            AbsDarc = RecRes * DelH * DT
!
            SorpFac = SorpFacParl(lay) + (1.d0-SorpFacParl(lay)) *
     &                ( 1.d0 -(dmax1((ThetaS(ic)-Theta(ic)),0.d0) /
     &                         (ThetaS(ic)-ThetaR(ic)))**SorpAlfa(Lay) )

            if (SwDarcy.eq.0) AbsDarc = 0.d0

            if (AbsSorp.gt.SorpFac*AbsDarc) then
               FlwOutMtxUnsDmCpPot(id,ic) = FrReduQ * AbsSorp
               FlSorp(id,ic) = .true.
               AbsDarc = 0.d0
!            tel = tel + 1
!            write(99,99) t, absdarc/dt, abssorp/dt, theta(ic), time,ic,tel
            else
!            tel = tel + 1
!            write(99,99) t, absdarc/dt, abssorp/dt, theta(ic), time,ic,tel
               FlwOutMtxUnsDmCpPot(id,ic) = FrReduQ * SorpFac*AbsDarc
!            FlwOutMtxUnsDmCpPot(id,ic) = AbsDarc
               AbsSorp = 0.d0
               FlSorp(id,ic) = .false.
            endif
!
            FlwOutMtxUnsDmPot= FlwOutMtxUnsDmPot +
     &                         FlwOutMtxUnsDmCpPot(id,ic)
!
! --- Set Flag to END SoRPtivity EVenT in MACROSTATE
            if (AbsSorp/DT.gt.1.d-7) then
               FlEndSrpEvt(id,ic)= .false.
            endif
         endif
!
 100  continue
!
      return
!
 2000 continue
! --- Determine derivatives
      do ic= ICpTpWaSrDm, ICpBtUnsMtxDm
         if (ic.gt.ICpBtPerZon .or. ic.lt.ICpTpPerZon) then
!
            if (QOutMtxUnsDmCp(id,ic).gt.1.d-7) then
               if (FlSorp(id,ic)) then
                  Lay   = Layer(ic)
                  Deriv = - QOutMtxUnsDmCp(id,ic) * SorpAlfa(Lay) /
     &                    (ThtSrpRefDmCp(id,ic) - Theta(ic))
                  Deriv = Deriv * moiscap(ic,H(ic),CofGen(1,ic),dt,
     &                                             swsophy,numtab,sptab)
               else
                  Deriv = 0.d0
                  if (abs(DelHDmCp(id,ic)).gt.1.d-14) then
                     Deriv = - QOutMtxUnsDmCp(id,ic) / DelHDmCp(id,ic)
                  endif
               endif
               dFdhMp(ic)= dFdhMp(ic) + Deriv
            endif
         endif
      end do
!
      return
      end

!=====================================================================
      SUBROUTINE RAPIDDRAIN(id,ICpBtDm,ICpTpWaSrDm,flDraTub,FrMpWalWet,
     &              FrReduQ,KDCrRlRef,VlMpDmCp,WaSrMpDm,ZBtDm,ZWaLevDm,
     &              FlwOutDrRapCpPot,FlwOutDrRapPot,VlMpUndrDrL)
! ----------------------------------------------------------------------
!     Date               : April 2008
!     Purpose            : To calculate rapid drainage from Main Bypass flow domain
! ----------------------------------------------------------------------
      use Variables
      implicit NONE

! --- global                                                          In
      integer ICpBtDm, ICpTpWaSrDm, id
      real*8  FrMpWalWet(MaDm,MaCp), FrReduQ, KDCrRlRef(MaDr)
      real*8  VlMpDmCp(MaDm,MaCp), WaSrMpDm, ZBtDm, ZWaLevDm
      logical flDraTub(Madr)
!     -                                                              Out
      real*8  FlwOutDrRapCpPot(MaCp), FlwOutDrRapPot, VlMpUndrDrL
      real*8  FacRes
! ----------------------------------------------------------------------
! --- local
      integer ic
      real*8  DelH, KDCrRl, KDCrRlCp(MaCp), RapDraRes
      real*8  VOLUNDR, WaSrDrainabl, WthCr
! ----------------------------------------------------------------------
!
! --- Check whether this subroutine is relevant for this macropore domain
!     Rapid drainage only in domain 1: Main Bypass flow domain
      if (id.eq.1) then
!   - Initialization
         FlwOutDrRapPot= 0.d0
         do 100 ic= 1, NumNod
            FlwOutDrRapCpPot(ic)= 0.d0
 100     continue
!
         if (SwDrRap.eq.1) then
            continue
         else
            return
         endif
      else
         return
      endif
!
! --- for drain tube rapid drainage only possible if ZBtDm .le. Zdrabas;
!   - for open drain, rapid drainage also possible when ZBtDm > Zdrabas
      if (ZBtDm.Lt.ZDraBas+1.0d-1 .or. .not.flDraTub(1)) then
! transmissivity kD
         KDCrRl= 0.d0
         do 200 ic= ICpTpWaSrDm, ICpBtDm
            WthCr= DiPoCp(ic) * (1.d0 - dsqrt(1.d0-VlMpDmCp(1,ic) /
     &             DZ(ic)))
            KDCrRlCp(ic)= ((WthCr**RapDraReaExp)/DiPoCp(ic)) * DZ(ic)
            if (ic.eq.ICpTpWaSrDm) KDCrRlCp(ic)=
     &                                FrMpWalWet(id,ic) * KDCrRlCp(ic)
            KDCrRl= KDCrRl + KDCrRlCp(ic)
 200     continue
!
         VlMpUndrDrL= 0.d0
         if (ZBtDm.lt.ZDraBas) VlMpUndrDrL= VOLUNDR(
     &                               1,ICpBtDm,ZDraBas,VlMpDmCp,ZBtDm)
         WaSrDrainabl= dmax1(0.d0,(WaSrMpDm - VlMpUndrDrL))
!
! for open drain: DelH =  zwalevdm - ditch water level or bottom of MB domain !
         DelH=  ZWaLevDm - dmax1(ZDraBas,ZBtDm)
         if (ZWaLevDm.gt.-1.d-7) DelH= DelH + Pond
! ! in this version: no infiltration
         DelH= dmax1(DelH,0.d0)
         if (KdCrRl.gt.1.d-10) then
            FacRes        = dmin1(KDCrRlRef(1)/KDCrRl,1.1d0)
            RapDraRes     = RapDraResRef(1) * FacRes
            FlwOutDrRapPot= FrReduQ * (DelH / RapDraRes) * DT
         else
            FlwOutDrRapPot= 0.d0
         endif
         FlwOutDrRapPot= dmin1(FlwOutDrRapPot,WaSrDrainabl)
! distribute over compartments according to kD
         do 300 ic= ICpTpWaSrDm, ICpBtDm
            if (KDCrRl.gt.1.d-15) then
               FlwOutDrRapCpPot(ic)= FlwOutDrRapPot * KDCrRlCp(ic) /
     &                               KDCrRl
            else
               FlwOutDrRapCpPot(ic)= 0.d0
               FlwOutDrRapPot      = 0.d0
            endif
 300     continue
      endif

      return
      END

!=======================================================================
      real*8 FUNCTION VOLUNDR(id,ICpBtDm,Level,VlMpDmCp,ZBtDm)
! ----------------------------------------------------------------------
!     Date               : september 2002
!     Purpose            : calculates volume of macropores under drain or
!                        : groundwater level
! ----------------------------------------------------------------------
      use Variables
      implicit NONE

! --- global
      integer ICpBtDm, id
      real*8  Level, VlMpDmCp(MaDm,MaCp), ZBtDm
! ----------------------------------------------------------------------
! --- local
      integer ic
      real*8  ZHlp
! ----------------------------------------------------------------------
! --- macropore volume under drain or groundwater level
      VOLUNDR= 0.d0
      ic= ICpBtDm + 1
      ZHlp= ZBtDm
      do 10 while(ZHlp.lt.Level .and. ic.gt.1)
         ic= ic - 1
         ZHlp= ZHlp + DZ(ic)
         VOLUNDR= VOLUNDR + VlMpDmCp(id,ic)
  10  continue
      if (VOLUNDR.gt.0.d0) VOLUNDR =
     &                     VOLUNDR - (ZHlp-Level)*VlMpDmCp(id,ic)/DZ(ic)
      return
      END

ext: for


text:
      FUNCTION rtsec(func,x1,x2,xacc)
      INTEGER MAXIT
      REAL rtsec,x1,x2,xacc,func
      EXTERNAL func
      PARAMETER (MAXIT=30)
      INTEGER j
      REAL dx,f,fl,swap,xl
      fl=func(x1)
      f=func(x2)
      if(abs(fl).lt.abs(f))then
        rtsec=x1
        xl=x2
        swap=fl
        fl=f
        f=swap
      else
        xl=x1
        rtsec=x2
      endif
      do 11 j=1,MAXIT
        dx=(xl-rtsec)*f/(f-fl)
        xl=rtsec
        fl=f
        rtsec=rtsec+dx
        f=func(rtsec)
        if(abs(dx).lt.xacc.or.f.eq.0.)return
11    continue
      pause 'rtsec exceed maximum iterations'
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/rtsec.for


text:
      subroutine fill_diff_coeff(lo,hi, &
                                 coefx,coefx_l1,coefx_l2,coefx_h1,coefx_h2,&
                                 coefy,coefy_l1,coefy_l2,coefy_h1,coefy_h2,dx)

      use probdata_module, only: diff_coeff

      implicit none
      integer         , intent(in   ) :: lo(2),hi(2)
      integer         , intent(in   ) :: coefx_l1,coefx_l2,coefx_h1,coefx_h2
      integer         , intent(in   ) :: coefy_l1,coefy_l2,coefy_h1,coefy_h2
      double precision, intent(inout) :: coefx(coefx_l1:coefx_h1,coefx_l2:coefx_h2)
      double precision, intent(inout) :: coefy(coefy_l1:coefy_h1,coefy_l2:coefy_h2)
      double precision, intent(in   ) :: dx(2)

      ! Local variables
      integer          :: i,j

      do j = lo(2),hi(2)
        do i = lo(1),hi(1)+1
          coefx(i,j) = diff_coeff
        end do
      end do

      do j = lo(2),hi(2)+1
        do i = lo(1),hi(1)
          coefy(i,j) = diff_coeff
        end do
      end do

      end subroutine fill_diff_coeff

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


      subroutine qpot(Nb,ddv,c,am,p,x,w,u,du)
      use cdat
      implicit real*8(a-h,o-z)
      integer*4,intent(IN) :: Nb
      real*8, intent(OUT) :: du(Nb),u(Nb)
      integer INFO
      real*8 :: w(Nb),x(Nb),c(nb),p(nb),u2(nb),du2(nb),r0(nb,nb)
      real*8 :: f(2),rx(Nb),ddv(nb),s(2,2),cr(2,1)

! defn = 0d0
      x1 = 0d0
      x2 = 0d0
      do i=1,nb
        an=an+w(i)
        x1=x1+w(i)*x(i)
        x2=x2+w(i)*x(i)**2
      enddo

      b1=-0.5d0*an*an/(an*x2-x1*x1)
      b2 = 0.5d0*x1*an/(an*x2-x1*x1)
      do i=1,Nb
        r=b1*x(i)+b2
        u(i) = -(r*r+b1)/2d0/am
        du(i) = -b1*r/am
      enddo
!ine c matrix
!        cr(1,1)=-0.5d0
!        cr(2,1)=0d0
!
!! quantum force
!        s=0d0
!        do i=1,Ntraj
!          f=(/x(i),1d0/)
!        do m=1,2
!          do n=1,2
!            s(m,n)=w(i)*f(m)*f(n)+s(m,n)
!          end do
!        end do
!        end do
!
!! calculate matrix c(t)
!        call DPOSV('U',2,1,s,2,cr,2,INFO)
!        if(INFO/=0) then
!        print *, "info=",info
!        print *, "QPOT: matrix fails"
!        stop
!        end if
!
!! the momentum operator r=cf
!        do i=1,Ntraj
!          rx(i)=cr(1,1)*x(i)+cr(2,1)
!
!! calculate quantum potential
!          u(i)=-rx(i)**2/(2d0*am)-cr(1,1)/(2d0*am)
!          du(i)=-rx(i)*cr(1,1)/am
!        end do
!      u = 0d0
!      du = 0d0
      call rep(nb,w,x,p,u2,du2)
!      u2 = 0d0
      du2 = 0d0
      do i=1,nb
!        u(i) = u(i)+u2(i)
        du(i) = du(i)+du2(i)
      enddo

      return
      end subroutine

!------------------------------------------------------
!------extra repulsion added like quantum potential
!----------------------------------------------------
      subroutine rep(nb,w,x,p,u2,du2)
      use cdat
      implicit real*8(a-h,o-z)
      real*8 x(nb),p(nb),dxp(nb),r0(nb,nb),w(nb),u2(nb),du2(nb)
      ww = 0d0
      do i=1,nb
        ww = w(i)+ww
      enddo

      c = 8d0
      u2 = 0d0
      d = 0.32d0

! --- parameters for repulsive force

      r = x(2)-x(1)
      du2(1) = d*exp(-c*r)/r

      r = x(nb)-x(nb-1)
      du2(nb) = -d*exp(-c*r)/r

      do i=2,nb-1
        r1 = x(i+1)-x(i)
        r2 = x(i)-x(i-1)
        du2(i) = -d*exp(-c*r2)/r2+d*exp(-c*r1)/r1
      enddo

!      c = 1d6
!      d = 1d3
!
!      r = x(2)-x(1)
!      du2(1) = d*exp(-c*r**2)
!
!      r = x(nb)-x(nb-1)
!      du2(nb) = -d*exp(-c*r**2)
!
!      do i=2,nb-1
!          r1 = x(i+1)-x(i)
!          r2 = x(i)-x(i-1)
!          du2(i) = -d*exp(-c*r2**2)+d*exp(-c*r1**2)
!      enddo

!      a = 0.0004d0
!      d = 16d1
!
!      r = x(2)-x(1)
!      du2(1) = d/cosh(r/a)**2
!
!      r = x(nb)-x(nb-1)
!      du2(nb) = -d/cosh(r/a)**2
!
!      do i=2,nb-1
!          r1 = x(i+1)-x(i)
!          r2 = x(i)-x(i-1)
!          du2(i) = -d/cosh(r2/a)**2+d/cosh(r1/a)**2
!      enddo

      return
      end subroutine


!***********************************/
!	Name:非構造格子用EulerSolverのための，要素中心界から面中心までの距離ベクトルを計算するプログラム
!	Alias:UCalcWidthCell2Edge
!	Description:セルごとに界面中心から面中心までの距離を求めると同時に，セルごとの平均距離を求める(Venkatakrishnanの流束制限関数用)
!	Type:UnstructuredGrid
!	Input:UG%GM%Dimension,UG%CD%Edge,UG%CD%Cell (xxxは使用する種類のセルデータ)
!	Output:UG%GM%Normal
!	Note:UG%GM%Widthは引数が3つもあってややこしいので気をつけること
!	Author:Akitaka Toyota
!	Date:2017.11.10
!	Update:2018.02.03 !delete AverageWidth
!	Other:現時点でハイブリッド格子・3次元格子への対応は未定(準備はある)
!***********************************/
subroutine UCalcWidthCell2Edge(UG)
    use StructVar_Mod
    use LoopVar_Mod
    use FrequentOperation
    implicit none
    type(UnstructuredGrid), intent(inout) :: UG

    if(UG%GM%Dimension == 2) then
        call WidthTriangle2Line
        if(UG%Quad%Total /= 0) then
            write(6,*) "coming soon...?"
        end if
    else
        write(6,*) "coming soon..."
        stop
    end if
return

contains
    subroutine WidthTriangle2Line
    implicit none
        !UG%GM%AverageWidth = 0.0d0
        do iCell=1, UG%Tri%Total
            do iLocalEdge=1,3
                iEdge = UG%Tri%Edge(iCell,iLocalEdge)
                UG%GM%Width(iCell,iLocalEdge,:) = UG%CD%Edge(iEdge,:) - UG%CD%Cell(iCell,:)
                !UG%GM%AverageWidth(iCell) = UG%GM%AverageWidth(iCell) + sqrt(sum(UG%GM%Width(iCell,iLocalEdge,:)**2)) !第2項はベクトルの絶対値
            end do
            !UG%GM%AverageWidth(iCell) = UG%GM%AverageWidth(iCell)/3.0d0
        end do
        return
    end subroutine WidthTriangle2Line

end subroutine UCalcWidthCell2Edge


text:
!=======================================================================
! This is part of the 2DECOMP&FFT library
!
! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil)
! decomposition. It also implements a highly scalable distributed
! three-dimensional Fast Fourier Transform (FFT).
!
! Copyright (C) 2009-2011 Ning Li, the Numerical Algorithms Group (NAG)
!
!=======================================================================

! This file contain common code to be included by subroutines
! 'write_every_...' in io.f90

  ! To write every few points of a 3D array to a file

    ! work out the distribution parameters, which may be different from
    ! the default distribution used by the decomposition library
    !  For exmample if nx=17 and p_row=4
    !    distribution is: 4 4 4 5

    ! If writing from the 1st element
    !  If saving every 3 points, then 5 points to be saved (17/3)
    !    default distribution would be 1 1 1 2
    !    However, 1st block (1-4) contains the 3rd point
    !             2nd block (5-8) contains the 6th point
    !             3rd block (9-12) contains the 9th and 12th point
    !             4th block (13-17) contains then 15th point
    !    giving a 1 1 2 1 distribution
    !    So cannot use the base decomposition library for such IO

    ! If writing from the n-th element (n=?skip)
    !  If saving every 3 points, then 6 points to be saved
    !    However, 1st block (1-4) contains the 1st & 4th point
    !             2nd block (5-8) contains the 7th point
    !             3rd block (9-12) contains the 10th point
    !             4th block (13-17) contains then 12th & 15th point
    !    giving a 1 2 2 1 distribution

    skip(1)=iskip
    skip(2)=jskip
    skip(3)=kskip

    do i=1,3
       if (from1) then
          xst(i) = (xstart(i)+skip(i)-1)/skip(i)
          if (mod(xstart(i)+skip(i)-1,skip(i))/=0) xst(i)=xst(i)+1
          xen(i) = (xend(i)+skip(i)-1)/skip(i)
       else
          xst(i) = xstart(i)/skip(i)
          if (mod(xstart(i),skip(i))/=0) xst(i)=xst(i)+1
          xen(i) = xend(i)/skip(i)
       end if
       xsz(i) = xen(i)-xst(i)+1
    end do

    do i=1,3
       if (from1) then
          yst(i) = (ystart(i)+skip(i)-1)/skip(i)
          if (mod(ystart(i)+skip(i)-1,skip(i))/=0) yst(i)=yst(i)+1
          yen(i) = (yend(i)+skip(i)-1)/skip(i)
       else
          yst(i) = ystart(i)/skip(i)
          if (mod(ystart(i),skip(i))/=0) yst(i)=yst(i)+1
          yen(i) = yend(i)/skip(i)
       end if
       ysz(i) = yen(i)-yst(i)+1
    end do

    do i=1,3
       if (from1) then
          zst(i) = (zstart(i)+skip(i)-1)/skip(i)
          if (mod(zstart(i)+skip(i)-1,skip(i))/=0) zst(i)=zst(i)+1
          zen(i) = (zend(i)+skip(i)-1)/skip(i)
       else
          zst(i) = zstart(i)/skip(i)
          if (mod(zstart(i),skip(i))/=0) zst(i)=zst(i)+1
          zen(i) = zend(i)/skip(i)
       end if
       zsz(i) = zen(i)-zst(i)+1
    end do

    ! if 'skip' value is large it is possible that some ranks do not
    ! contain any points to be written. Subarray constructor requires
    ! nonzero size so it is not possible to use MPI_COMM_WORLD for IO.
    ! Create a sub communicator for this...
    color = 1
    key = 0  ! rank order doesn't matter
    if (ipencil==1) then
       if (xsz(1)==0 .or. xsz(2)==0 .or. xsz(3)==0) then
          color = 2
       end if
    else if (ipencil==2) then
       if (ysz(1)==0 .or. ysz(2)==0 .or. ysz(3)==0) then
          color = 2
       end if
    else if (ipencil==3) then
       if (zsz(1)==0 .or. zsz(2)==0 .or. zsz(3)==0) then
          color = 2
       end if
    end if
    call MPI_COMM_SPLIT(MPI_COMM_WORLD,color,key,newcomm,ierror)

    if (color==1) then ! only ranks in this group do IO collectively

       ! generate subarray information
       sizes(1) = xsz(1)
       sizes(2) = ysz(2)
       sizes(3) = zsz(3)
       if (ipencil==1) then
          subsizes(1) = xsz(1)
          subsizes(2) = xsz(2)
          subsizes(3) = xsz(3)
          starts(1) = xst(1)-1
          starts(2) = xst(2)-1
          starts(3) = xst(3)-1
       else if (ipencil==2) then
          subsizes(1) = ysz(1)
          subsizes(2) = ysz(2)
          subsizes(3) = ysz(3)
          starts(1) = yst(1)-1
          starts(2) = yst(2)-1
          starts(3) = yst(3)-1
       else if (ipencil==3) then
          subsizes(1) = zsz(1)
          subsizes(2) = zsz(2)
          subsizes(3) = zsz(3)
          starts(1) = zst(1)-1
          starts(2) = zst(2)-1
          starts(3) = zst(3)-1
       end if

       ! copy data from original array
       ! needs a copy of original array in global coordinate
       if (ipencil==1) then
          allocate(wk(xst(1):xen(1),xst(2):xen(2),xst(3):xen(3)))
          allocate(wk2(xstart(1):xend(1),xstart(2):xend(2),xstart(3):xend(3)))
          wk2=var
          if (from1) then
             do k=xst(3),xen(3)
                do j=xst(2),xen(2)
                   do i=xst(1),xen(1)
                      wk(i,j,k) = wk2((i-1)*iskip+1,(j-1)*jskip+1,(k-1)*kskip+1)
                   end do
                end do
             end do
          else
             do k=xst(3),xen(3)
                do j=xst(2),xen(2)
                   do i=xst(1),xen(1)
                      wk(i,j,k) = wk2(i*iskip,j*jskip,k*kskip)
                   end do
                end do
             end do
          end if
       else if (ipencil==2) then
          allocate(wk(yst(1):yen(1),yst(2):yen(2),yst(3):yen(3)))
          allocate(wk2(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
          wk2=var
          if (from1) then
             do k=yst(3),yen(3)
                do j=yst(2),yen(2)
                   do i=yst(1),yen(1)
                      wk(i,j,k) = wk2((i-1)*iskip+1,(j-1)*jskip+1,(k-1)*kskip+1)
                   end do
                end do
             end do
          else
             do k=yst(3),yen(3)
                do j=yst(2),yen(2)
                   do i=yst(1),yen(1)
                      wk(i,j,k) = wk2(i*iskip,j*jskip,k*kskip)
                   end do
                end do
             end do
          end if
       else if (ipencil==3) then
          allocate(wk(zst(1):zen(1),zst(2):zen(2),zst(3):zen(3)))
          allocate(wk2(zstart(1):zend(1),zstart(2):zend(2),zstart(3):zend(3)))
          wk2=var
          if (from1) then
             do k=zst(3),zen(3)
                do j=zst(2),zen(2)
                   do i=zst(1),zen(1)
                      wk(i,j,k) = wk2((i-1)*iskip+1,(j-1)*jskip+1,(k-1)*kskip+1)
                   end do
                end do
             end do
          else
             do k=zst(3),zen(3)
                do j=zst(2),zen(2)
                   do i=zst(1),zen(1)
                      wk(i,j,k) = wk2(i*iskip,j*jskip,k*kskip)
                   end do
                end do
             end do
          end if
       end if
       deallocate(wk2)

       ! MPI-IO
       call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts,  &
            MPI_ORDER_FORTRAN, data_type, newtype, ierror)
       call MPI_TYPE_COMMIT(newtype,ierror)
       call MPI_FILE_OPEN(newcomm, filename, &
            MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, &
            fh, ierror)
       filesize = 0_MPI_OFFSET_KIND
       call MPI_FILE_SET_SIZE(fh,filesize,ierror)  ! guarantee overwriting
       disp = 0_MPI_OFFSET_KIND
       call MPI_FILE_SET_VIEW(fh,disp,data_type, &
            newtype,'native',MPI_INFO_NULL,ierror)
       call MPI_FILE_WRITE_ALL(fh, wk, &
            subsizes(1)*subsizes(2)*subsizes(3), &
            data_type, MPI_STATUS_IGNORE, ierror)
       call MPI_FILE_CLOSE(fh,ierror)
       call MPI_TYPE_FREE(newtype,ierror)

       deallocate(wk)

    end if ! color==1

    call MPI_BARRIER(MPI_COMM_WORLD, ierror)

max_stars_repo_name: ykawazura/calliope

max_issues_repo_name: ykawazura/calliope

max_forks_repo_name: ykawazura/calliope