subroutine zher2m (uplo,transa,transb,m,n,alpha,a,lda,b,ldb,beta,c,ldc)
!
! Author: Dieter Kvasnicka
!         Institute of Physical and Theoretical Chemistry
!         Vienna University of Technology, Austria
!         2000-2002
!
!  Purpose
!  =======
!
!  zher2m performs one of the matrix-matrix operations
!
!     op1(C) := alpha*op2( A )*op2( B ) + beta*op1(C),
!
!  where  op1( X ) means that only the UPPER or only the LOWER
!  part of X is accessed
!
!  and  op2( X ) is one of
!
!     op2( X ) = X   or   op2( X ) = X',
!
!  alpha and beta are scalars, and A, B and C are matrices, with
!  op2( A ) an m by n matrix,
!  op2( B ) an n by m matrix, and
!  C an m by m matrix.
!
!     uplo: 'U'pper or 'L'ower
!     transa,transb: op2( X ) = X   or   op2( X ) = X', with X=A,B
!     m: number of rows of op2(A) and op2(B) and of
!        rows and columns of C
!     n: number of columns of op2(A) and op2(B)
!     alpha
!     op2(A): mXn matrix
!     lda: leading dimension of a
!     op2(B): mXn matrix
!     ldb: leading dimension of b
!     beta
!     c: symmetric mXm matrix,
!        stored either above (uplo='U') or below (uplo='L') the
!        diagonal
!     ldc: leading dimension of c
!
  implicit none
  character*1 uplo,transa,transb
  integer m,n,lda,ldb,ldc
  complex*16 alpha,beta
  complex*16 a(lda,*), b(ldb,*), c(ldc,m)

  integer nblock
  parameter (nblock=128)
  complex*16 one,zero
  parameter (one=(1.0d0,0.0d0), zero=(0.0d0,0.0d0))

  integer j,ii,jj,imax,jmax
  logical upper,atransC,btransC,atransT,btransT,atransN,btransN

  logical lsame
  external lsame

  upper=lsame(uplo, 'U')
  atransC=lsame(transa, 'C')
  btransC=lsame(transb, 'C')
  atransT=lsame(transa, 'T')
  btransT=lsame(transb, 'T')
  atransN=lsame(transa, 'N')
  btransN=lsame(transb, 'N')

  if(.not.upper .and. atransN .and. (btransC.or.btransT)) then
     do ii=1,m,nblock
        imax=min(m,ii+nblock-1)
        do j=ii,imax
           call zgemv('n', imax-j+1, n, alpha, a(j,1), lda,b(j,1), ldb, beta, c(j,j), 1)
        end do
        if(imax.lt.m) call zgemm('n',transb,m-imax,imax-ii+1,n,alpha,a(imax+1,1),lda,b(ii,1),ldb,beta,c(imax+1,ii),ldc)
     end do
  else if(.not.upper .and. (atransC.or.atransT) .and. btransN) then
     do ii=1,m,nblock
        imax=min(m,ii+nblock-1)
        do j=ii,imax
           call zgemv(transa, n, imax-j+1, alpha, a(1,j), lda,b(1,j), 1, beta, c(j,j), 1)
        end do
        if(imax.lt.m) call zgemm(transa,'n',m-imax,imax-ii+1,n,alpha,a(1,imax+1),lda,b(1,ii),ldb,beta,c(imax+1,ii),ldc)
     end do
  else if(upper .and. (atransC.or.atransT) .and. btransN) then
     do jj=1,m,nblock
        jmax=min(m,jj+nblock-1)
        do j=jj,jmax
           call zgemv(transa, n, j-jj+1, alpha, a(1,jj), lda,b(1,j), 1, beta, c(jj,j), 1)
        end do
        if(jj.ne.1) call zgemm(transa,'n',jj-1,jmax-jj+1,n,alpha,a(1,1),lda,b(1,jj),ldb,beta,c(1,jj),ldc)
     end do
  else
     stop 'wrong parameter in zher2m'
  end if
end subroutine zher2m


text:
      program hwtest
      real*8 r1, r2 ,s
      r1 = 1.0
      r2 = 0.0
      s = hw1(r1, r2)
      write(*,*) 'hw1, result:',s
      write(*,*) 'hw2, result:'
      call hw2(r1, r2)
      call hw3(r1, r2, s)
      write(*,*) 'hw3, result:', s
      end

      real*8 function hw1(r1, r2)
      real*8 r1, r2
      hw1 = sin(r1 + r2)
      return
      end

      subroutine hw2(r1, r2)
      real*8 r1, r2, s
      s = sin(r1 + r2)
      write(*,1000) 'Hello, World! sin(',r1+r2,')=',s
 1000 format(A,F6.3,A,F8.6)
      return
      end

C     special version of hw1 where the result is
C     returned as an argument:

      subroutine hw3_v1(r1, r2, s)
      real*8 r1, r2, s
      s = sin(r1 + r2)
      return
      end

C     F2py treats s as input arg. in hw3_v1; fix this:

      subroutine hw3(r1, r2, s)
      real*8 r1, r2, s
Cf2py intent(out) s
      s = sin(r1 + r2)
      return
      end

C     test case sensitivity:

      subroutine Hw4(R1, R2, s)
      real*8 R1, r2, S
Cf2py intent(out) s
      s = SIN(r1 + r2)
      ReTuRn
      end

C end of F77 file

ext: f

max_stars_repo_name: sniemi/SamPy

max_issues_repo_name: sniemi/SamPy

max_forks_repo_name: sniemi/SamPy


text:
      double precision facscale
      common/facscale/facscale

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


! { dg-do run }
! A basic functional test of derived type extension.
!
! Contributed by Paul Thomas  <pault@gcc.gnu.org>
!
module persons
  type :: person
    character(24) :: name = ""
    integer :: ss = 1
  end type person
end module persons

module person_education
  use persons
  type, extends(person) :: education
    integer ::  attainment = 0
    character(24) :: institution = ""
  end type education
end module person_education

  use person_education
  type, extends(education) :: service
    integer :: personnel_number = 0
    character(24) :: department = ""
  end type service

  type, extends(service) :: person_record
    type (person_record), pointer :: supervisor => NULL ()
  end type person_record

  type(person_record), pointer :: recruit, supervisor

! Check that references by ultimate component work

  allocate (supervisor)
  supervisor%name = "Joe Honcho"
  supervisor%ss = 123455
  supervisor%attainment = 100
  supervisor%institution = "Celestial University"
  supervisor%personnel_number = 1
  supervisor%department = "Directorate"

  recruit => entry ("John Smith", 123456, 1, "Bog Hill High School", &
                    99, "Records", supervisor)

  if (trim (recruit%name) /= "John Smith") STOP 1
  if (recruit%name /= recruit%service%name) STOP 2
  if (recruit%supervisor%ss /= 123455) STOP 3
  if (recruit%supervisor%ss /= supervisor%person%ss) STOP 4

  deallocate (supervisor)
  deallocate (recruit)
contains
  function entry (name, ss, attainment, institution, &
                  personnel_number, department, supervisor) result (new_person)
    integer :: ss, attainment, personnel_number
    character (*) :: name, institution, department
    type (person_record), pointer :: supervisor, new_person

    allocate (new_person)

! Check mixtures of references
    new_person%person%name = name
    new_person%service%education%person%ss = ss
    new_person%service%attainment = attainment
    new_person%education%institution = institution
    new_person%personnel_number = personnel_number
    new_person%service%department = department
    new_person%supervisor => supervisor
  end function
end


text:
Davis nice climate provides little challenge for some introduced species.  While not nearly as problematic as quackgrass or the notorious zebra mussel, a thriving introduced species in Davis is the feral cat.  A feral cat is a descendant of a domesticated cat that has returned to the wild. It is distinguished from a stray cat, which is a pet cat that has been lost or abandoned, while feral cats are born in the wild; however, the offspring of a stray cat can be considered feral if born in the wild.

Life for feral cats can be notsonice. Dangers include: contracting viruses, attacks by other cats or other animals, and injuries caused by cars.

Controlling the Feral Cat Population: Spaying and Neutering

An important way for you to help control the feral cat population (and help both cats and native wildlife) is to have your own cat spayed or neutered. Low income pet owners can take advantage of the annual free Spay Day at UC Davis, or seek financial aid through the http://www.saacvoucher.org/saacvouchers/form.php Sacramento Area Animal Coalition.  The rest of you can see your Veterinary Hospitals vet.

The Feral Feline Organization is also dedicated to helping improve the lives of feral cats in Davis.  They conduct a Feral Feline Organization TNR program to keep the colonies from multiplying out of control.  If you are aware of feral cats in your neighborhood (or apartment complex) that have not been spayed or neutered, contact the Feral Feline Organization FFO or the Yolo County SPCA SPCA.

Generally, feral cats that have been spayed or neutered get their ears marked (so it doesnt happen a second time).  The markings can take the form of a snipped ear, a triangular notch, or even a holepunch.  A female would typically have her right ear marked, and a male would have his left ear marked.

For more info on other local wild animals see Town Wildlife, Cats or Lost Pets.

20050829 14:22:49 nbsp One vomited in my yard last weekend. Last year, one left a dead bird in my garden! Users/AlphaDog

20061206 06:39:38 nbsp A year ago a frightened 6 month old abandoned, and thus semiferal kitten showed up in my yard. My cats took him in and treated him with kindness, as I did. Its taken over a year, but I can finally pet himand pick him up. Granted, he takes me for a Bed and Breakfast: shows up for dinner, sleeps inside, has breakfast and leaves. But he is now family, and has been neutered, vaccinated, and will soon get a chip put in him. He keeps loosing collars, but Id rather them be a little loose than too tight. Give these guys a chance. Their almost all abandoned or runaways. Youll be surprised how grateful they are if you show them a little bit of respect and tenderness.  Users/MonicaBallyurban

20070823 12:05:45 nbsp there are some feral cats at parkside Users/StevenDaubert

20101020 11:22:54 nbsp Not only are there feral cats in my neighborhood, but also stray cats with collars.  Given that cats in the United States kill approximately 1,000,000 birds PER DAY, owners letting their cats run wild seems flagrantly inappopriate.  Additionally, the owner is reponsible if their cat kills any noninvasive bird species, and is subject to prosecution under both state and federal laws.  KEEP YOUR CATS INSIDE. Users/701Davis
 Um, what?!? Cats deserve freedom too! What cats and birds do together is their own business. My family had 2 outdoor cats at our old house for almost 20 years and they would only kill about 1 bird a year. Why should I believe 1000000 birds a day is very many? Do you have any sort of per capita statistics that would put it in context? Anyway, there doesnt seem to be any particular reason why birds are better to have around than cats in general, although I like seeing both out of doors very often. Users/NickSchmalenberger
       Cats and rodents > http://en.wikipedia.org/wiki/Toxoplasmosis Toxoplasmosis.  This is roughly where that old saying about pregnant women should stay away from cats/litter boxes comes from.  Cats are a primary host, and they typically get it from rodents they catch (see transmission).  I also think that what 701Davis is saying is that if the cats are killing that many birds, its messing with the natural food chain and the local ecosystem.  Those birds go on to eat other things (bugs, seeds, whatever) and ultimately are responsible for something or many things, right? Everything has its place.  More humans > more outdoor cats > throws off the balance, which I think is his/her summary.  And to respond to your below post, urban/feral cats seem to go hand in hand with human overpopulation.  More people, more outdoor cats, more cats breeding, and then too many wild cats.  Pretty much anyway you slice it, the natural balance is off.   users/edwins ES

20101020 12:14:25 nbsp As this article http://www.nytimes.com/2007/12/02/magazine/02catsvbirdst.html says, habitat destruction is the largest cause of bird population decline and cats, power poles and car traffic are secondary causes associated with the habitat destruction because they come with the people moving in. As with so many environmental problems, human overpopulation is the underlying cause. Without confronting that directly, not much will happen except sad impositions like keeping cats indoors. Users/NickSchmalenberger

20101106 11:08:50 nbsp I recently went trapping with the Yolo County SPCA at Sudwerk. We didnt get any cats, because most of them had already been spayed. But it was a great experience, anyone interested should definitely look into it!! Users/Churro615

20110402 13:21:57 nbsp The habitat destruction argument as an excuse not to do something about feral cats is a complete red herring.  Combined land use in the US for urban and rural residential lands is 8%.  Feral cats live only in a fraction of that 8% while habitat destruction occurs widely over the remaining 92%.  The impact on the bird may not be as great sources, but the impact BY the cat is enormous in small areas.  TNR programs are not effective at removing the impact cats have.  They do not reduce population sizes (the vacuum effect operates whether cats are trapped and removed or die naturally) and thus do not reduce the impact cats have.  Im all about animal welfare, Im just not about the biased feralcentric view of what that means.  We are fortunate to have a low incidence of Feline leukemia and Feline Immunodeficiency Virus in Davis, but its becoming a major problem in the foothills.  It wont be long before pets are catching it from ferals. Users/oob

20110514 08:52:09 nbsp Hi, Oob, I think most people would agree that something should be done about feral (I prefer to call them community) cats. The question is WHAT should be done? There are a couple of choices. 1. Do nothing. 2, We can trap and kill through the shelter. 3. We can shoot/drown/poison them. 4. Or we can trap and return them. I prefer the latter and Ill tell you why. TNR is more humane, more cost effective, doesnt cost taxpayers dollars (unlike trap and kill). Its also more effective because hundreds (thousands?) of people will volunteer their Saturday nights to trap cats that they know will survive after a spay surgery whereas very few people will spend their evenings trapping cats to take them to a shelter to get them killed. Budget strapped animal control agencies are no longer going out and rounding up cats. Trapped cats that are brought to AC are killed at taxpayer expense at approximately $50.00 to $100.00 per cat. Property owners and volunteers who choose TNR spend $10.00$20.00 per cat to get them spayed and vaccinated. Trap and Kill has not worked; indicated by the increasing number of cats at the shelters that we pay for. Finally, every year our cities consider whether they are going to contract with Animal Control, and its quite likely that some day soon, certain cities will NOT. The option to bring unowned/stray/homeless/found cats (and dogs) to the pound will not be an option. That limits your options to 1. do nothing, 2. Poison/shoot/drown or 3. TNR.

Davis has a comparatively small number of community cats. Estimates of cat numbers are that, nationwide, there is 1 community cat for every 6 humans. Davis is well below that number. And yes, while human populations have introduced cats, I believe increased numbers of free roaming cats are more a result of economics/demographics than population: more affluent communities tend to have fewer unowned cats. The rural and less affluent areas around Davis (Woodland farms/West Sac mobile home parks) have many more cats.

A managed colony has less impact on environment than a nonmanaged colony. Ive seen community cats catch rodents, never a bird (not saying it doesnt happen because Im sure it does). I have seen the birds eat the cat food :). TNRd cats who test positive for Felv/FIV are almost always euthanized. All others are vaccinated for both rabies and FVRCP. Spayed/neutered cats tend not to spread FIV because the two primary modes of transmission, fighting and sex, are reduced or eliminated. And all pet cat owners can and should vaccinate their cats to significantly reduce the chances of their cats getting these diseases. Better yet, keep them inside only.

If someone wants to do something about community cats, they can spend their free nights finding community cats, renting traps, setting them, sitting out in their car for hours, loading the cat into their car, holding it overnight at their house until animal control opens, loading it back into their car, calling in late for work, taking it to the shelter, loading the trap (which now has significant amounts of pee/poop in it) back into their car and doing again another night (as long as AC continues to take them). Or they can donate $15.00 and let other people do all that, understanding, of course, that those people probably love the cats and will spay, vaccinate, and return them.

Lastly, AC is considering more cost effective, humane methods of community cat control. Please support them when they consider these options. I think most of us would rather have our tax dollars spent supporting police/fire/schools/wildlife habitat preservation than killing cats.

 Users/SharonEM

20131222 08:16:55 nbsp Something paradoxical about the relationship between humans and feral cats in America(I say America because I have been to countries were feral cats are not much of an issue).  Cats were domesticated for or as a byproduct of them behaving in a way we would today in this country regard as a feral cat.  They roamed freely around farms, silos, barns, etc eating pests that would eat crops and stored food.  They were part of the pest management and would probably be given a treat now and then to reward their efforts and keep them around.  Regarding the cat as a helpless creature to be babied was usually something done only by upper classes.  There are still places on the world, especially those highly dependent on agriculture(I know because I visited one) were feral cats are a normal part of the towns landscape and pets cats are sometimes indistinguishable.  Thinking back on this I think there are exaggerated and legitimate reasons for controlling feral cat populations.  The more concerning things about feral cats would be disease spreading and competing with native wildlife, however many TNR activists(including the author of messybeast.com a popular cat site) vehemently deny cats kill millions of native wildlife.  At the same time however a lot of important questions are not being asked.  Such questions include what portion of the local bird or small animal population is being killed by cats, is this number concerning enough, how would one compare this to native predators, and what about feral cats in places full of nonnative wildlife anyway?(starlings, house sparrows, street rats, and pigeons(rock doves) are not native to America and among the most abundant urban wildlife).  I think there are a lot of gaps in understanding the full nature of this problem that is only muddied by people who vehemently see feral cats to be babied and protected by being fed and see them all like poor homeless orphans who will be all better once their gonads are removed or bird maniacs that see feral cats as the worst enemy and need to be killed.  When this politics occurs the right questions never get asked and instead implied to be assumed in one concrete way or another.  One of these is comparing the effectiveness of TNR vs. THR(TrapHysterectomyRelease), THR is meant to address unaltered feral cats moving into places were altered feral cats are, by instead of spaying and neutering to remove the uteruses of the females in a colony so the colony can resume defending their territory.  On the other hand it still poses a risk to pet cats that are let outdoors or run away.  On the other hand cats need to be handled more like dogs.  Most dog owners would never let their dogs roam the neighborhood!  Cats are seen as too aloof and independent to be trained but this is not true and cats can be leash trained and can benefit from catproof fencing.  I understand that keeping a car indoors can be a challenge but we seem to have a toolaid back approach to how pet cats are treated and disciplined in the first place.  Another thing:  some people claim that feral cats have no natural predators and thus flourish in numbers, obviously they ever heard of the threat coyotes pose to pet cats, or vehicles and what about those short miserable lives feral cats supposedly lead?(short and miserable compared to pet cats, but keep in mind that the exact same thing can be said when you compare the lifespan of a wild animal vs. the same wild animal in captivity.  Also if they are so miserable and starving what about those millions of birds killed?).   Users/KatherineLiu

20131222 08:23:31 nbsp I know there are typos there I am using a mobile device.   Users/KatherineLiu

20131222 08:44:03 nbsp Another thing Id like to add: feral cats first came to America as early as the first ships from Europe landed.  Feral cats on Australia are beginning to evolve and adapt to the climate the way dingoes have.  Makes me wonder if there are some feral cats that have been wild for many generations, much like feral horses and perhaps even become part of the ecosystem in some places.  Also a lot about the benefits of spaying and neutering is propaganda because nobody seems to talk about the bad side effects of them, or even the risk anesthesia poses especially done at a large and mass scale in spayneuter clinics.  Nobody mentions the cardiac, bone, and urinary issues that fixing increases.  Nobody mentions all the pets that dont survive their surgeries or dont wake up quite the same.  Nobody also mentions alternate methods of behavior and birth control on animals. Users/KatherineLiu

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      MODULE MOD_DLIOP

      CONTAINS

      SUBROUTINE DLIOP (I,ENTOTL,DOPAL,DETAL,VDOP,RSTAR,N,
     $                  EINST,WEIGHT,ELEVEL,LASTIND,INDLOW,INDNUP)
C*******************************************************************************
C***  DERIVATIVES OF ALL LINE OPACITIES AND EMISSIVITIES WITH RESPECT TO EN(I)
C***  AT CURRENT DEPTH POINT
C*******************************************************************************

      implicit real*8(a-h,o-z)

      DIMENSION EINST(N,N),WEIGHT(N),ELEVEL(N)
      DIMENSION INDLOW(LASTIND),INDNUP(LASTIND)
      DIMENSION DETAL(LASTIND),DOPAL(LASTIND)

C***  C3 = 4 * PI / H / C   ( CGS UNITS )
      DATA C3 /6.3268d16 /
C***  PI8 = 8*PI
      DATA PI8 /25.1327d0 /

C***  LOOP OVER ALL LINE TRANSITIONS  ----------------------------------
      DO 1 IND = 1, LASTIND

      LOW = INDLOW(IND)
      NUP = INDNUP(IND)

      IF (EINST(LOW, NUP) .EQ. -2.d0) GOTO 1

      DETAL(IND) = 0.0D0
      DOPAL(IND) = 0.0D0

      IF (I .NE. LOW .AND. I .NE. NUP) GOTO 1

      XLAMCM = 1.0D0 / (ELEVEL(NUP) - ELEVEL(LOW))

!     DND = DELTA-NUE-DOPPLER (HERTZ)
!     ********************************************
!     RINAT TAGIROV:
!     The description of VDOP is given in LIOP.for
!     ********************************************
      DND = VDOP * 1.0D5 / XLAMCM

      EMINDU=EINST(NUP,LOW)*XLAMCM*XLAMCM/PI8*RSTAR
      IF (I .EQ. LOW) THEN
C***  DERIVATIVE WITH RESPECT TO LOWER LEVEL
      DETAL(IND)=.0
      ABSORP=EMINDU*WEIGHT(NUP)/WEIGHT(LOW)
      DOPAL(IND)=ENTOTL*ABSORP/DND
      ELSE
C***  DERIVATIVE WITH RESPECT TO UPPER LEVEL
      EMSPON=EINST(NUP,LOW)/XLAMCM/C3*RSTAR
      DETAL(IND)=ENTOTL*EMSPON/DND
      DOPAL(IND)=-ENTOTL*EMINDU/DND
      ENDIF
    1 CONTINUE
C***  ENDLOOP  ---------------------------------------------------------

      RETURN

      END SUBROUTINE

      END MODULE

ext: for

max_stars_repo_name: rtagirov/nessy

max_issues_repo_name: rtagirov/nessy

max_forks_repo_name: rtagirov/nessy-src


text:
#define int INTEGER
       PROGRAM PREPRO
       int f
       END

max_stars_repo_name: eWert-Online/esy-cmake

max_issues_repo_name: eWert-Online/esy-cmake

max_forks_repo_name: eWert-Online/esy-cmake


text:
      LOGICAL FUNCTION DIGIT(STR)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    DIGIT
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1994-01-06
C
C ABSTRACT: THIS LOGICAL FUNCTION TESTS THE CHARACTERS IN A STRING TO
C   DETERMINE IF THEY ARE ALL DIGITS ('0','1','2','3','4','5','6','7',
C   '8' or '9').
C
C PROGRAM HISTORY LOG:
C 1994-01-06  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED
C                           DOCUMENTATION (INCLUDING HISTORY)
C 2007-01-19  J. ATOR    -- SIMPLIFIED LOGIC
C
C USAGE:    DIGIT (STR)
C   INPUT ARGUMENT LIST:
C     STR      - CHARACTER*(*): STRING
C
C   OUTPUT ARGUMENT LIST:
C     DIGIT    - LOGICAL: TRUE IF ALL CHARACTERS IN STR ARE DIGITS
C                ('0' - '9'), OTHERWISE FALSE
C
C REMARKS:
C    THIS ROUTINE CALLS:        None
C    THIS ROUTINE IS CALLED BY: CKTABA   NUMBCK   RDBFDX   RDUSDX
C                               Normally not called by any application
C                               programs but it could be.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      CHARACTER*(*) STR
      DIGIT = .FALSE.
      DO I=1,LEN(STR)
        IF( LLT(STR(1:1),'0') .OR. LGT(STR(1:1),'9') ) GOTO 100
      ENDDO
      DIGIT = .TRUE.

C  EXIT
C  ----

100   RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!----------------------------------------------------------------------------------
!:sdoc+:
!
! CSEM_LandMW_SfcOptics
!
! Module to compute the surface optical properties of Land surfaces at microwave
! frequencies, which is required for determining the surface contribution
! to the overall radiative transfer.
!
! This module provide a general template for developers to integrate individual
! models into the model repository, and a generic interface for
! the upper-level applications to access the model options in the model repository.
!
! This module implementes the FWD, TL and AD functions for variational data
! assimilation application and retrieval applications. It replacess the similar
! functions as those in the original CRTM surface module "CRTM_LandMW_SfcOptics",
! which was Written by Paul van Delst, 23-Jun-2005
!
!
! CREATION HISTORY:
!       Written by:     Ming Chen, 28-Feb-2015
!                       ming.chen@noaa.gov
!:sdoc-:
!-----------------------------------------------------------------------------------------------
!

MODULE CSEM_LandMW_SfcOptics

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE CSEM_Type_Kinds, ONLY: fp => CSEM_fp
  USE CSEM_Define
  USE CSEM_Model_Manager
  USE CSEM_Exception_Handler

  USE NESDIS_LandMW_PhyModel, ONLY: phy_iVar =>  iVar_type, &
                                    NESDIS_LandMW_Emiss,    &
                                    NESDIS_LandMW_Emiss_TL, &
                                    NESDIS_LandMW_Emiss_AD
  USE NESDIS_LandEM_Module,   ONLY: NESDIS_LandEM_213 =>  NESDIS_LandEM_213, &
                                    NESDIS_LandEM_205 =>  NESDIS_LandEM_OLD
  USE TELSEM_Atlas_Module
  USE TELSEM2_Atlas_Module
  USE CNRM_Atlas_Module

  !Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Data types
  ! Science routines
  PUBLIC :: CSEM_Compute_LandMW_SfcOptics
  PUBLIC :: CSEM_Compute_LandMW_SfcOptics_TL
  PUBLIC :: CSEM_Compute_LandMW_SfcOptics_AD
  PUBLIC :: iVar_type
  PUBLIC :: Get_Ref_Index

  REAL(fp), PARAMETER  ::  ONE = 1.0_fp, ZERO = 0.0_fp
  INTEGER,  PARAMETER  ::  MAX_N_ANGLES = 16
  ! Valid type indices for the microwave land emissivity model
  ! ...The soil types
  INTEGER, PARAMETER :: N_VALID_SOIL_TYPES = 8
  INTEGER, PARAMETER :: INVALID_SOIL    =  0
  INTEGER, PARAMETER :: COARSE          =  1
  INTEGER, PARAMETER :: MEDIUM          =  2
  INTEGER, PARAMETER :: FINE            =  3
  INTEGER, PARAMETER :: COARSE_MEDIUM   =  4
  INTEGER, PARAMETER :: COARSE_FINE     =  5
  INTEGER, PARAMETER :: MEDIUM_FINE     =  6
  INTEGER, PARAMETER :: COARSE_MED_FINE =  7
  INTEGER, PARAMETER :: ORGANIC         =  8
  ! ...The vegetation types
  INTEGER, PARAMETER :: N_VALID_VEGETATION_TYPES       = 12
  INTEGER, PARAMETER :: INVALID_VEGETATION             =  0
  INTEGER, PARAMETER :: BROADLEAF_EVERGREEN_TREES      =  1
  INTEGER, PARAMETER :: BROADLEAF_DECIDUOUS_TREES      =  2
  INTEGER, PARAMETER :: BROADLEAF_NEEDLELEAF_TREES     =  3
  INTEGER, PARAMETER :: NEEDLELEAF_EVERGREEN_TREES     =  4
  INTEGER, PARAMETER :: NEEDLELEAF_DECIDUOUS_TREES     =  5
  INTEGER, PARAMETER :: BROADLEAF_TREES_GROUNDCOVER    =  6
  INTEGER, PARAMETER :: GROUNDCOVER                    =  7
  INTEGER, PARAMETER :: GROADLEAF_SHRUBS_GROUNDCOVER   =  8
  INTEGER, PARAMETER :: BROADLEAF_SHRUBS_BARE_SOIL     =  9
  INTEGER, PARAMETER :: DWARF_TREES_SHRUBS_GROUNDCOVER = 10
  INTEGER, PARAMETER :: BARE_SOIL                      = 11
  INTEGER, PARAMETER :: CULTIVATIONS                   = 12
  ! --------------------------------------
  ! Structure definition to hold forward
  ! variables across FWD, TL, and AD calls
  ! --------------------------------------
  TYPE :: iVar_type
    PRIVATE
    TYPE(phy_iVar) :: phyVar(MAX_N_ANGLES)
    LOGICAL :: FREQ_CUTOFF = .FALSE.
  END TYPE iVar_type

  TYPE(CSEM_Model_ID), SAVE  :: MODEL

CONTAINS

!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CSEM_Compute_LandMW_SfcOptics
!
! PURPOSE:
!       Function to compute the surface emissivity and reflectivity at microwave
!       frequencies over a Land surface.
!
!       This function is a wrapper for third party code.
!
! PURPOSE:
!       Function to compute the surface emissivity  at microwave
!       frequencies over a Land surface.
!
!       This function encapsulates all available  CSEM LandMW emissivity models,and
!       provides a genereic interface for the upper-level user application.
!
! CALLING SEQUENCE:
!       IO_Status = CSEM_Compute_LandMW_SfcOptics(
!                                        Surface                  ,&  ! Input
!                                        SfcOptics                ,&  ! Output
!                                        Options                  ,&  ! OptionalInput
!                                        iVar)                     &  ! output
!
! INPUTS:
!
!       Surface:         CSEM Land Surface structure containing the surface state
!                        inputs.
!                        UNITS:      N/A
!                        TYPE:       CSEM_Land_Surface
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Options:         CSEM input extension structure containing additional inputs,
!                        e.g., Gelocation/Time Metadata,  Sensor observations, et al,
!                        which may be used by emprirical and semi-empirical models.
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       iVar:            Structure containing internal variables required for
!                        subsequent tangent-linear or adjoint model calls.
!                        The contents of this structure are NOT accessible
!                        outside of the CRTM_LandMW_SfcOptics module.
!                        UNITS:      N/A
!                        TYPE:       iVar_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(OUT)
!
!
! OUTPUTS:
!       SfcOptics:       CSEM output structure containing the surface optical
!                        properties required for the radiative transfer calculation,
!                        including wavelength and geometry structure as the inputs
!                        UNITS:      N/A
!                        TYPE:       CSEM_SfcOptics_Type
!                        DIMENSION:  Assumed-shape Array
!                        ATTRIBUTES: INTENT(IN OUT)
! FUNCTION RESULT:
!       IO_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the CSEM_Exception_Handler module.
!                        If == SUCCESS the computation was sucessful
!                           == FAILURE an unrecoverable error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!
! COMMENTS:
!                       CSEM infrastructure data structures are defined in CSEM_Define.f90
!
!
!:sdoc-:
!-----------------------------------------------------------------------------------------------

  FUNCTION CSEM_Compute_LandMW_SfcOptics( &
      Surface,                  &  ! Input
      SfcOptics,                &  ! Output
      Options,                  &  ! OptionalInput
      iVar)                     &  ! Output
    RESULT (IO_Status)

    TYPE(CSEM_Land_Surface),   INTENT(IN)     :: Surface
    TYPE(CSEM_Options_Type),   INTENT(IN)     :: Options
    TYPE(CSEM_SfcOptics_Type), INTENT(INOUT)  :: SfcOptics
    TYPE(iVar_type),           INTENT(OUT)    :: iVar

    ! Function result
    INTEGER :: IO_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CSEM_Compute_LandMW_SfcOptics'
    REAL(fp),     PARAMETER :: FREQUENCY_CUTOFF   = 200.0_fp  ! GHz
    REAL(fp),     PARAMETER :: DEFAULT_EMISSIVITY = 0.95_fp
    CHARACTER(LEN=256):: DB_PATH
    ! Local variables
    INTEGER :: i
    INTEGER :: imonth
    INTEGER :: N_Angles
    CHARACTER(256) :: msg

    IO_Status = SUCCESS
    N_Angles = SfcOptics%N_Angles
    SfcOptics%Reflectivity = ZERO
    SfcOptics%Emissivity   = ZERO

    MODEL = Inq_Model_Option("MW_LAND")
    ! use telsem_atlas_path in CSEM alg registor file
    DB_PATH = TRIM(MODEL%DATA_PATH)
    IF(TRIM(MODEL%NAME) .EQ. 'TELSEM_ATLAS' ) THEN

        imonth = Options%GeoInfo%Month
        IF (.NOT. (TELSEM_Atlas_Initialized(imonth))) THEN
          IO_Status = TELSEM_Atlas_Setup(imonth, path=TRIM(DB_PATH))
        ENDIF
        DO i = 1, N_Angles
          IO_Status = TELSEM_Atlas_Emiss(        &
                SfcOptics%Frequency,                &  ! Input
                SfcOptics%Angle(i),                 &  ! Input
                Options%GeoInfo%Latitude,           &  ! Input
                Options%GeoInfo%Longitude,          &  ! Input
                Options%GeoInfo%Month,              &  ! Input
                SfcOptics%Emissivity(i,2),          &  ! Output
                SfcOptics%Emissivity(i,1))             ! Output
                SfcOptics%Emissivity(i,3:4) = ZERO
        END DO

    ELSE IF(TRIM(MODEL%NAME) .EQ. 'TELSEM2_ATLAS' ) THEN

        imonth = Options%GeoInfo%Month
        IF (.NOT. (TELSEM_Atlas_Initialized(imonth))) THEN
          IO_Status = TELSEM2_Atlas_Setup(imonth, path=TRIM(DB_PATH))
        ENDIF
        DO i = 1, N_Angles
          IO_Status = TELSEM2_Atlas_Emiss(       &
                SfcOptics%Frequency,                &  ! Input
                SfcOptics%Angle(i),                 &  ! Input
                Options%GeoInfo%Latitude,           &  ! Input
                Options%GeoInfo%Longitude,          &  ! Input
                Options%GeoInfo%Month,              &  ! Input
                SfcOptics%Emissivity(i,2),          &  ! Output
                SfcOptics%Emissivity(i,1))             ! Output
                SfcOptics%Emissivity(i,3:4) = ZERO
        END DO
    ELSE IF(TRIM(MODEL%NAME) .EQ. 'CNRM_AMSUA_ATLAS' .AND. &
       INDEX(TRIM(Options%SensorObs%Sensor_Id),'amsu') > 0) THEN

        imonth = Options%GeoInfo%Month
        IF (.NOT. (CNRM_Atlas_Initialized(imonth))) THEN
          IO_Status = CNRM_Atlas_Setup(imonth, path=TRIM(DB_PATH))
        ENDIF

        DO i = 1, N_Angles
          IO_Status = CNRM_Atlas_Emiss(         &
                SfcOptics%Frequency,               &  ! Input
                SfcOptics%Angle(i),                &  ! Input
                Options%GeoInfo%Latitude,          &  ! Input
                Options%GeoInfo%Longitude,         &  ! Input
                Options%GeoInfo%Month,             &  ! Input
                SfcOptics%Emissivity(i,2),         &  ! Output
                SfcOptics%Emissivity(i,1))            ! Output
          SfcOptics%Emissivity(i,3:4) = ZERO
       END DO

    ELSE IF(TRIM(MODEL%NAME) .EQ. 'NESDIS_Land_MW213' ) THEN
       IF ( Surface%Soil_Type < 1 .OR. &
         Surface%Soil_Type > N_VALID_SOIL_TYPES ) THEN
         SfcOptics%Emissivity   = ZERO
         SfcOptics%Reflectivity = ZERO
         IO_Status = FAILURE
         msg = 'Invalid soil type index specified'
         WRITE(*,*) TRIM(msg);  RETURN
       END IF
       ! ...and the vegetation type
       IF ( Surface%Vegetation_Type < 1 .OR. &
         Surface%Vegetation_Type > N_VALID_VEGETATION_TYPES ) THEN
         SfcOptics%Emissivity   = ZERO
         SfcOptics%Reflectivity = ZERO
         IO_Status = FAILURE
         msg = 'Invalid vegetation type index specified'
         WRITE(*,*) TRIM(msg);  RETURN
       END IF

       IF(SfcOptics%Frequency < 80.0_fp) THEN

         DO i = 1, N_Angles
           CALL  NESDIS_LandEM_213(                &  ! Input
                 SfcOptics%Angle(i),               &  ! Input
                 SfcOptics%Frequency,              &  ! Input
                 Surface%Top_Soil_Moisture,        &  ! Input
                 Surface%Vegetation_Fraction,      &  ! Input
                 Surface%Top_Soil_Temperature,     &  ! Input
                 Surface%Land_Skin_Temperature,    &  ! Input
                 Surface%LAI,                      &  ! Input
                 Surface%Soil_Type,                &  ! Input
                 Surface%Vegetation_Type,          &  ! Input
                 SfcOptics%Emissivity(i,2),        &  ! Output
                 SfcOptics%Emissivity(i,1))            ! Output

           ! Assume specular surface
           SfcOptics%Reflectivity(i,1,i,1) = ONE-SfcOptics%Emissivity(i,1)
           SfcOptics%Reflectivity(i,2,i,2) = ONE-SfcOptics%Emissivity(i,2)
         END DO
       ELSE
         ! Frequency is too high for model. Use default.
         DO i = 1, SfcOptics%n_Angles
            SfcOptics%Emissivity(i,1:2)         = DEFAULT_EMISSIVITY
            SfcOptics%Reflectivity(i,1:2,i,1:2) = ONE-DEFAULT_EMISSIVITY
         END DO
       END IF

    ELSE IF(TRIM(MODEL%NAME) .EQ. 'NESDIS_Land_MW' ) THEN
       IF(SfcOptics%Frequency < 80.0_fp) THEN

         DO i = 1, N_Angles
           IO_Status = NESDIS_LandMW_Emiss(    &
                SfcOptics%Frequency,              &  ! Input
                SfcOptics%Angle(i),               &  ! Input
                Surface%Land_Skin_Temperature,    &  ! Input
                Surface%Top_Soil_Temperature,     &  ! Input
                Surface%Top_Soil_Moisture,        &  ! Input
                Surface%Vegetation_Fraction,      &  ! Input
                Surface%LAI,                      &  ! Input
                Surface%Vegetation_Type,          &  ! Input
                Surface%Soil_Type,                &  ! Input
                SfcOptics%Emissivity(i,2),        &  ! Output
                SfcOptics%Emissivity(i,1),        &  ! Output
                iVar%phyVar(i))                      ! Output
            ! Assume specular surface
            SfcOptics%Reflectivity(i,1,i,1) = ONE-SfcOptics%Emissivity(i,1)
            SfcOptics%Reflectivity(i,2,i,2) = ONE-SfcOptics%Emissivity(i,2)
        END DO
       ELSE
         ! Frequency is too high for model. Use default.
         ivar%FREQ_CUTOFF = .TRUE.
         DO i = 1, SfcOptics%n_Angles
            SfcOptics%Emissivity(i,1:2)         = DEFAULT_EMISSIVITY
            SfcOptics%Reflectivity(i,1:2,i,1:2) = ONE-DEFAULT_EMISSIVITY
         END DO

       ENDIF
    ELSE
       IO_Status = FAILURE
        msg = 'Invalid MW_LAND model name'
       WRITE(*,*) TRIM(msg);  RETURN

    END IF

    IO_Status = SUCCESS

  END FUNCTION CSEM_Compute_LandMW_SfcOptics

  SUBROUTINE Get_Ref_Index(Frequency, Polarization, i_ref_h,i_ref_v)

    REAL(fp):: Frequency(:)
    INTEGER :: Polarization(:)
    INTEGER :: i, i_ref_h,i_ref_v
    !INTEGER, PARAMETER :: N_POLARIZATION_TYPES    = 12
    !INTEGER, PARAMETER :: INVALID_POLARIZATION    = 0
    !INTEGER, PARAMETER :: UNPOLARIZED             = 1
    !INTEGER, PARAMETER :: INTENSITY               = UNPOLARIZED
    !INTEGER, PARAMETER :: FIRST_STOKES_COMPONENT  = UNPOLARIZED
    !INTEGER, PARAMETER :: SECOND_STOKES_COMPONENT = 2
    !INTEGER, PARAMETER :: THIRD_STOKES_COMPONENT  = 3
    !INTEGER, PARAMETER :: FOURTH_STOKES_COMPONENT = 4
    !INTEGER, PARAMETER :: VL_POLARIZATION         = 5
    !INTEGER, PARAMETER :: HL_POLARIZATION         = 6
    !INTEGER, PARAMETER :: plus45L_POLARIZATION    = 7
    !INTEGER, PARAMETER :: minus45L_POLARIZATION   = 8
    !INTEGER, PARAMETER :: VL_MIXED_POLARIZATION   = 9
    !INTEGER, PARAMETER :: HL_MIXED_POLARIZATION   = 10
    !INTEGER, PARAMETER :: RC_POLARIZATION         = 11
    !INTEGER, PARAMETER :: LC_POLARIZATION         = 12

    i_ref_h = -1 ; i_ref_v = -1
    DO i = 1, size(Frequency)
       IF(Frequency(i) > 15.0_fp .AND. Frequency(i) <= 20.0_fp ) THEN
           IF(Polarization(i) == 5 .OR. Polarization(i) == 9)  i_ref_v = i
           IF(Polarization(i) == 6 .OR. Polarization(i) == 10) i_ref_h = i
       ENDIF
    ENDDO

    IF(i_ref_v > 0 .AND. i_ref_h > 0)  RETURN
    i_ref_h = -1 ; i_ref_v = -1
    DO i = 1, size(Frequency)
       IF(Frequency(i) > 20.0_fp .AND. Frequency(i) <= 25.0_fp ) THEN
          i_ref_v = i
          i_ref_h = i
          EXIT
       ENDIF
    ENDDO

  END SUBROUTINE Get_Ref_Index

!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CSEM_Compute_LandMW_SfcOptics_TL
!
! PURPOSE:
!       Function to compute the tangent-linear surface emissivity and
!       reflectivity at microwave frequencies over a Land surface.
!
!       This function is a wrapper for third party code.
!
!
! CALLING SEQUENCE:
!       IO_Status = CSEM_Compute_LandMW_SfcOptics_TL(
!                                    Surface_TL,      &
!                                    SfcOptics_TL,    &
!                                    iVar)
!
! INPUTS:
!      Surface_TL:   Structure containing the tangent-linear surface
!                        state inputs required for the tangent-
!                        linear radiative transfer calculation.
!                        UNITS:      N/A
!                        TYPE:       CSEM_Land_Surface
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       iVar:            Structure containing internal variables required for
!                        subsequent tangent-linear or adjoint model calls.
!                        The contents of this structure are NOT accessible
!                        outside of the CRTM_LandMW_SfcOptics module.
!                        UNITS:      N/A
!                        TYPE:       iVar_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!      SfcOptics_TL: Structure containing the tangent-linear surface
!                        optical properties required for the tangent-
!                        linear radiative transfer calculation.
!                        UNITS:      N/A
!                        TYPE:       CSEM_SfcOptics_Type
!                        DIMENSION:  Assumed-shape Array
!                        ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
!       IO_Status:    The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the computation was sucessful
!                           == FAILURE an unrecoverable error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
! COMMENTS:
!
!:sdoc-:
!----------------------------------------------------------------------------------

  FUNCTION CSEM_Compute_LandMW_SfcOptics_TL( &
      Surface_TL,               &  ! Input
      SfcOptics_TL,             &  ! Output
      iVar        )             &  ! Internal variable input
   RESULT ( IO_Status  )
    ! Arguments
    TYPE(CSEM_Land_Surface),   INTENT(IN)     :: Surface_TL
    TYPE(CSEM_SfcOptics_Type), INTENT(INOUT)  :: SfcOptics_TL
    TYPE(iVar_type),           INTENT(IN)     :: iVar
    ! Function result
    INTEGER :: IO_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CSEM_Compute_LandMW_SfcOptics_TL'
    ! Local variables
    REAL(fp) :: Emissivity_H_TL, Emissivity_V_TL
    INTEGER :: i, N_Angles

    IO_Status = SUCCESS
    N_Angles = SfcOptics_TL%N_Angles

    SfcOptics_TL%Direct_Reflectivity = ZERO
    SfcOptics_TL%Reflectivity = ZERO
    SfcOptics_TL%Emissivity   = ZERO

    MODEL = Inq_Model_Option("MW_LAND")

    IF(TRIM(MODEL%NAME) .EQ. 'NESDIS_Land_MW' .AND. &
       .NOT. ivar%FREQ_CUTOFF) THEN
        DO i = 1, N_Angles
          IO_Status = NESDIS_LandMW_Emiss_TL(        &
                Surface_TL%Land_Skin_Temperature,       &  ! Input
                Surface_TL%Top_Soil_Temperature,        &  ! Input
                Surface_TL%Top_Soil_Moisture,           &  ! Input
                Surface_TL%Vegetation_Fraction,         &  ! Input
                Emissivity_H_TL,                        &  ! Output
                Emissivity_V_TL,                        &  ! Output
                iVar%phyVar(i))                            ! Internal variable output
          SfcOptics_TL%Emissivity(i,1)     = Emissivity_V_TL
          SfcOptics_TL%Emissivity(i,2)     = Emissivity_H_TL
          SfcOptics_TL%Reflectivity(i,1,i,1) = -SfcOptics_TL%Emissivity(i,1)
          SfcOptics_TL%Reflectivity(i,2,i,2) = -SfcOptics_TL%Emissivity(i,2)

        END DO

      ELSE
          SfcOptics_TL%Reflectivity = ZERO
          SfcOptics_TL%Emissivity   = ZERO
    END IF

  END FUNCTION CSEM_Compute_LandMW_SfcOptics_TL

!----------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!      CSEM_Compute_LandMW_SfcOptics_AD
!
! PURPOSE:
!       Function to compute the adjoint surface emissivity and
!       reflectivity at microwave frequencies over a land surface.
!
!       This function is a wrapper for third party code.
!
!
! CALLING SEQUENCE:
!       IO_Status = CSEM_Compute_LandMW_SfcOptics_AD(
!                                  SfcOptics_AD,         &
!                                  Surface_AD,           &
!                                  iVar)
!
! INPUTS:
!       SfcOptics_AD:    Structure containing the adjoint surface optical
!                        properties required for the adjoint radiative
!                        transfer calculation.
!                        *** COMPONENTS MODIFIED UPON OUTPUT ***
!                        UNITS:      N/A
!                        TYPE:       CSEM_SfcOptics_Type
!                        DIMENSION:  Assumed-shape Array
!                        ATTRIBUTES: INTENT(IN OUT)
!
!       iVar:            Structure containing internal variables required for
!                        subsequent tangent-linear or adjoint model calls.
!                        The contents of this structure are NOT accessible
!                        outside of the CRTM_LandMW_SfcOptics module.
!                        UNITS:      N/A
!                        TYPE:       iVar_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
! OUTPUTS:
!
!       Surface_AD:     Structure containing the adjoint surface state
!                        properties required for the adjoint radiative
!                        transfer calculation.
!                        *** COMPONENTS MODIFIED UPON OUTPUT ***
!                        UNITS:      N/A
!                        TYPE:       CSEM_Land_Surface
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!
! FUNCTION RESULT:
!       IO_Status:       The return value is an integer defining the error status.
!                        The error codes are defined in the Message_Handler module.
!                        If == SUCCESS the computation was sucessful
!                           == FAILURE an unrecoverable error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
! COMMENTS:
!       Note the INTENT on the input adjoint arguments are IN OUT regardless
!       of their specification as "input" or "output". This is because these
!       arguments may contain information on input, or need to be zeroed on
!       output (or both).
!
!:sdoc-:
!----------------------------------------------------------------------------------

  FUNCTION CSEM_Compute_LandMW_SfcOptics_AD( &
      SfcOptics_AD,             &  ! Input
      Surface_AD,               &  ! Output
      iVar        )             &  ! Internal variable input
   RESULT ( IO_Status  )
    ! Arguments
    TYPE(CSEM_Land_Surface),   INTENT(INOUT) :: Surface_AD
    TYPE(CSEM_SfcOptics_Type), INTENT(INOUT) :: SfcOptics_AD
    TYPE(iVar_type),           INTENT(IN)    :: iVar

    ! Function result
    INTEGER :: IO_Status
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CSEM_Compute_LandMW_SfcOptics_AD'
    ! Local variables

    ! Set up
    REAL(fp) :: Emissivity_H_AD, Emissivity_V_AD
    INTEGER :: i, N_Angles

    IO_Status = SUCCESS
    N_Angles = SfcOptics_AD%N_Angles

    MODEL = Inq_Model_Option("MW_LAND")
    IF(TRIM(MODEL%NAME) .EQ. 'NESDIS_Land_MW' .AND. &
        .NOT. ivar%FREQ_CUTOFF) THEN

        DO i = 1, N_Angles
           Emissivity_V_AD = SfcOptics_AD%Emissivity(i,1)
           Emissivity_H_AD = SfcOptics_AD%Emissivity(i,2)
           Emissivity_V_AD = Emissivity_V_AD - SfcOptics_AD%Reflectivity(i,1,i,1)
           Emissivity_H_AD = Emissivity_H_AD - SfcOptics_AD%Reflectivity(i,2,i,2)
           IO_Status = NESDIS_LandMW_Emiss_AD(       &
                Surface_AD%Land_Skin_Temperature,       &  ! Input
                Surface_AD%Top_Soil_Temperature,        &  ! Input
                Surface_AD%Top_Soil_Moisture,           &  ! Input
                Surface_AD%Vegetation_Fraction,         &  ! Input
                Emissivity_H_AD,                        &  ! Output
                Emissivity_V_AD,                        &  ! Output
                iVar%phyVar(i))                            ! Internal variable output

        END DO

      ELSE
        Surface_AD%Land_Skin_Temperature     =  Surface_AD%Land_Skin_Temperature     + ZERO
        Surface_AD%Top_Soil_Temperature      =  Surface_AD%Top_Soil_Temperature      + ZERO
        Surface_AD%Top_Soil_Moisture         =  Surface_AD%Top_Soil_Moisture         + ZERO
        Surface_AD%Vegetation_Fraction       =  Surface_AD%Vegetation_Fraction       + ZERO

    END IF

    SfcOptics_AD%Direct_Reflectivity = ZERO
    SfcOptics_AD%Reflectivity = ZERO
    SfcOptics_AD%Emissivity   = ZERO

  END FUNCTION CSEM_Compute_LandMW_SfcOptics_AD

END MODULE CSEM_LandMW_SfcOptics


text:
module mod_person
  type :: Person
    character(len=10) :: name
  contains
    procedure, pass(self) :: greet
  end type Person
contains
  impure elemental subroutine greet(self)
    class(Person), intent(in) :: self
    print *, 'Hello, my name is ' // trim(self % name) // '!'
  end subroutine greet
end module mod_person

program hello_derived_types
  use mod_person, only: Person
  implicit none
  type(Person) :: people(3) = [Person('Jill'), Person('James'), Person('Allison')]
  call people % greet()
end program hello_derived_types

max_stars_repo_name: wilsonify/c-consumer

max_issues_repo_name: wilsonify/c-consumer

max_forks_repo_name: wilsonify/c-consumer


text:
      subroutine tprn2(abeta,acphi,alamn2,amun3,in2pr,ipbtmx,
     $ jpfcmx,natmax,ncvn2,nerr,nn2pr,nn3tr,noutpt,npx2mx,npx2t,
     $ nttyo,nwarn,pcvn2,qpdn2,uaqsp,upair)
c
c     Test and process Pitzer data read from the data file that
c     pertain to a single neutral species (nn and nnn combinations,
c     such as CO2(aq)-CO2(aq) and CO2(aq)-CO2(aq)-CO2(aq). Find and
c     flag errors, such as duplication of data (e.g., two data
c     blocks for the same pair composed of a repeated neutral).
c     The conventional primitive Pitzer parameters are identical
c     to the observable parameters read from the data file:
c
c       lambda(nn) -> lambda(nn)
c       mu(nnn)    -> mu(nnn)
c
c     Check the coverage of entered Pitzer data against all possible
c     pairs that can be composed by repeating a neutral species present
c     on the data file.
c
c     This subroutine is called by:
c
c       EQPT/eqpt.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       nat    = the number of aqueous species
c       uaqsp  = array of names of aqueous species
c
c     Principal output:
c
c       nerr   = cumulative error counter
c       nwarn  = cumulative warning counter
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ipbtmx,jpfcmx,natmax,nn2pr,nn3tr,npx2mx
c
      integer noutpt,nttyo
c
      integer in2pr(nn2pr)
c
      integer ncvn2,nerr,npx2t,nwarn
c
      logical qpdn2(nn2pr)
c
      character(len=24) uaqsp(natmax),upair(2,npx2mx)
c
      real*8 abeta(jpfcmx,0:ipbtmx,npx2mx),acphi(jpfcmx,npx2mx),
     $ alamn2(jpfcmx,0:ipbtmx,nn2pr),amun3(jpfcmx,nn3tr)
c
      real*8 pcvn2
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,j,jp,jpair,j2,j3,j4,j5,n,ncount,ndupl,nlistl,
     $ nn,nodatc
c
      integer ilnobl
c
      character*56 ustr56
      character*24 unam1
      character*8 ux8
c
c-----------------------------------------------------------------------
c
c     Limit on the list of pairs for which no data were found.
c
      data nlistl / 20 /
c
c-----------------------------------------------------------------------
c
c     Initialize the data arrays.
c
      do n = 1,nn2pr
        qpdn2(n) = .false.
        do j = 1,jpfcmx
          do i = 0,ipbtmx
            alamn2(j,i,n) = 0.
          enddo
        enddo
      enddo
c
c     Note: nn3tr = nn2pr.
c
      do n = 1,nn3tr
        do j = 1,jpfcmx
          amun3(j,n) = 0.
        enddo
      enddo
c
c     Check the entered data for nn pairs.
c
      nodatc = 0
c
      do n = 1,nn2pr
        i = in2pr(n)
        unam1 = uaqsp(i)
c
c       Search for unam1 in the upair array.
c
        call srch22(jpair,unam1,unam1,upair,npx2mx,npx2t)
c
        if (jpair .gt. 0) then
c
c         Have found an entry.
c
          qpdn2(n) = .true.
c
c         Store the data.
c
          do j = 1,jpfcmx
            do i = 0,ipbtmx
              alamn2(j,i,n) = abeta(j,i,jpair)
            enddo
          enddo
c
          do j = 1,jpfcmx
            amun3(j,n)   = acphi(j,jpair)
          enddo
c
c         Check for duplicate data sets.
c
          ndupl = 0
          do jp = jpair + 1,npx2t
            if (unam1(1:24) .eq. upair(1,jp)(1:24)) then
              if (unam1(1:24) .eq. upair(2,jp)(1:24)) then
                ndupl = ndupl + 1
              endif
            endif
          enddo
c
          if (ndupl .gt. 0) then
            j2 = ilnobl(unam1)
            if (ndupl .eq. 1) then
              write (noutpt,1010) unam1(1:j2),unam1(1:j2)
              write (nttyo,1010) unam1(1:j2),unam1(1:j2)
 1010         format(/' * Error - (EQPT/tprn2) Have found a',
     $        ' duplicate data block on the DATA0 file',
     $        /7x,'for the nn pair ',a,', ',a,'.')
            else
              ux8 = ' '
              write (ux8,'(i5)') ndupl
              call lejust(ux8)
              j5 = ilnobl(ux8)
              write (noutpt,1020) ux8(1:j5),unam1(1:j2),unam1(1:j2)
              write (nttyo,1020) ux8(1:j5),unam1(1:j2),unam1(1:j2)
 1020         format(/' * Error - (EQPT/tprn2) Have found ',a,
     $        ' duplicate data blocks on the DATA0 file',
     $        /7x,'for the nn pair ',a,', ',a,'.')
            endif
            nerr = nerr + ndupl
          endif
c
        else
c
c         No data block was found on the DATA0 file.
c         Note that qpdn2(n) is left with a value of .false.
c
          nodatc = nodatc + 1
        endif
c
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (nodatc .gt. 0) then
        write (noutpt,1050)
        write (nttyo,1050)
 1050   format(/' * Warning - (EQPT/tprn2) Did not find a data',
     $  ' block on the DATA0 file',/7x,'for any of the following',
     $  ' nn pairs:',/)
c
        ncount = 0
        do n = 1,nn2pr
          if (.not.qpdn2(n)) then
            ncount = ncount + 1
            i = in2pr(n)
            unam1 = uaqsp(i)
            j2 = ilnobl(unam1)
            ustr56 = unam1(1:j2) // ', ' // unam1(1:j2)
            j4 = ilnobl(ustr56)
            write (noutpt,1060) ustr56(1:j4)
            write (nttyo,1060) ustr56(1:j4)
 1060       format(9x,a)
            if (ncount .eq. nlistl) go to 200
          endif
        enddo
  200   continue
c
        nn = nodatc - ncount
        if (nn .gt. 0) then
          write (ux8,'(i5)') nn
          call lejust(ux8)
          j3 = ilnobl(ux8)
          write (noutpt,1070) ux8(1:j3)
          write (nttyo,1070) ux8(1:j3)
 1070     format(/9x,'plus ',a,' others')
        endif
        write (noutpt,1080)
        write (nttyo,1080)
 1080   format(1x)
        nwarn = nwarn + 1
      endif
c
      ncvn2 = nn2pr - nodatc
      pcvn2 = (100.*ncvn2)/float(nn2pr)
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      end


text:
      MODULE cnc_C
      USE vast_kind_param, ONLY:  DOUBLE
!...Created by Pacific-Sierra Research 77to90  4.3E  06:33:54  12/28/06
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      REAL(DOUBLE), DIMENSION(2:5,2:4) :: CNC5C
      REAL(DOUBLE), DIMENSION(6,2:6) :: CNC6C
      END MODULE cnc_C

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


SUBROUTINE p41(input_file,output_file)
    !-------------------------------------------------------------------------
    ! Program 4.1 One dimensional analysis of axially loaded elastic rods
    !             using 2-node rod elements.
    !-------------------------------------------------------------------------
    USE main
    USE geom
    IMPLICIT NONE
    CHARACTER(len=60),INTENT(IN) :: input_file
    CHARACTER(len=60),INTENT(OUT) :: output_file
    INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
    INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,ndof=2,nels,neq,nlen,nod=2, &
    nodof=1,nn,nprops=1,np_types,nr
    REAL(iwp)::penalty=1.0e20_iwp,zero=0.0_iwp
    ! CHARACTER(LEN=15)::argv
    CHARACTER(LEN=30)::file_p41
    !-----------------------dynamic arrays------------------------------------
    INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),kdiag(:),nf(:,:),no(:),      &
    node(:),num(:)
    REAL(iwp),ALLOCATABLE::action(:),eld(:),ell(:),km(:,:),kv(:),loads(:),   &
    prop(:,:),value(:)
    !-----------------------input and initialisation--------------------------
    !CALL getname(argv,nlen)
    !OPEN(10,FILE=input_file)
    !OPEN(11,FILE=output_file)
    !OPEN(10, FILE='.//subordinates//Coupled-Problems//p41_1.dat')
    !OPEN(11, FILE='.//subordinates//Coupled-Problems//p41_1.res')
    OPEN(10, FILE=input_file)
    OPEN(11, FILE=output_file)
    READ(10,*)nels,np_types
    nn=nels+1
    ALLOCATE(g(ndof),num(nod),nf(nodof,nn),etype(nels),ell(nels),eld(ndof),  &
    km(ndof,ndof),action(ndof),g_g(ndof,nels),prop(nprops,np_types))
    READ(10,*)prop
    etype=1
    IF(np_types>1)READ(10,*)etype
    READ(10,*)ell
    nf=1
    READ(10,*)nr,(k,nf(:,k),i=1,nr)
    CALL formnf(nf)
    neq=MAXVAL(nf)
    ALLOCATE(kdiag(neq),loads(0:neq))
    kdiag=0
    !-----------------------loop the elements to find global arrays sizes-----
    elements_1: DO iel=1,nels
    num=(/iel,iel+1/)
    CALL num_to_g(num,nf,g)
    g_g(:,iel)=g
    CALL fkdiag(kdiag,g)
    END DO elements_1
    DO i=2,neq
    kdiag(i)=kdiag(i)+kdiag(i-1)
    END DO
    ALLOCATE(kv(kdiag(neq)))
    WRITE(11,'(2(A,I5))')                                                    &
    " There are",neq," equations and the skyline storage is",kdiag(neq)
    !-----------------------global stiffness matrix assembly------------------
    kv=zero
    elements_2: DO iel=1,nels
    CALL rod_km(km,prop(1,etype(iel)),ell(iel))
    g=g_g(:,iel)
    CALL fsparv(kv,km,g,kdiag)
    END DO elements_2
    !-----------------------read loads and/or displacements-------------------
    loads=zero
    READ(10,*)loaded_nodes,(k,loads(nf(:,k)),i=1,loaded_nodes)
    READ(10,*)fixed_freedoms
    IF(fixed_freedoms/=0)THEN
    ALLOCATE(node(fixed_freedoms),no(fixed_freedoms),value(fixed_freedoms))
    READ(10,*)(node(i),value(i),i=1,fixed_freedoms)
    DO i=1,fixed_freedoms
        no(i)=nf(1,node(i))
    END DO
    kv(kdiag(no))=kv(kdiag(no))+penalty
    loads(no)=kv(kdiag(no))*value
    END IF
    !-----------------------equation solution --------------------------------
    CALL sparin(kv,kdiag)
    CALL spabac(kv,loads,kdiag)
    loads(0)=zero
    WRITE(11,'(/A)')"  Node   Disp"
    DO k=1,nn
    WRITE(11,'(I5,2E12.4)')k,loads(nf(:,k))
    END DO
    !-----------------------retrieve element end actions----------------------
    WRITE(11,'(/A)')" Element Actions"
    elements_3: DO iel=1,nels
    CALL rod_km(km,prop(1,etype(iel)),ell(iel))
    g=g_g(:,iel)
    eld=loads(g)
    action=MATMUL(km,eld)
    WRITE(11,'(I5,2E12.4)')iel,action
    END DO elements_3
END SUBROUTINE


text:
program test
  real(kind=8), parameter :: g = 2.345

  print *, g, log(g)

end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE GG1112(IK1,IK2,BK1,BK2,ID1,ID2,BD1,BD2,K1,  &
                        QM1,QM2,QM3,QM4,WW)
!                                                                  *
! ----------------  SECTION METWO    SUBPROGRAM 15  -------------  *
!                                                                  *
!     THIS PACKAGE DETERMINES THE VALUES OF FOLLOWING MATRIX       *
!     ELEMENT:                                                     *
!                                                                  *
!       N1       (k1) (j1)(j2)  N1'       N2     (j2)   N2'     +- *
!     (j Q J ::[W(11)*A(1)]  ::j  Q'J')*(j Q J ::A(2)::j Q'J')  -+ *
!       1 1 1                   1  1 1    2 2 2         2 2 2   ++ *
!                                                               -- *
!     SUBROUTINE CALLED: C0T5S,RMEAJJ,WAP1                         *
!                                                                  *
!   Written by  G. Gaigalas                                        *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE CONS_C,          ONLY: ZERO, TENTH, HALF, EPS
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE mes_I
      USE c0t5s_I
      USE rmeajj_I
      USE wap1_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER,      INTENT(IN)               :: K1
      INTEGER,      INTENT(IN), DIMENSION(7) :: IK1, IK2, ID1, ID2
      REAL(DOUBLE), INTENT(IN)               :: QM1, QM2, QM3, QM4
      REAL(DOUBLE), INTENT(IN), DIMENSION(3) :: BK1, BK2, BD1, BD2
      REAL(DOUBLE), INTENT(OUT)              :: WW
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IQMM1, IQMM2, IQMM3, IQMM4, IQMM23, KK1
      REAL(DOUBLE) :: A1, S, BK, WA
!-----------------------------------------------
      WW=ZERO
      IF(IK1(3) > 9) THEN
        IF(IK1(4) > 2) CALL MES(34)
        IF(ID1(4) > 2) CALL MES(34)
      ENDIF
      IF(IK2(3) > 9) THEN
        IF(IK2(4) > 2) CALL MES(34)
        IF(ID2(4) > 2) CALL MES(34)
      ENDIF
      IQMM1=QM1+QM1+TENTH*QM1
      IQMM2=QM2+QM2+TENTH*QM2
      IQMM3=QM3+QM3+TENTH*QM3
      IQMM23=IQMM1+IQMM2+IQMM3
      IF(IK1(4) /= (ID1(4)+IQMM23))RETURN
      IQMM4=QM4+QM4+TENTH*QM4
      IF(IK2(4) /= (ID2(4)+IQMM4))RETURN
      KK1=K1*2
      CALL C0T5S(BD2(1),BD2(3),QM4,BK2(1),BK2(3),A1)
      IF(ABS(A1) < EPS)RETURN
!GG      CALL SJJ(IK2(3),IK2(1),IK2(7),IK2(6),ID2(1),ID2(7),ID2(6),S)
      CALL RMEAJJ(IK2(3),IK2(1),IK2(7),IK2(6),ID2(1),ID2(7),ID2(6),S)
      IF(ABS(S) < EPS)RETURN
      BK=HALF*DBLE(IK2(3))
      CALL WAP1(IK1,BK1,ID1,BD1,K1,BK,QM1,QM2,QM3,WA)
      IF(ABS(WA) < EPS)RETURN
      WW=-A1*WA*S/SQRT(DBLE(IK2(7)+1))
      RETURN
      END SUBROUTINE GG1112


text:
!  *********************************************************************
!  *                                                                   *
!  *                         subroutine print4                         *
!  *                                                                   *
!  *********************************************************************
!  Single Precision Version 1.2
!  Written by Gordon A. Fenton, TUNS, Sat Sep  9 13:21:15 1995
!  Latest Update: Oct 4, 2001
!
!  PURPOSE  prints four values according to a C style format
!
!  This routine takes as input a format, which is a character string
!  containing C style format descriptors, along with four additional
!  arguments and prints them to the given unit. For example, say you wish
!  to print the values of the variables 'i', `d', 'j', and 'f' to unit 6
!  so that the output appears as
!
!	d(3) = 16.49, f(12) = -0.03
!
!  (assuming i = 3, d = 16.49, j = 12, and f = -0.03), then this routine
!  would be called as follows;
!
!	call print4( 6, 'd(%i) = %f, f(%i) = %f%n', i, d, j, f )
!
!  Note that if the trailing '%n' characters are absent from the format string,
!  then the carriage return is inhibited in the output.
!
!  Arguments to this routine are as follows;
!
!       k	unit number to which output is directed. (input)
!
!     fmt	format string. (input)
!
!    val1	the first input argument. val1 can be either real or integer,
!		it's type is determined by the presence of a corresponding
!	 	%f or %e (for real) or %i (for integer) descriptor. See notes
!		below on the optional format specification which may follow
!		any of these format descriptors.
!		(input)
!
!    val2	see val1. (input)
!
!    val3	see val1. (input)
!
!    val4	see val1. (input)
!
!   The following % formats are recognized by this routine; only the %f, %e,
!   and %i formats actually consume an argument (ie, make use of one of
!   the val? arguments).
!
!	   %<space> to produce a space
!
!	   %%	to produce a % character
!
!	   %b	to produce a '\' character
!
!	   %t	to produce a tab character
!
!	   %n	to produce a newline character. This is NOT added to the
!		end of the output string by default. It must be explicitly
!		added (as in C).
!
!	   %f	to format a real value in the form '12.236'. If the value
!		is outside the range (1.e-6, 1.e+7) then the %e format is
!		automatically used. The %f format can be optionally written
!		in the following forms;
!
!		%w.df	specifies that a field width of `w' is to be used
!			with d digits to the right of the decimal place.
!			The decimal place takes up one character in the
!			total field width. `w' is actually a minimum width
!			and will be increased if more digits are required to
!			the left of the decimal point at the expense of
!			digits to the right of the decimal, if any. If the
!			number to be printed occupies less than w digits,
!			it is right justified. Digits to the right of the
!			7'th significant digit are left blank. For example,
!			%6.2f implies that the number 12.345 will be
!			printed as ' 12.35'.
!
!		%wf	specifies a minimum field width. The number of
!			decimal places is determined as necessary. Smaller
!			numbers are right justified. At most 7 significant
!			digits will be displayed. For example, %8f implies
!			that the number 12.34500 will be printed as
!			'  12.345'.
!
!		%.df	specifies a maximum number of digits to display
!			to the right of the decimal place.
!			Trailing zeroes are not shown. For example, %.2f
!			implies that the number 12.345 will be printed
!			as '12.35'.
!
!		%.f	at the moment, this is only specifies that zero is
!			printed as '0.'. The default is '0'.
!
!	   %e	to format a real value in scientific notion. For example
!		the number 12.345 is written as 1.2345e+01 in scientific
!		notation (where the e+01 means 'times 10 to the power 1').
!		The %e format can be optionally written in the following
!		forms, where;
!
!			y = 6 if the value being printed is non-negative
!			y = 7 if the value being printed is negative
!
!		%w.de	specifies that a field width of `w' is to be used
!			with d digits to the right of the decimal place.
!			w must be at least (d+y) to accomodate the leading
!			digit, the decimal, and the trailing 'e+dd' characters.
!			If w < (d+y) it is internally set to (d+y). If w
!			is greater than (d+y), the number is right justified
!			with leading blanks. For example, %10.2e implies that
!			the number 12.345 will be printed as '  1.23e+01'.
!
!		%we	specifies the field width only. The number of decimal
!			digits is determined as d = (w-y), so w must be at
!			least y (if it isn't, it is set to y). Since this is
!			a single precision routine, values of w greater than
!			12 result in a 12 character right justified number
!			(at most 7 significant digits are shown). For example,
!			%10e implies that the number 12.345 will be printed
!			as '1.2340e+01'.
!
!		%.de	specifies the number of digits to display to the right
!			of the decimal place. w is computed as (d+y). If d is
!			greater than 6, it is set internally to 6 to provide
!			at most 7 significant digits. For example, %.2e implies
!			that the number 12.345 will be printed as '1.23e+01'.
!
!	   %i	to format an integer value. Normally, the minimum field width
!		required to represent the number is used. However, the %i
!		format can be optionally written in the following form;
!
!		%wi	specifies that a field width of at least 'w' is used.
!			For example, %10i implies that the number 12 will be
!			printed as '        12'.
!
!  NOTES:
!   1)	if more than 4 occurrences of the formats %f, %e or %i appear in
!	the format string, the last argument value is repeatedly used (this may
!	lead to problems if it's type doesn't match the repeated format
!	descriptor).
!
!   2)	trailing white space in the format string is not printed in the output.
!	To actually get it in the output, it must be escaped, as in % % % ...
!
!  REVISION HISTORY
!  1.1	added %fw specification, where w is interpreted as the desired max.
!	number of significant digits to show on real output (Oct 9, 1995)
!  1.11	added comment about trailing blanks in fmt not being printed (10/14/96)
!  1.12	modified above writeup for new prints routines (Jun 9/97)
!  1.2	modified above writeup to reflect new formats and escapes (Oct 4/01)
!-------------------------------------------------------------------------
      subroutine print4( k, fmt, val1, val2, val3, val4 )
      parameter (MP = 4)
      character*(*) fmt
      real val1, val2, val3, val4, tmp(MP)

      tmp(1) = val1
      tmp(2) = val2
      tmp(3) = val3
      tmp(4) = val4

      call printv( k, fmt, tmp, MP )

      return
      end

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


!------------------------------------------------------------------------------
!
! `data` Source File
!
! data.f90 source file. This contains subroutines for initializing the
! positions, velocities, and masses for atoms in the simulation.
!
!------------------------------------------------------------------------------

MODULE data

  USE prms
  USE tools
  USE parallel

  IMPLICIT none

  real, allocatable, dimension(:,:)   :: X, X_old, X_new     ! atomic positions
  real, allocatable, dimension(:,:)   :: Xrand ! random atomic positions of the ions
  real, allocatable, dimension(:,:)   :: Xi    ! initial atomic positions
  real, allocatable, dimension(:,:)   :: V     ! atomic velocities
  real, allocatable, dimension(:)     :: mass,ken  ! atomic mass
  integer, allocatable, dimension(:)  :: atype,pp  ! atom type
  integer                             :: Nt0

CONTAINS

!This function initializes array storage for data stored for each atom
!X: position, V: velocities, mass:mass, atype:integer flag for type of atom
!Xrand: randomized starting ion positions
  SUBROUTINE initdata

    integer   :: i,l,j,jo
    integer   :: ix,iy,iz
    integer   :: Ncs

    real      :: ran1

    allocate (X(Natm,3),Xi(Natm,3),V(Natm,3),mass(Natm),atype(Natm),Xrand(Nrand,2), pp(Natm),ken(Natm))
    call ic

  END SUBROUTINE initdata

!determines which initial conditions to start the simulation with
!4 starts from scratch, 0 starts from a previously generated restart file
!-2 starts from a .dat file containing atomic positions and atom types
  SUBROUTINE ic

    integer   :: i,ierr, reason
    real      :: ran1
    character(100) :: cmd, fn

!for a from-scratch run, set timestep, time, and ion count to 0
    if (whichic.gt.0) then
        Nt0 = 0
        ions = 0
        restart_time = 0.0
    end if

    select case (whichic)
    case (4) ! added by Kallol, from-scratch run
        call initlat(X,Xi,atype)
        call readrandimp(Xrand)
        call assignmass
        if(amorcrys.eq.0) call randvel
        !call constrainmom(V,mass)

    case(0)
        !from a .restart file generated in a previous run
        !call initlat(X,Xi,atype)
        call readlat(Xi,atype)
        call MPI_BCAST(Xi,Natm*3,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
        if(startimpact .eq. 0) then
            open(res_unit,file='D/restart.in_old',form='UNFORMATTED',status='OLD')
        else
            write(cmd,"('ls D/restart.in',I6.6, '* > F/restart_list_',I3.3,'.dat')") startimpact, myid
            call system (cmd)
            write(fn,"('F/restart_list_',I3.3,'.dat')")myid
            open(33, file=fn)
            read(33,"(A28)",IOSTAT=reason)fn
            close (33)
            open(res_unit,file=fn,form='UNFORMATTED',status='OLD')
            if (myid .eq. 0) print*, "-----Restarting from ", startimpact, "th impact with restart file ", trim(fn)," -----"
        end if
        read(res_unit)Nt0,restart_time
        read(res_unit)X,V,mass,atype
        read(res_unit)ions
        close(res_unit)
        Xi(Nsg+1:Natm,:) = X(Nsg+1:Natm,:)
        V(Nsg+ions+1:Natm,:) = 0.0
        call readrandimp(Xrand)

    case(-2)
        call readlat(Xi,atype)
        call MPI_BCAST(Xi,Natm*3,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
        Nt0 = 4684000; ions = 400; restart_time = 0.3482100000E-09;
        open(1,file='D/atoms004684000.dat')
        do i = 1,Natm
            read(1,*)X(i,:),atype(i),pp(i),ken(i), V(i,:)
        end do
        call assignmass
        Xi(Nsg+1:Natm,:) = X(Nsg+1:Natm,:)
        V(Nsg+ions+1:Natm,:) = 0.0
        call readrandimp(Xrand)

    case default
        stop 'bad whichic'
    end select

  END SUBROUTINE ic

!initlat(X,Xi,atyp) reads in initial lattice conditions as X, assigns to Xi
  SUBROUTINE initlat(X,Xi,atype)
    integer, dimension(Natm) :: atype
    real, dimension(Natm,3)     :: X,Xi

    call readlat(X,atype)
    call initIons(X,atype)
    Xi = X

  END SUBROUTINE initlat

!readlat(X,atype) reads positions and types from si.dat to X and atype
  SUBROUTINE readlat(X,atype)
    integer, dimension(Natm) :: atype
    real, dimension(Natm,3)     :: X
    integer                  :: ierr
    integer                  :: i
    real                     :: ran1

    if (myid.eq.0) then
        if (amorcrys .eq. 0) then
            open(1,file='si.dat')
            do i = 1,Nsg
                read(1,*)X(i,:),atype(i)
            end do
        else
            open(1,file='5nm_relaxed_done.dat')
            do i = 1,Nsg
                read(1,*)X(i,:),atype(i), V(i,:)
            end do
        end if
        close(1)
    end if

  END SUBROUTINE readlat

!readrandimp(Xrand) reads randomized ion positions from fwhm.dat
  SUBROUTINE readrandimp(Xrand)
    real, dimension(Nrand,2)     :: Xrand
    integer                  :: ierr
    integer                  :: i
    real                     :: ran1

    if (myid.eq.0) then
        open(1,file='fwhm.dat')
        do i = 1,Nrand
            read(1,*)Xrand(i,:)
        end do
        close(1)
    end if
    call MPI_BCAST(Xrand,Nrand*2,MPI_REAL8,0,MPI_COMM_WORLD,ierr)

  END SUBROUTINE readrandimp

!randvel assigns initially random velocities to atoms, sets net momentum to zero
  SUBROUTINE randvel
    integer :: i, ierr
    real    :: ran1
    real, dimension(3)	:: mmtm

    !added by Kallol
    if (myid.eq.0) then
        mmtm = 0.0
        if (amorcrys .eq.0) then
            do i = 1,Nsg
                V(i,1) = Uo*(ran1(ranseed)-0.5)
                V(i,2) = Uo*(ran1(ranseed)-0.5)
                V(i,3) = Uo*(ran1(ranseed)-0.5)
                mmtm(:) = mmtm(:)+mass(i)*V(i,:)
            end do
      			mmtm = mmtm/Nsg

      			!subtract average momentum per particle from each particle (to make net zero)
      			do i = 1,Nsg
      				V(i,:) = V(i,:) - mmtm/mass(i)
      			end do

            do i = Nsg+1, Natm
                V(i,:) = 0.0
            end do

        else
            do i = Nsg+1, Natm
                V(i,:) = 0.0
            end do
        end if

    end if
    call MPI_BCAST(V,Natm*3,MPI_REAL8,0,MPI_COMM_WORLD,ierr)

  END SUBROUTINE randvel

!assignmass initializes mass to each atom, according to its atype value
  SUBROUTINE assignmass
    integer :: i

    do i = 1,Natm
        if (atype(i).eq.1) then
            mass(i) = massSi
        else if (atype(i).eq.2) then
            mass(i) = massGe
        else if (atype(i).eq.3) then
            mass(i) = massGa
        end if
    end do

  END SUBROUTINE assignmass

!initIons(X,atype) initializes ion positions in the pre-firing layer, into X
  SUBROUTINE initIons(X,atype)
    integer ,dimension(Natm)	:: atype
    integer                    	:: ii,i
    integer                     :: j, ierr
    real,dimension(Natm,3)      :: X
    real,dimension(3)           :: xxij
    real                        :: ran1
    real                        :: rij

    if(myid .eq. 0) then
        do i=Nsg+1,Natm
            X(i,1) = Lb(1)*ran1(ranseed)
            X(i,2) = Lb(2)*ran1(ranseed)
            X(i,3) = ionz*Lb(3)/10.0
            !V(i,:) = 0.
            atype(i) = 3
        end do

        !  Check Ion positions and reassign if necessary
        do i=Nsg+1,Natm
            do j=Nsg+1,Natm
                if (i .gt. j) then
79                  xxij=X(i,:)-X(j,:)
                    rij = sqrt(xxij(1)**2 + xxij(2)**2)
                    if (rij .lt. sigmaSi) then
                        X(i,1) = Lb(1)*ran1(ranseed)
                        X(i,2) = Lb(2)*ran1(ranseed)
                        goto 79
                    end if
                end if
            end do
        end do
    end if

    call MPI_BCAST(X,Natm*3,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(atype,Natm,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

  END SUBROUTINE initIons

!initranions creates fwhm.dat, a random 2D gaussian distribution for ion impact locations
  SUBROUTINE initranions
    real          :: rion1, rion2, rion3, rion4, ran1
    integer       :: i
    character(30) :: fn

    write(fn,"('fwhm.dat')")
    open(ran_unit,file=fn)

    do i=1,Nrand
      rion1 = ran1(ranseed)
      rion2 = ran1(ranseed)

      !convert to a normal distribution
      !also scale to a fwhm given in nm
      rion3 = sqrt(-2*LOG(rion1))*COS(2*Pi*rion2)*fwhm/2.35*1.E-9
      rion4 = sqrt(-2*LOG(rion2))*SIN(2*Pi*rion1)*fwhm/2.35*1.E-9

      write(ran_unit,"(E12.6,A1,E12.6)")rion3,' ',rion4
    end do

    close(ran_unit)

  END SUBROUTINE initranions

END MODULE data


text:
module mod_combine
    use mod_types
    contains
        subroutine combine_params(params, np, indp, &
                constants, nc, indc, allparams)

            ! Inputs
            integer(kind=i2), intent(in) :: np, nc, indp(np), indc(nc)
            real(kind=r2), intent(in) :: params(np), constants(nc)

            ! Internals
            real(kind=r2) :: allparams(np+nc)

            allparams(indp) = params
            if(nc > 0) allparams(indc) = constants
            return
        end subroutine
end module

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_licenses: NCSA Unlicense