text:
      subroutine rburn(plant)
c
c     + + + PURPOSE + + +
c     Updates range vegetation attributes as the result of burning.
c     (Burning may only be done once a year.)
c
c     Called from RANGE
c     Author(s): Weltz
c     Reference in User Guide: Chapter 8
c
c     Changes:
c           1) Dereferenced the following common blocks:
c              CRPVR3, CRPVR5, COVER, CRPOUT, CRPPRM, ENDS,
c              PARAME, RINPT2, RINPT6, UPDATE.
c           2) Dereferenced: PTILTY.INC & PMXPND.INC.
c           3) Changed:
c                nowres = 1
c                rmogt(nowres,iplane) = rmogt(nowres,iplane) * reduce(iplane)
c              to:
c                rmogt(1,iplane) = rmogt(1,iplane) * reduce(iplane)
c           4) Removed generic SAVE of all local variables.
c           5) Changed PLIVE(PLANT) to PLIVE(PLANT,IPLANE). -- CRM -- 5/14/93.
c           6) PROOT was altered analogous to PLIVE. -- CRM -- 5/14/93.
c
c     Version: This module recoded from WEPP Version 92.25.
c     Date recoded: 01/28/93.
c     Recoded by: Charles R. Meyer.
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
      include 'pntype.inc'
      include 'pmxelm.inc'
      include 'pmxnsl.inc'
      include 'pmxpln.inc'
      include 'pmxres.inc'
      include 'pmxgrz.inc'
      include 'pmxcrp.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer plant
c
c     + + + ARGUMENT DEFINITIONS + + +
c     plant  - current plant type; set to itype(nowcrp,iplane) in RANGE
c
c     + + + COMMON BLOCKS + + +
      include 'ccrpvr1.inc'
c      modify: rmogt(mxres,mxplan), rmagt(mxplan)
c
      include 'ccrpvr2.inc'
c       write: vdmt(mxplan)
c
      include 'chydrol.inc'
c        read: rain(mxplan)
c
      include 'crinpt1.inc'
c      modify: plive(ntype,mxplan), proot(ntype,iplane)
c       write: access(mxplan)
c
      include 'crinpt3.inc'
c        read: alter(mxplan),burned(mxplan),change(mxplan),
c              hurt(mxplan),reduce(mxplan)
c      modify: jfdate(mxplan),shgt(ntype),spop(ntype)
c
      include 'crinpt4.inc'
c      modify: xlive(mxplan), dwood(mxplan)
c
      include 'crinpt5.inc'
c      modify: tpop(ntype),thgt(ntype)
c
      include 'crout.inc'
c      modify: tlive(mxplan)
c
      include 'cstruc.inc'
c        read: iplane
c
      include 'cclim.inc'
c        read: am2(mxplan)
c
c     + + + END SPECIFICATIONS + + +
c
c
c      *** L0 IF ***
c      If rainfall is <= 7.5 mm, AND five day antecedent moisture for
c      submerged residue <= 25 mm, proceed with burning.
c
c XXX NOTE - may want to change "rain(iplane)" to "rain+melt+irdept"
c XXX        so that impact of snow melt and irrigation can also be
c XXX        accounted for.   dcf  12/19/94
c
      if (rain(iplane).le.0.0075.and.am2(iplane).le.0.025) then
c
c       *** L1 IF ***
c       ------ (EVERGREENS PRESENT)
        if (xlive(iplane).gt.0.0) then
c         -------- (MIXED POPULATION)
          if (tlive(iplane).gt.xlive(iplane)) then
c           **************************************************************
c           * Note: The amount of live material is reduced by the same   *
c           *       fractions as the dead material.  The accuracy to be  *
c           *       gained, did not seem worth the complexity it would   *
c           *       cost to handle them separately. -- CRM -- 3/17/93.   *
c           **************************************************************
c           -----------total live plant material
c           (WEPP Equation 8.6.14)
            tlive(iplane) = ((tlive(iplane)-xlive(iplane))*
     1          reduce(iplane)) + (xlive(iplane)*hurt(iplane))
c         -------- (ONLY EVERGREENS)
          else
            xlive(iplane) = xlive(iplane) * hurt(iplane)
c           (WEPP Equation 8.6.13)
            tlive(iplane) = tlive(iplane) * hurt(iplane)
          end if
c
c       *** L1 ELSE ***
c       ------ (NO EVERGREENS PRESENT)
        else
          tlive(iplane) = tlive(iplane) * reduce(iplane)
c       *** L1 ENDIF ***
        end if
c
c       ------ update on-ground residue by non-evergreen reduction fraction
c       **************************************************************
c       * Note: If present, evergreen residue is reduced by the same *
c       *       fraction as non-evergreen.  The accuracy to be       *
c       *       gained, did not seem worth the complexity required   *
c       *       to handle them separately. -- CRM -- 3/17/93.        *
c       **************************************************************
c       -- XXX -- We are assuming evergreen residue (if present) is reduced
c       at the same rate as non-evergreens. -- CRM 1/28/93.
c       (WEPP Equation 8.6.16)
        rmogt(1,iplane) = rmogt(1,iplane) * reduce(iplane)
c       ------ update standing herbaceous residue by non-evergreen reduction fraction
c       (WEPP Equation 8.6.15)
        rmagt(iplane) = rmagt(iplane) * reduce(iplane)
c       ------ standing woody biomass
c       (WEPP Equation 8.6.10)
        dwood(iplane) = dwood(iplane) * burned(iplane)
c       ------ update fraction of forage available for consumption, to coefficient
c       of increase in availability
c       ------ replace the fraction of forage available, with fraction available
c       after burning
        access(iplane) = alter(iplane)
c       ------ update maximum [potential] standing live biomass, using fraction
c       change in above and below ground biomass
c       (WEPP Equation 8.6.13)
        plive(plant,iplane) = plive(plant,iplane) * change(iplane)
c       ------ update peak root biomass, using fraction change in above and
c       below ground biomass
c       (WEPP Equation 8.6.12)
        proot(plant,iplane) = proot(plant,iplane) * change(iplane)
c       ------ update total above ground dry weight (VDMT), using updated:
c       total live plant material on day of simulation (TLIVE),
c       surface residue mass above ground today (RMAGT), and standing
c       dead biomass left after burning or herbicide (DWOOD).
        vdmt(iplane) = tlive(iplane) + rmagt(iplane) + dwood(iplane)
c       **************************************************************
c       * Note: The average shrub & tree heights, and the average    *
c       *       shrub & tree populations are reduced by the factor   *
c       *       for reduction of woody biomass (BURNED).  It was     *
c       *       decided that this simplification would give adequate *
c       *       accuracy.  -- CRM -- 3/17/93.                        *
c       **************************************************************
c       ------ update average shrub height, using reduction in standing dead
c       wood after burning
c       (WEPP Equation 8.6.11)
        shgt(plant) = shgt(plant) * burned(iplane)
c       ------ average tree height
        thgt(plant) = thgt(plant) * burned(iplane)
c       ------ average number of shrubs along a 100 m transect
        spop(plant) = spop(plant) * burned(iplane)
c       ------ average number of trees along a 100 m transect
        tpop(plant) = tpop(plant) * burned(iplane)
c
c     *** L0 ELSE ***
c     If rainfall is > 7.5 mm, or five day antecedent moisture for
c     submerged residue > 25 mm, set day of burning 1 day later; ie,
c     delay burning 1 day.
c
      else
        jfdate(iplane) = jfdate(iplane) + 1
c
c     *** L0 ENDIF ***
      end if
c
      return
      end

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
module funcs_mod
implicit none
private
public :: wp, f_real_real, square, cube
integer, parameter :: wp = kind(1.0d0)
!
interface ! defines a function signature
   pure function f_real_real(x) result(y)
   import wp
   implicit none
   real(kind=wp), intent(in) :: x
   real(kind=wp)             :: y
   end function f_real_real
end interface
!
contains
!
pure function square(x) result(y)
real(kind=wp), intent(in) :: x
real(kind=wp)             :: y
y = x**2
end function square
!
pure function cube(x) result(y)
real(kind=wp), intent(in) :: x
real(kind=wp)             :: y
y = x**3
end function cube
!
end module funcs_mod

module integral_mod
use funcs_mod, only: wp, f_real_real
implicit none
private
public :: integral
contains
!
pure function integral(f,a,b,n) result(y)
procedure(f_real_real)    :: f
! f can be any function with f_real_real signature
real(kind=wp), intent(in) :: a,b
integer      , intent(in) :: n
real(kind=wp)             :: y
integer                   :: i
real(kind=wp)             :: x,h,fsum
h = (b-a)/(n-1)
x = a
fsum = 0.0
do i=1,n
   if (i == 1 .or. i == n) then
      fsum = fsum + f(x)/2
   else
      fsum = fsum + f(x)
   end if
   x = x + h
end do
y = h*fsum
end function integral
!
end module integral_mod
!
program test_procedure_arg
use integral_mod, only: integral
use funcs_mod   , only: wp, square, cube
implicit none
real(kind=wp) :: a,b
integer       :: n
n = 10**6
a = 0.0_wp
b = 10.0_wp
print "(*(f12.6))", integral(square,a,b,n),integral(cube,a,b,n)
end program test_procedure_arg
! output:
!  333.333333 2500.000000

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
*=======================================================================
*
* WCSLIB 4.24 - an implementation of the FITS WCS standard.
* Copyright (C) 1995-2014, Mark Calabretta
*
* This file is part of WCSLIB.
*
* WCSLIB is free software: you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* WCSLIB is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with WCSLIB.  If not, see http://www.gnu.org/licenses.
*
* Direct correspondence concerning WCSLIB to mark@calabretta.id.au
*
* Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
* http://www.atnf.csiro.au/people/Mark.Calabretta
* $Id: tspc.f,v 4.24 2014/09/18 15:25:02 mcalabre Exp $
*=======================================================================

      PROGRAM TSPC
*-----------------------------------------------------------------------
*
* TSPC tests the spectral transformation driver routines for closure.
*
*-----------------------------------------------------------------------
*     Maximum length of spectral axis - see CLOSURE.
      INTEGER   NSPEC
      PARAMETER (NSPEC = 10001)

      INTEGER   CLOSURE, NAXISJ, NFAIL
      DOUBLE PRECISION C, CDELTX, CRPIXJ, CRVALX, MARS(0:6), RESTFRQ,
     :          RESTWAV, X1, X2

      COMMON /SPECTRO/ MARS

      DATA C /2.99792458D8/

*     KPNO MARS spectrograph grism parameters.
      DATA MARS /4.5D5, 1D0, 27D0, 1.765D0, -1.077D6, 3D0, 5D0/
*-----------------------------------------------------------------------
      WRITE (*, 10)
 10   FORMAT ('Testing closure of WCSLIB spectral transformation ',
     :        'routines (tspc.f)',/,
     :        '--------------------------------------------------',
     :        '-----------------')

      NFAIL = 0

*     PGPLOT initialization.
      CALL PGBEG (0, '/xwindow', 1, 1)

      NAXISJ = NSPEC
      CRPIXJ = NAXISJ/2 + 1

      RESTFRQ = 1420.40595D6
      RESTWAV = C/RESTFRQ
      X1 = 1D9
      X2 = 2D9
      CDELTX = (X2 - X1)/(NAXISJ - 1)
      CRVALX = X1 + (CRPIXJ - 1.0)*CDELTX
      WRITE (*, 20) X1*1D-9, X2*1D-9, CDELTX*1D-3
 20   FORMAT (/,'Linear frequency axis, span:',F4.1,' to',F4.1,
     :        ' (GHz), step:',F8.3,' (kHz)',/,'---------------------',
     :        '-----------------------------------------------------')
      NFAIL = NFAIL + CLOSURE('WAVE-F2W', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VOPT-F2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ZOPT-F2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AWAV-F2A', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VELO-F2V', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('BETA-F2V', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)

      RESTWAV = 700D-9
      RESTFRQ = C/RESTWAV
      X1 = 300D-9
      X2 = 900D-9
      CDELTX = (X2 - X1)/(NAXISJ - 1)
      CRVALX = X1 + (CRPIXJ - 1D0)*CDELTX
      WRITE (*, 30) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9
 30   FORMAT (/,'Linear vacuum wavelength axis, span:',I4,' to',I4,
     :        ' (nm), step:',F9.6,' (nm)',/,'----------------------',
     :        '----------------------------------------------------')
      NFAIL = NFAIL + CLOSURE('FREQ-W2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AFRQ-W2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ENER-W2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVN-W2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VRAD-W2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AWAV-W2A', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VELO-W2V', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('BETA-W2V', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)

      WRITE (*, 40) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9
 40   FORMAT (/,'Linear air wavelength axis, span:',I4,' to',I4,
     :        ' (nm), step:',F9.6,' (nm)',/,'----------------------',
     :        '----------------------------------------------------')
      NFAIL = NFAIL + CLOSURE('FREQ-A2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AFRQ-A2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ENER-A2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVN-A2F', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VRAD-A2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVE-A2W', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VOPT-A2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ZOPT-A2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VELO-A2V', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('BETA-A2V', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)

      RESTFRQ = 1420.40595D6
      RESTWAV = C/RESTFRQ
      X1 = -0.96D0*C
      X2 =  0.96D0*C
      CDELTX = (X2 - X1)/(NAXISJ - 1)
      CRVALX = X1 + (CRPIXJ - 1D0)*CDELTX
      WRITE (*, 50) INT(X1), INT(X2), INT(CDELTX)
 50   FORMAT (/,'Linear velocity axis, span:',I11,' to',I10,
     :        ' m/s, step:',I6,' (m/s)',/,'-----------------------',
     :        '---------------------------------------------------')
      NFAIL = NFAIL + CLOSURE('FREQ-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AFRQ-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ENER-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVN-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VRAD-V2F', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVE-V2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VOPT-V2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ZOPT-V2W', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AWAV-V2A', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)

      RESTWAV = 650D-9
      RESTFRQ = C/RESTWAV
      X1 =  300D-9
      X2 = 1000D-9
      CDELTX = (X2 - X1)/(NAXISJ - 1)
      CRVALX = X1 + (CRPIXJ - 1D0)*CDELTX
      WRITE (*, 60) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9
 60   FORMAT (/,'Vacuum wavelength grism axis, span:',I4,' to',I5,
     :        ' (nm), step:',F9.6,' (nm)',/,'----------------------',
     :        '----------------------------------------------------')
      NFAIL = NFAIL + CLOSURE('FREQ-GRI', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AFRQ-GRI', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ENER-GRI', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVN-GRI', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VRAD-GRI', RESTFRQ, 0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('WAVE-GRI', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VOPT-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('ZOPT-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('AWAV-GRI', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VELO-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('BETA-GRI', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)

*     Reproduce Fig. 5 of Paper III.
      NAXISJ = 1700
      CRPIXJ = 719.8D0
      CRVALX = 7245.2D-10
      CDELTX = 2.956D-10
      RESTWAV = 8500D-10
      RESTFRQ = C/RESTWAV
      X1 = CRVALX + (1 - CRPIXJ)*CDELTX
      X2 = CRVALX + (NAXISJ - CRPIXJ)*CDELTX
      MARS(5) = 0D0
      MARS(6) = 0D0
      WRITE (*, 70) INT(X1*1D9), INT(X2*1D9), CDELTX*1D9
 70   FORMAT (/,'Air wavelength grism axis, span:',I4,' to',I5,
     :        ' (nm), step:',F9.6,' (nm)',/,'----------------------',
     :        '----------------------------------------------------')
      NFAIL = NFAIL + CLOSURE('AWAV-GRA', 0D0,     0D0, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)
      NFAIL = NFAIL + CLOSURE('VELO-GRA', 0D0, RESTWAV, NAXISJ, CRPIXJ,
     :                        CDELTX, CRVALX)

      CALL PGASK(0)
      CALL PGEND()

      IF (NFAIL.NE.0) THEN
        WRITE (*, 80) NFAIL
 80     FORMAT (/,'FAIL:',I5,' closure residuals exceed reporting ',
     :    'tolerance.')
      ELSE
        WRITE (*, 90)
 90     FORMAT (/,'PASS: All closure residuals are within reporting ',
     :    'tolerance.')
      END IF

      END

*=======================================================================

      INTEGER FUNCTION CLOSURE (CTYPES, RESTFRQ, RESTWAV, NAXISJ,
     :  CRPIXJ, CDELTX, CRVALX)

      INTEGER   NSPEC
      PARAMETER (NSPEC = 10001)

      INTEGER   J, NAXISJ, NFAIL, RESTREQ, STAT1(NSPEC), STAT2(NSPEC),
     :          STATUS
      REAL      TMP, X(NSPEC), XMIN, XMAX, Y(NSPEC), YMAX, YMIN
      DOUBLE PRECISION CDELTS, CDELTX, CLOS(NSPEC), CRPIXJ, CRVALS,
     :          CRVALX, DSDX, MARS(0:6), RESID, RESIDMAX, RESTFRQ,
     :          RESTWAV, SPEC1(NSPEC), SPEC2(NSPEC), TOL
      CHARACTER CTYPES*8, PTYPE, SCODE*3, SNAME*21, STYPE*4, TITLE*80,
     :          UNITS*7, XTYPE, YLAB*80

*     On some systems, such as Sun Sparc, the struct MUST be aligned
*     on a double precision boundary, done here using an equivalence.
*     Failure to do this may result in mysterious "bus errors".
      INCLUDE 'spx.inc'
      INCLUDE 'spc.inc'
      INTEGER   SPC(SPCLEN)
      DOUBLE PRECISION DUMMY
      EQUIVALENCE (SPC,DUMMY)

      COMMON /SPECTRO/ MARS

      DATA TOL /1D-11/
*-----------------------------------------------------------------------
*     Get keyvalues for the required spectral axis type.
      STATUS = SPCXPS (CTYPES, CRVALX, RESTFRQ, RESTWAV, PTYPE, XTYPE,
     :                 RESTREQ, CRVALS, DSDX)
      IF (STATUS.NE.0) THEN
        WRITE (*, 5) STATUS, CTYPES
 5      FORMAT ('ERROR',I2,' from SPCXPS for',A,'.')
        RETURN
      END IF
      CDELTS = CDELTX * DSDX

      STATUS = SPCINI(SPC)

      IF (CTYPES(6:6).EQ.'G') THEN
*       KPNO MARS spectrograph grism parameters.
        STATUS = SPCPUT (SPC, SPC_PV, MARS(0), 0)
        STATUS = SPCPUT (SPC, SPC_PV, MARS(1), 1)
        STATUS = SPCPUT (SPC, SPC_PV, MARS(2), 2)
        STATUS = SPCPUT (SPC, SPC_PV, MARS(3), 3)
        STATUS = SPCPUT (SPC, SPC_PV, MARS(4), 4)
        STATUS = SPCPUT (SPC, SPC_PV, MARS(5), 5)
        STATUS = SPCPUT (SPC, SPC_PV, MARS(6), 6)
      END IF

*     Construct the axis.
      DO 10 J = 1, NAXISJ
        SPEC1(J) = (J - CRPIXJ)*CDELTS
 10   CONTINUE

      WRITE (*, 20) CTYPES, CRVALS+SPEC1(1), CRVALS+SPEC1(NAXISJ),
     :              CDELTS
 20   FORMAT (A,' (CRVALk+w) range: ',1PE13.6,' to ',1PE13.6,', step: ',
     :        1PE13.6)

*     Initialize.
      STATUS = SPCPUT (SPC, SPC_FLAG, 0, 0)
      STATUS = SPCPUT (SPC, SPC_CRVAL, CRVALS, 0)
      STATUS = SPCPUT (SPC, SPC_RESTFRQ, RESTFRQ, 0)
      STATUS = SPCPUT (SPC, SPC_RESTWAV, RESTWAV, 0)
      STATUS = SPCPUT (SPC, SPC_TYPE, CTYPES, 0)
      STATUS = SPCPUT (SPC, SPC_CODE, CTYPES(6:8), 0)

*     Convert the first to the second.
      STATUS = SPCX2S(SPC, NAXISJ, 1, 1, SPEC1, SPEC2, STAT1)
      IF (STATUS.NE.0) THEN
        WRITE (*, 30) STATUS
 30     FORMAT ('SPCX2S ERROR',I2,'.')
        RETURN
      END IF

*     Convert the second back to the first.
      STATUS = SPCS2X(SPC, NAXISJ, 1, 1, SPEC2, CLOS, STAT2)
      IF (STATUS.NE.0) THEN
        WRITE (*, 40) STATUS
 40     FORMAT ('SPCS2X ERROR',I2,'.')
        RETURN
      END IF

*     Test closure.
      NFAIL = 0
      RESIDMAX = 0D0
      STATUS = SPCGET (SPC, SPC_TYPE, STYPE, 0)
      DO 80 J = 1, NAXISJ
        IF (STAT1(J).NE.0) THEN
          WRITE (*, 50) CTYPES, SPEC1(J), STYPE, STAT1(J)
 50       FORMAT (A,': w =',1PE20.12,' -> ',A,' = ???, stat = ',I2)
          GO TO 80
        END IF

        IF (STAT2(J).NE.0) THEN
          WRITE (*, 60) CTYPES, SPEC1(J), STYPE, SPEC2(J), STAT2(J)
 60       FORMAT (A,': w =',1PE20.12,' -> ',A,' =',1PE20.12,
     :            ' -> w = ???, stat = ',I2)
          GO TO 80
        END IF

        RESID = ABS((CLOS(J) - SPEC1(J))/CDELTS)
        IF (RESID.GT.RESIDMAX) RESIDMAX = RESID

        IF (RESID.GT.TOL) THEN
          NFAIL = NFAIL + 1
          WRITE (*, 70) CTYPES, SPEC1(J), STYPE, SPEC2(J), CLOS(J),
     :                  RESID
 70       FORMAT (A,': w =',1PE20.12,' -> ',A,' =',1PE20.12,' ->',/,
     :           '          w =',1PE20.12,',  resid =',1PE8.1)
        END IF
 80   CONTINUE

      WRITE (*, 90) CTYPES, RESIDMAX
 90   FORMAT (A,': Maximum closure residual =',1PE8.1,' pixel')

*     Draw graph.
      CALL PGBBUF()
      CALL PGERAS()

      XMIN = REAL(CRVALS + SPEC1(1))
      XMAX = REAL(CRVALS + SPEC1(NAXISJ))
      YMIN = REAL(SPEC2(1)) - XMIN
      YMAX = YMIN
      DO 100 J = 1, NAXISJ
        X(J) = REAL(J)
        Y(J) = REAL(SPEC2(J) - (CRVALS + SPEC1(J)))
        IF (Y(J).GT.YMAX) YMAX = Y(J)
        IF (Y(J).LT.YMIN) YMIN = Y(J)
 100  CONTINUE

      J = INT(CRPIXJ+1)
      IF (Y(J).LT.0D0) then
        TMP  = YMIN
        YMIN = YMAX
        YMAX = TMP
      END IF

      CALL PGASK(0)
      CALL PGENV(1.0, REAL(NAXISJ), YMIN, YMAX, 0, -1)

      CALL PGSCI(1)
      CALL PGBOX('ABNTS', 0.0, 0, 'BNTS', 0.0, 0)

      STATUS = SPCTYP (CTYPES, STYPE, SCODE, SNAME, UNITS, PTYPE, XTYPE,
     :                 RESTREQ)
      DO 110 J = 21, 1, -1
        IF (SNAME(J:J).NE.' ') GO TO 120
 110  CONTINUE
 120  YLAB  = SNAME(:J) // ' - correction ' // UNITS
      TITLE = CTYPES // ':  CRVALk + w ' // UNITS
      CALL PGLAB('Pixel coordinate', YLAB, TITLE)

      CALL PGAXIS('N', 0.0, YMAX, REAL(NAXISJ), YMAX, XMIN, XMAX, 0.0,
     :            0, -0.5, 0.0, 0.5, -0.5, 0.0)

      CALL PGAXIS('N', REAL(NAXISJ), ymin, REAL(NAXISJ), YMAX,
     :            REAL(YMIN/CDELTS), REAL(YMAX/CDELTS), 0.0, 0, 0.5,
     :            0.0, 0.5, 0.1, 0.0)
      CALL PGMTXT('R', 2.2, 0.5, 0.5, 'Pixel offset')

      CALL PGLINE(NAXISJ, X, Y)
      CALL PGSCI(7)
      CALL PGPT1(REAL(CRPIXJ), 0.0, 24)
      CALL PGEBUF()

      WRITE (*, '(A,$)') 'Type <RETURN> for next page: '
      READ (*, *, END=130)
 130  WRITE (*, *)

      CLOSURE = NFAIL

      END

ext: f

max_stars_repo_name: kapteyn-astro/kapteyn

max_issues_repo_name: kernsuite-debian/kapteyn

max_forks_repo_name: kernsuite-debian/kapteyn


text:
C     LAST UPDATE 16/03/98
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE BGCSYM(B,BB,SMOO,TENS,PC1,PC2)
      IMPLICIT NONE
C
C Purpose: Fits the background by the following method. The pixels are
C          assigned to a radial bin. The values held in this bin are
C          sorted and pixel values in the range PC1 to PC2 are averaged
C          to form a background value for this bin. When the binning is
C          complete, the radial bin values are used to calculate spline
C          coefficients. Background values for the 2D image are then
C          interpolated using these coefficients.
C
C Calls   2: SORT , CURVS
C Called by: LSQINT
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Parameters:
C
      INTEGER MAXBIN,MAXDIM,MAXDAT,MAXNUM
      PARAMETER (MAXBIN=256,MAXDIM=512,MAXDAT=MAXDIM**2,MAXNUM=2048)
C
C Scalar arguments:
C
      REAL SMOO,TENS,PC1,PC2
C
C Array arguments:
C
      REAL B(MAXDAT),BB(MAXDAT)
C
C Common Arrays:
C
      REAL AD(MAXDAT),PSIG(MAXDAT)
C
C Common scalars:
C
      REAL DELR,RMIN,DELZ,ZMIN
      INTEGER IFPIX,ILPIX,IFRAST,ILRAST,NPIX,NRAST,IBACK
      LOGICAL SIGMA,POLAR
C
C Local arrays:
C
      REAL CSYB(MAXBIN),CSYD(MAXBIN),YS(MAXBIN),YSP(MAXBIN),
     &     WRK(9*MAXBIN),PIXBIN(MAXNUM)
      INTEGER INDEX(MAXNUM)
C
C Local Scalars:
C
      INTEGER NBIN,I,J,M,NB,IER,IBAD,NP,IP,IW,IBCK,IW1,IW2,I1,I2
      REAL RS,ZS,DINC,DS,DMAX,RMAX,ZMAX,S,EPS,D1,D2,BCK
C
C Common blocks:
C
      COMMON /INSPEC/ DELR,RMIN,DELZ,ZMIN,POLAR
      COMMON /IMDATA/ NPIX,NRAST,IBACK,AD,PSIG,SIGMA
      COMMON /IMAGE / IFPIX,ILPIX,IFRAST,ILRAST
C
C External function:
C
      REAL CURV2
      EXTERNAL CURV2
C
C-----------------------------------------------------------------------
C
      WRITE(6,1000)
      WRITE(4,1000)
C
C========Set B to zero
C
      DO 20 J=1,NRAST
         DO 10 I=1,NPIX
            M = (J-1)*NPIX + I
            B(M) = 0.0
 10      CONTINUE
 20   CONTINUE
C
C========Initialize limits and increment for binning
C
      IF(POLAR)THEN
         DMAX = MAX(ABS(RMIN),FLOAT(NPIX+1)*DELR+RMIN)
         NBIN = NPIX
      ELSE
         RMAX = MAX(ABS(RMIN),FLOAT(NPIX+1)*DELR+RMIN)
         ZMAX = MAX(ABS(ZMIN),FLOAT(NRAST+1)*DELZ+ZMIN)
         DMAX = SQRT(RMAX*RMAX+ZMAX*ZMAX)
         NBIN = NINT(SQRT(FLOAT(NPIX*NPIX+NRAST*NRAST)))
      ENDIF
      IF(NBIN.GT.MAXBIN)NBIN = MAXBIN
      DINC = DMAX/FLOAT(NBIN)
C
C========Loop over pixels and bin into d*
C
      IBAD = 0
      IBCK = 0
      DO 70 NB=1,NBIN
         D1 = FLOAT(NB-1)*DINC
         D2 = D1 + DINC
         NP = 0
         DO 50 J=1,NRAST
            IF(POLAR)THEN
               I1 = MAX(INT((D1-RMIN)/DELR)-1,1)
               I2 = MIN(INT((D2-RMIN)/DELR)+2,NPIX)
            ELSE
               ZS = (FLOAT(J)-0.5)*DELZ + ZMIN
               I1 = MAX(INT((SQRT(D1*D1-ZS*ZS)-RMIN)/DELR)-1,1)
               I2 = MIN(INT((SQRT(D2*D2-ZS*ZS)-RMIN)/DELR)+2,NPIX)
            ENDIF
            DO 40 I=I1,I2
               IF(POLAR)THEN
                  DS = (FLOAT(I)-0.5)*DELR + RMIN
               ELSE
                  RS = (FLOAT(I)-0.5)*DELR + RMIN
                  DS = SQRT(RS*RS+ZS*ZS)
               ENDIF
               IF(DS.GT.D1.AND.DS.LE.D2)THEN
                  M = (J-1)*NPIX + I
                  IF(AD(M).GT.-0.99E+30)THEN
                     IF(NP.LT.MAXNUM)THEN
                        NP = NP + 1
                        PIXBIN(NP) = AD(M)
                     ENDIF
                  ENDIF
               ENDIF
 40         CONTINUE
 50      CONTINUE
C
C========Sort pixel values and sum over desired interval
C
         IF(NP.GT.0)THEN
            CALL SORT(PIXBIN,NP,INDEX)
            IP = 0
            BCK = 0.0
            IW1 = MIN(INT(PC1*FLOAT(NP))+1,NP)
            IW2 = MIN(INT(PC2*FLOAT(NP))+1,NP)
            DO 60 IW=IW1,IW2
               IP = IP + 1
               BCK = BCK + PIXBIN(INDEX(IW))
 60         CONTINUE
C
C========Form average
C
            IBCK = IBCK + 1
            CSYB(IBCK) = BCK/FLOAT(IP)
            CSYD(IBCK) = (FLOAT(NB)-0.5)*DINC
         ELSE
            IBAD = IBAD + 1
         ENDIF
 70   CONTINUE
      WRITE(6,1005)IBCK,DINC,CSYD(1),CSYD(IBCK)
      WRITE(4,1005)IBCK,DINC,CSYD(1),CSYD(IBCK)
C
C========Calculate spline coefficients for interpolating radial bins
C
      IF(IBCK.GT.1)THEN
         EPS = SQRT(2.0/FLOAT(IBCK))
         S = SMOO*FLOAT(IBCK)
         CALL CURVS(IBCK,CSYD,CSYB,1.0,1,S,EPS,YS,YSP,TENS,WRK,IER)
         IF(IER.NE.0)THEN
            WRITE(6,2000)IER
            WRITE(4,2000)IER
            STOP 'Fatal error'
         ENDIF
C
C========Loop over image pixels to interpolate background values
C
         DO 90 J=IFRAST,ILRAST
            IF(.NOT.POLAR)ZS = (FLOAT(J)-0.5)*DELZ + ZMIN
            DO 80 I=IFPIX,ILPIX
               M = (J-1)*NPIX + I
               IF(POLAR)THEN
                  DS = (FLOAT(I)-0.5)*DELR + RMIN
               ELSE
                  RS = (FLOAT(I)-0.5)*DELR + RMIN
                  DS = SQRT(RS*RS+ZS*ZS)
               ENDIF
               IF(AD(M).GT.-0.99E+30)THEN
                  B(M) = CURV2(DS,IBCK,CSYD,YS,YSP,TENS)
               ENDIF
 80         CONTINUE
 90      CONTINUE
      ENDIF
C
C========Pass A back with the background subtracted
C
      DO 110 J=1,NRAST
         DO 100 I=1,NPIX
            M = (J-1)*NPIX + I
            IF(AD(M).GT.-0.99E+30)THEN
               BB(M) = BB(M) + B(M)
               AD(M) = AD(M) - B(M)
            ENDIF
 100     CONTINUE
 110  CONTINUE
      WRITE(6,1010)
      WRITE(4,1010)
      RETURN
C
 1000 FORMAT(/1X,'BGCSYM: Circularly-symmetric background fitting...'/)
 1005 FORMAT(1X,'Number of radial bins: ',I4,2X,'increment: ',G12.5/
     &       1X,'over d* range: ',G12.5,' to ',G12.5)
 1010 FORMAT(/1X,'Circularly-symmetric background fitting done'/)
 2000 FORMAT(1X,'***Error in spline fitting: FITPACK error ',I3)
      END


text:
!  DLAASRT
!
!  Purpose:
!  ========
!
!  Sort the array D(1:N) in absolute-value descending order and
!  record the permutation in array P(1:N) s.t. new D(j) = old D(P(j))
!
!  Arguments
!  =========
!
!  N        (input) INTEGER
!           The length of the array D.
!
!  D        (input/output) DOUBLE PRECISION array, dimension (N)
!           On entry, the array to be sorted.
!           On exit, D has been sorted into absolute-value descending
!           order.
!
!  P        (output) INTEGER array, dimension (N)
!           P contains the permutation vector.
!
      SUBROUTINE DLAASRT( N, D, P )
!
      INTEGER, INTENT(IN) :: N
      DOUBLE PRECISION, INTENT(INOUT) :: D(N)
      INTEGER, INTENT(OUT) :: P(N)
!
      INTEGER :: I, J, TMP2
      DOUBLE PRECISION :: TMP1
      LOGICAL :: SWAPPED
!
      SWAPPED = .TRUE.
      J = 0
      DO I = 1,N
         P(I) = I
      ENDDO
      DO WHILE (SWAPPED)
        SWAPPED = .FALSE.
        J = J+1
        DO I = 1, N-J
           IF(ABS(D(I))<ABS(D(I+1))) THEN
              TMP1 = D(I)
              TMP2 = P(I)
              D(I) = D(I+1)
              P(I) = P(I+1)
              D(I+1) = TMP1
              P(I+1) = TMP2
              SWAPPED = .TRUE.
           END IF
        END DO
      END DO
      END

max_stars_repo_name: jbringewatt/PolyTimeAlgorithmsforAQC

max_issues_repo_name: jbringewatt/PolyTimeAlgorithmsforAQC

max_forks_repo_name: jbringewatt/PolyTimeAlgorithmsforAQC


text:
c
c =======================================================================
      subroutine outval(val,nvar,mitot,mjtot,mktot,mptr,outgrd,naux,aux)
c =======================================================================
c
      use amr_module
      implicit double precision (a-h,o-z)

      dimension  val(nvar,mitot,mjtot,mktot)
      dimension  aux(naux,mitot,mjtot,mktot)
      logical    outgrd

c ::::::::::::::::::::::OUTVAL :::::::::::::::::::::::::::::::
c print solution and aux. variables to output.
c if cell outside domain, don't print soln. value - nothing
c currently in ghost cells.
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
      if (.not. outgrd) go to 99
      level = node(nestlevel,mptr)
      hx    = hxposs(level)
      hy    = hyposs(level)
      hz    = hzposs(level)
      cornx = rnode(cornxlo,mptr) - nghost*hx
      corny = rnode(cornylo,mptr) - nghost*hy
      cornz = rnode(cornzlo,mptr) - nghost*hz
c
      do 30 k=nghost+1,mktot-nghost
      do 25 j=nghost+1,mjtot-nghost
      do 20 i=nghost+1,mitot-nghost

         x  = cornx + hx*(dble(i)-.5d0)
         y  = corny + hy*(dble(j)-.5d0)
         z  = cornz + hz*(dble(k)-.5d0)
         write(outunit,107) x,y,z,i,j,k,(val(ivar,i,j,k),ivar=1,nvar)
 107     format(2hx=,f7.3,3h y=,f7.3,3h z=,f7.3,
     *          4h, i=,i3,4h, j=,i3 ,4h, k=,i3 ,
     *          ' a= ',5(e11.4,1x))
         if (naux.gt.0) write(outunit,108) (aux(iaux,i,j,k),iaux=1,naux)
 108     format(1x,'aux = ',7(e10.3,1x))

 20   continue
 25   continue
 30   continue

 99   return
      end

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text:
! the network module provides the information about the species we are
! advecting:
!
! nspec      -- the number of species
!
! aion       -- atomic number
! zion       -- proton number
!
! aion_inv   -- 1/aion
!
! spec_names -- the name of the isotope
! short_spec_names -- the abbreviated name of the isotope
!
! This module contains the following routines:
!
!  network_init          -- initialize the isotope properties
!
!  network_species_index -- return the index of the species given its name
!
!  network_finalize      -- do any network cleanup

module network

  use bl_types, only : dp_t
  use actual_network

  implicit none

  private :: dp_t

  logical :: network_initialized = .false.

  ! this will be computed here, not in the actual network
  real(kind=dp_t) :: aion_inv(nspec)

  !$acc declare create(aion_inv)

contains

  subroutine network_init()

    use bl_error_module, only : bl_error
    use bl_constants_module, only : ONE

    implicit none

    ! First, we call the specific network initialization.
    ! This should set the number of species and number of
    ! aux variables, and the components of the species.
    ! Note that the network MUST define nspec and naux
    ! as parameters, or else the compiler will throw an error.

    call actual_network_init()

    ! Check to make sure, and if not, throw an error.

    if ( nspec .le. 0 ) then
       call bl_error("Network cannot have a negative number of species.")
    endif

    if ( naux .lt. 0 ) then
       call bl_error("Network cannot have a negative number of auxiliary variables.")
    endif

    aion_inv(:) = ONE/aion(:)

    !$acc update device(aion_inv)

    network_initialized = .true.

  end subroutine network_init

  function network_species_index(name) result(r)

    character(len=*) :: name
    integer :: r, n

    r = -1

    do n = 1, nspec
       if (name == spec_names(n) .or. name == short_spec_names(n)) then
          r = n
          return
       endif
    enddo

  end function network_species_index

  subroutine network_finalize()
    implicit none

    call actual_network_finalize()

  end subroutine network_finalize

end module network

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
This page contains Cascading Style Sheets (CSS) that control the look of your wiki.  You should be familiar with CSS before changing these files.  To change the look of your wiki, change the following CSS files on this page:

   File(style.css)  The CSS file for controlling colors, fonts, and the like.
   File(layout.css)  The CSS file which controls the margins, paddings and the like.
   File(common.css)  The somewhat general CSS file that you most likely wont need to change.

See the wikispot page on wiki:wikispot:customizing wikis for more information.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


       subroutine pl_write_parms_cal

       use maximum_data_module
       use calibration_data_module
       use hydrograph_module
       use hru_module, only : hru
       use input_file_module
       use plant_module

       implicit none

       character (len=80) :: titldum   !           |title of file
       character (len=80) :: header    !           |header of file
       integer :: eof                  !           |end of file
       logical :: i_exist              !none       |check to determine if file exists
       integer :: mreg                 !none       |end of loop
       integer :: i                    !none       |counter
       integer :: ilum
       integer :: ilum_mx
       integer :: isp                  !none       |counter
       integer :: ielem1               !none       |counter
       integer :: iihru                !none       |counter
       integer :: ihru                 !none       |counter
       integer :: nspu                 !           |
       integer :: ipl                  !           |

       mreg = 0
       eof = 0

       open (107,file="plant_parms.cal",recl = 1500)

        write (107,*) "Calibrated plant parameter file"
        write (107,*) db_mx%plcal_reg
        write (107,*) "header"

      do i = 1, db_mx%plcal_reg

        write (107,*) pl_prms(i)%name, pl_prms(i)%lum_num, pl_prms(i)%parms, "  0"
        write (107,*) "header"

        !! read landscape soft calibration data for each land use and parameter
        if (pl_prms(i)%lum_num > 0) then
          ilum_mx = pl_prms(i)%lum_num * pl_prms(i)%parms
          do ilum = 1, ilum_mx
            select case (pl_prms(i)%prm(ilum)%var)
              case ("pest_stress")
                if (plcal(i)%lum(ilum)%prm%pest_stress > 0.) then
                  pl_prms(i)%prm(ilum)%init_val = plcal(i)%lum(ilum)%prm%pest_stress
                end if
              case ("lai_pot")
                pl_prms(i)%prm(ilum)%init_val = pl_prms(i)%prm(ilum)%init_val + plcal(i)%lum(ilum-pl_prms(i)%lum_num)%prm%lai_pot
                pl_prms(i)%prm(ilum)%init_val = amin1 (pl_prms(i)%prm(ilum)%init_val, pl_prms(i)%prm(ilum)%up)
                pl_prms(i)%prm(ilum)%init_val = Max (pl_prms(i)%prm(ilum)%init_val, pl_prms(i)%prm(ilum)%lo)
              case ("harv_idx")
                pl_prms(i)%prm(ilum)%init_val = pl_prms(i)%prm(ilum)%init_val + plcal(i)%lum(ilum-2*pl_prms(i)%lum_num)%prm%harv_idx
                pl_prms(i)%prm(ilum)%init_val = amin1 (pl_prms(i)%prm(ilum)%init_val, pl_prms(i)%prm(ilum)%up)
                pl_prms(i)%prm(ilum)%init_val = Max (pl_prms(i)%prm(ilum)%init_val, pl_prms(i)%prm(ilum)%lo)
            end select
            !! write parameters to file
            write (107,*) pl_prms(i)%prm(ilum)
          end do
        end if

      end do    !mreg

       close(107)
       return
       end subroutine pl_write_parms_cal


text:
!
! Copyright (c) 2000, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
! test f2008 hypot intrinsic

program p
 use ISO_C_BINDING
 use check_mod

  interface
    subroutine get_expected_f( src1, src2, expct, n ) bind(C)
     use ISO_C_BINDING
      type(C_PTR), value :: src1
      type(C_PTR), value :: src2
      type(C_PTR), value :: expct
      integer(C_INT), value :: n
    end subroutine

    subroutine get_expected_d( src1, src2, expct, n ) bind(C)
     use ISO_C_BINDING
      type(C_PTR), value :: src1
      type(C_PTR), value :: src2
      type(C_PTR), value :: expct
      integer(C_INT), value :: n
    end subroutine
  end interface

  integer, parameter :: N=20
  real*4, target,  dimension(N) :: r_src1
  real*4, target,  dimension(N) :: r_src2
  real*4, target,  dimension(N) :: r_rslt
  real*4, target,  dimension(N) :: r_expct

  real*8, target,  dimension(N) :: d_src1
  real*8, target,  dimension(N) :: d_src2
  real*8, target,  dimension(N) :: d_rslt
  real*8, target,  dimension(N) :: d_expct

  do i =  0,N-1
    r_src1(i+1) = .14*i +  i+1
    r_src2(i+1) = .14*i - i+N
    d_src1(i+1) = .14*i + i+2
    d_src2(i+1) = .14*i - i+N
  enddo

  r_rslt = hypot(r_src1, r_src2)
  d_rslt = hypot(d_src1, d_src2)

  call get_expected_f(C_LOC(r_src1), C_LOC(r_src2), C_LOC(r_expct), N)
  call get_expected_d(C_LOC(d_src1), C_LOC(d_src2), C_LOC(d_expct), N)

  call checkr4( r_rslt, r_expct, N, rtoler=0.0000003)
  call checkr8( d_rslt, d_expct, N, rtoler=0.0000003_8)

!   print *, "r_expct:"
!   print *, r_expct
!   print *, "r_rslt:"
!   print *, r_rslt

!   print *, "d_expct:"
!   print *, d_expct
!   print *, "d_rslt:"
!   print *, d_rslt
end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
!    EmbryoMaker software (General Node Model)
!    Computational model to simulate morphogenetic processes in living organs and tissues.
!    Copyright (C) 2014 Miquel Marin-Riera, Miguel Brun-Usan, Roland Zimm, Tommi Välikangas & Isaac Salazar-Ciudad

!    This program is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    any later version.

!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.

!    You should have received a copy of the GNU General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.

module nexus

  use general
  use genetic
  use growth  !>>>> Miquel 17-6-13
  use death    !>>>> Miquel 18-6-13
  use ecm !>>>> Is 31-8-13
  use single_node
  use mitosis

contains

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

! Declare the interface for POSIX fsync function

subroutine nexe

!interface
!  function fsync (fd) bind(c,name="fsync")
!  use iso_c_binding, only: c_int
!     integer(c_int), value :: fd
!     integer(c_int) :: fsync
!   end function fsync
!end interface

  integer i,j,k,ii,jj,kk,ik,ikk,ick,ret
  real*8 a
  real*8 differ
  real*8 kplast,kvol,dvol
  character*300 cx

! GENETIC REGULATION OF CELL BEHAVIOURS

  ! extracellular matrix secretion

  if (npag(nparam_per_node+4)>0) then ; call should_I_secrete   ; end if

  if(ffu(1)==0)then

    ! cell polarization  THIS ONE SHOULD BE THE FIRST IN HERE
    if (npag(nparam_per_node+8)>0) then ; call polarization ; end if
    !nparam_per_node+9 is to tell that growth is polarized
    ! cell growth
    if (npag(nparam_per_node+1)>0) then ; call should_I_grow        ; end if  !It considers also polar growth in it, with nparam_per_node+9

    !cell division
    !if (npag(nparam_per_node+2)>0) then ;
    call should_I_divide ;
    !end if
    !if (ffu(9)==1) then                                        !>>> Is 5-2-14                                        !>>> Is 5-2-14
    !  do ick=1,ncels                                           !>>> Is 5-2-14
    !    if (cels(ick)%nunodes>cels(ick)%maxsize_for_div) then  !>>> Is 5-2-14 If the cell has too many nodes it divides si o si
    !       call division(ick)                                   !>>> Is 5-2-14
    !      cels(ick)%fase=cels(ick)%fase-1d0                    !>>> Is 5-2-14
    !    end if                                                 !>>> Is 5-2-14
    !  end do                                                   !>>> Is 5-2-14
    !end if                                                     !>>> Is 5-2-14

    ! cell apoptosis
    if (npag(nparam_per_node+3)>0.or.npag(nparam_per_node+14)>0) then ; call should_I_die  ; end if

    ! change the size of the cell required for division
    if (npag(nparam_per_node+10)>0) then; call change_minsize_for_div ; end if

    !nparam_per_node+11 is to orient division according to the chemical polarization and not according to the physical (hertwig) one

    !nparam_per_node+12 the larger the most asymetric (in mass) is the plane of division

    ! change the maximal number of nodes per cell before the cell divides
    if (npag(nparam_per_node+15)>0) then; call change_maxsize_for_div ; end if  !>>> Is 23-3-14

  else
    if (ffu(25)==1) then ; if (npag(nparam_per_node+1)>0) then ; call should_I_grow ; end if  ; end if  !>>> Is 18-4-15 !It considers also polar growth in it, with nparam_per_node+9 ACHTUNG
    if (npag(nparam_per_node+8)>0) then ; call polarization_single ; end if
    if (npag(nparam_per_node+2)>0) then ; call should_I_divide_single ; end if
    ! cell apoptosis
    if (npag(nparam_per_node+3)>0.or.npag(nparam_per_node+14)>0) then ; call should_I_die  ; end if  !>>Miquel20-3-14
    if (npag(nparam_per_node+13)>0) then; call emt ; end if   ! >>> Is 4-10-14
   !call emt_single
  end if

  ! epithelial-mesenchymal transition
  if (npag(nparam_per_node+13)>0) then; call emt ; end if

  if (ffu(5)==1) then ; agex(:nd,1)=abs(node(:nd)%x) ;end if; !external source trick >>> Is 29-6-14
  if (ffu(5)==2) then ; agex(:nd,1)=node(:nd)%z**4    ;end if!((1.1d0+maxval(node(:nd)%z)+node(:nd)%z)**4) !external source trick miguel >>> Is 29-6-14

  if(ffu(5)==3)then !horizontal gradient (x-y plane) centered on the origin of coordinates
    do i=1,nd
      agex(i,6)=1/(1+(node(i)%x**2+node(i)%y**2)) !; print*,"i",i,"gex",agex(i,5)
    end do
  end if

  if (ffu(15)==1)then; !especial conditions, gene 1 is always expressing in the borders !>>Miquel12-5-14
    do i=1,nd
      if(node(i)%hold==1)then
        agex(i,1)=1d0
        !if(node(i)%tipus<3)then
        !  agex(i,1)=1d0
        !else
        !  agex(i,4)=1d0
        !end if
      end if
    end do
  end if

  if (ffu(8)==1.and.nd>2) then   ! IS 23-4-13 this eliminates the nodes that get alone for too long
    ik=1
    do while(ik<=nd) !;print*,"node(",ik,")%talone=",node(ik)%talone
      if (node(ik)%talone>ttalone) then !;print*,ik,"entra mort",node(ik)%tipus
        if (node(ik)%tipus>2) then  !we only delete an epithelial node if its altre is also lonely
          ikk=node(ik)%icel !;print*,"entra mesenq"
          call apoptosis(ik)
          if (ikk>0) then   !notice that if ikk is negative it means that is not a celular node
            if(cels(ikk)%nunodes==0) call eliminate_cell(ikk)  !we eliminate the cell if it has no nodes left
          end if
        else
          if (node(node(ik)%altre)%talone>ttalone) then
            ikk=node(ik)%icel
            call apoptosis(ik)
            if (ikk>0) then   !notice that if ikk is negative it means that is not a celular node
              if(cels(ikk)%nunodes==0) call eliminate_cell(ikk)  !we eliminate the cell if it has no nodes left
            end if
          end if
        end if
        ik=ik+1 ! Is it? >>> Is 16-1-14
      else
        ik=ik+1   ! I know, it is kind of funky but it should be this way, a loop with nd wont do because nd decreases because of apoptosis
      end if
    end do
!    call cellbreak ! to see if the cell is split in two
  end if

! GENETIC REGULATION OF NODE PROPERTIES
  do i=1,nd         ! we update that parameter in each cell that expresses the gene
                    ! WE ONLY UPDATES DE NODES IN WHICH THE GENE IS EXPRESSED, OTHERWISE WE LEAVE IT IS AS IT WAS
    if(node(i)%hold==1) cycle  !we don't want the border cells to perform behaviours because that would alter and possibly break the border  !>>>>Miquel9-1-14
    differ=1-node(i)%diffe !;print*,"differ",differ

    ! DIFFERENTIATION
    if (npag(25)>0) then
      ii=25;a=0.0;
      do k=1,npag(ii)
        kk=whonpag(ii,k) ;
        if (gex(i,kk)>0.0d0) then ;
          a=a+gex(i,kk)*gen(kk)%wa(ii)
        endif
      end do
      node(i)%diffe=node(i)%diffe+a*delta
      if (node(i)%diffe>1.0) node(i)%diffe=1.0
      if (node(i)%diffe<0.0d0) node(i)%diffe=0.0

      nodeo(i)%diffe=node(i)%diffe !>>Miquel17-9-14
    end if

    if(node(i)%tipus<2)then  !this for epithelial nodes  !>>Miquel12-5-14
      j=node(i)%altre
      if(ffu(11)==1) then        !plastic deformation   !>>Miquel5-2-14
        !kplast=node(i)%kplast
        !if(fmeanl(i)<epsilod)then;ki=1/node(i)%rep;else;ki=1/node(i)%you;endif
        !if(fmeanl(j)<epsilod)then;kj=1/node(j)%rep;else;kj=1/node(j)%you;endif
        node(i)%reqp=node(i)%reqp+(node(i)%kplast*fmeanl(i))*delta
        node(j)%reqp=node(j)%reqp+(node(j)%kplast*fmeanl(j))*delta !; print*,"fmeanl",fmeanl(i),fmeanl(j)

        nodeo(i)%reqp=node(i)%reqp  !>>Miquel17-9-14
        nodeo(j)%reqp=node(j)%reqp  !>>Miquel17-9-14

      else
        node(i)%reqp=0
        node(j)%reqp=0
      end if

      if (npag(21)>0) then  ! Contraction by genes
        ii=21 ; a=0 ; b=0
        do k=1,npag(ii) ;
          kk=whonpag(ii,k) !;print*,"ij",i,j,"gex",gex(i,kk),gex(j,kk),"kk",kk
          if (gex(i,kk)>0.0d0.or.gex(j,kk)>0.0d0) then ;
            a=a+gex(i,kk)*gen(kk)%wa(ii)
            b=b+gex(j,kk)*gen(kk)%wa(ii)
          endif
        enddo
        node(i)%reqc=nodeo(i)%reqc+a*differ  !;print*,"a",a,"b",b
        node(j)%reqc=nodeo(j)%reqc+b*differ
        !;print*,"reqc",node(i)%reqc,node(j)%reqc
      else
        node(i)%reqc=0 ;node(j)%reqc=0
      end if

      if(ffu(17)==1) then        !volume conservation   !>>Miquel6-5-14 ! >>> Is 24-5-14
        !kvol=node(i)%kvol
        dvol=0.5*(node(i)%reqcr+node(j)%reqcr-(node(i)%req+node(j)%req))
        node(i)%reqv=node(i)%reqv+node(i)%kvol*dvol*delta
        node(j)%reqv=node(j)%reqv+node(j)%kvol*dvol*delta
        nodeo(i)%reqv=node(i)%reqv  !>>Miquel17-9-14
        nodeo(j)%reqv=node(j)%reqv  !>>Miquel17-9-14

      else
        node(i)%reqv=0 ; node(j)%reqv=0
      end if

      if (ffu(18)==1) call diffusion_of_reqcr       !diffusion of reqcr ! >>> Is 25-5-14

      a=node(i)%da-node(i)%req
      node(i)%req=node(i)%reqcr+node(i)%reqc+node(i)%reqp+node(i)%reqv  !now req is the sum of the req components: growth/apoptosis and contraction/deformation
      if(node(i)%req>df_reqmax) node(i)%req=df_reqmax !put an upper an lower boundary on how much  !>>Miquel28-7-14
      if(node(i)%req<reqmin) node(i)%req=reqmin !the req can be deformed
      node(i)%da=node(i)%req+a

      b=node(j)%da-node(j)%req
      node(j)%req=node(j)%reqcr+node(j)%reqc+node(j)%reqp+node(j)%reqv  !now req is the sum of the req components: growth/apoptosis and contraction/deformation
      if(node(j)%req>df_reqmax) node(j)%req=df_reqmax !put an upper an lower boundary on how much  !>>Miquel28-7-14
      if(node(j)%req<reqmin) node(j)%req=reqmin !the req can be deformed
      node(j)%da=node(j)%req+b
      if (npag(6)>0) then;ii=6;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in req-space units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%da=nodeo(i)%da+a*differ;node(j)%da=nodeo(j)%da+b*differ;
      if (node(i)%da<0) node(i)%da=0.0;if (node(j)%da<0) node(j)%da=0.0;end if
      if (npag(7)>0) then;ii=7;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in req-space units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%you=nodeo(i)%you+a*differ;node(j)%you=nodeo(j)%you+b*differ;
      if (node(i)%you<0) node(i)%you=0.0;if (node(j)%you<0) node(j)%you=0.0;end if
      if (npag(8)>0) then;ii=8;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%adh=nodeo(i)%adh+a*differ;node(j)%adh=nodeo(j)%adh+b*differ;
      if (node(i)%adh<0) node(i)%adh=0.0;if (node(j)%adh<0) node(j)%adh=0.0;end if
      if (npag(9)>0) then;ii=9;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%rep=nodeo(i)%rep+a*differ;node(j)%rep=nodeo(j)%rep+b*differ;
      if (node(i)%rep<0) node(i)%rep=0.0;if (node(j)%rep<0) node(j)%rep=0.0;end if
      if (npag(10)>0) then;ii=10;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%repcel=nodeo(i)%repcel+a*differ;node(j)%repcel=nodeo(j)%repcel+b*differ;
      if (node(i)%repcel<0) node(i)%repcel=0.0;if (node(j)%repcel<0) node(j)%repcel=0.0;end if
      if (npag(11)>0) then;ii=11;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%tor=nodeo(i)%tor+a*differ;node(j)%tor=nodeo(j)%tor+b*differ;
      if (node(i)%tor<0) node(i)%tor=0.0;if (node(j)%tor<0) node(j)%tor=0.0;end if
      if (npag(12)>0) then;ii=12;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%stor=nodeo(i)%stor+a*differ;node(j)%stor=nodeo(j)%stor+b*differ;
      if (node(i)%stor<0) node(i)%stor=0.0;if (node(j)%stor<0) node(j)%stor=0.0;end if
      if (npag(13)>0) then;ii=13;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%reqs=nodeo(i)%reqs+a*differ;node(j)%reqs=nodeo(j)%reqs+b*differ;
      if (node(i)%reqs<0) node(i)%reqs=0.0;if (node(j)%reqs<0) node(j)%reqs=0.0;end if
      if (npag(14)>0) then;ii=14;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%ke=nodeo(i)%ke+a*differ;node(j)%ke=nodeo(j)%ke+b*differ;
      if (node(i)%ke<0) node(i)%ke=0.0;if (node(j)%ke<0) node(j)%ke=0.0;end if
      if (npag(15)>0) then;ii=15;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in force units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%mo=nodeo(i)%mo+a*differ;node(j)%mo=nodeo(j)%mo+b*differ;
      if (node(i)%mo<0) node(i)%mo=0.0;if (node(j)%mo<0) node(j)%mo=0.0;end if
      if (npag(16)>0) then;ii=16;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in req-space units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%dmo=nodeo(i)%dmo+a*differ;node(j)%dmo=nodeo(j)%dmo+b*differ;
      if (node(i)%dmo<0) node(i)%dmo=0.0;if (node(j)%dmo<0) node(j)%dmo=0.0;end if
      if (npag(27)>0) then; ii=27;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in req-space units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%kplast=nodeo(i)%kplast+a*differ;node(j)%kplast=nodeo(j)%kplast+b*differ;
      if (node(i)%kplast<0) node(i)%kplast=0.0;if (node(j)%kplast<0) node(j)%kplast=0.0;end if
      if (npag(28)>0) then;ii=28;a=0 ; b=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; !wa in req-space units
      if (gex(i,kk)>0.0d0) then ; a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;
      if (gex(j,kk)>0.0d0) then ; b=b+gex(j,kk)*gen(kk)%wa(ii) ;endif;enddo;
      node(i)%kvol=nodeo(i)%kvol+a*differ;node(j)%kvol=nodeo(j)%kvol+b*differ;
      if (node(i)%kvol<0) node(i)%kvol=0.0;if (node(j)%kvol<0) node(j)%kvol=0.0;end if

    else if (node(i)%tipus>2)then !this for mesenchyme and ECM

      if (npag(21)>0) then
        ii=21 ; a=0
        do k=1,npag(ii) ;
          kk=whonpag(ii,k)
          if (gex(i,kk)>0.0d0) then ;
            a=a+gex(i,kk)*gen(kk)%wa(ii)
          endif
        enddo
        node(i)%reqc=nodeo(i)%reqc+a*differ*delta !wa in req-space units
      else
        node(i)%reqc=0
      end if

      a=node(i)%da-node(i)%req
      node(i)%req=node(i)%reqcr+node(i)%reqc  !now req is the sum of the req components: growth/apoptosis and contraction/deformation
      if(node(i)%req>df_reqmax) node(i)%req=df_reqmax !put an upper an lower boundary on how much  !>>Miquel28-7-14
      if(node(i)%req<reqmin) node(i)%req=reqmin !the req can be deformed
      node(i)%da=node(i)%req+a

      if (npag(6)>0) then;ii=6;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%da=nodeo(i)%da+a*differ;if (node(i)%da<0) node(i)%da=0.0;end if
      if (npag(7)>0) then;ii=7;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%you=nodeo(i)%you+a*differ;if (node(i)%you<0) node(i)%you=0.0;end if
      if (npag(8)>0) then;ii=8;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%adh=nodeo(i)%adh+a*differ;if (node(i)%adh<0) node(i)%adh=0.0;end if
      if (npag(9)>0) then;ii=9;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%rep=nodeo(i)%rep+a*differ;if (node(i)%rep<0) node(i)%rep=0.0;end if
      if (npag(10)>0) then;ii=10;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%repcel=nodeo(i)%repcel+a;if (node(i)%repcel<0) node(i)%repcel=0.0;end if
      if (npag(11)>0) then;ii=11;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%tor=nodeo(i)%tor+a*differ;if (node(i)%tor<0) node(i)%tor=0.0;end if
      if (npag(12)>0) then;ii=12;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%stor=nodeo(i)%stor+a*differ;if (node(i)%stor<0) node(i)%stor=0.0;end if
      if (npag(13)>0) then;ii=13;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%reqs=nodeo(i)%reqs+a*differ;if (node(i)%reqs<0) node(i)%reqs=0.0;end if
      if (npag(14)>0) then;ii=14;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%ke=nodeo(i)%ke+a*differ;if (node(i)%ke<0) node(i)%ke=0.0;end if
      if (npag(15)>0) then;ii=15;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii) ;endif;enddo;node(i)%mo=nodeo(i)%mo+a*differ;if (node(i)%mo<0) node(i)%mo=0.0;end if
      if (npag(16)>0) then;ii=16;a=0 ; do k=1,npag(ii) ; kk=whonpag(ii,k) ; if (gex(i,kk)>0.0d0) then ;
      a=a+gex(i,kk)*gen(kk)%wa(ii);endif;enddo;node(i)%dmo=nodeo(i)%dmo+a*differ;if (node(i)%dmo<0) node(i)%dmo=0.0;end if

    end if
  end do

  ! if a node gets 3 times the original size of node 1 we fucking kill the program
  if (mod(getot,100).eq.0) then !>>> Is 4-2-14
    if (ffu(7)==1) then  ! if a node gets twice its original size we fucking kill the program ! >>> Is 4-2-14
      do i=1,nd                                                                               ! >>> Is 4-2-14
        if (node(i)%da>ramax) then                        ! >>> Is 4-2-14
          node(i)%da=ramax
        end if                                                                                ! >>> Is 4-2-14
        !if (node(i)%req>ramax) then   !maximal req now is controlled above in the same sub. !>>Miquel7-8-14
        !  node(i)%req=ramax
        !end if
      end do                                                                                ! >>> Is 4-2-14
    end if                                                                                ! >>> Is 4-2-14
  end if                                                                                ! >>> Is 4-2-14

! We check in here if all cells are already differentiated, if they are we stop the simulations

!>>> Is 4-4-14
if (ffu(10)==1) then
  do i=1,nd
    if (node(i)%tipus==4) cycle !pfh, ignore ecm
    if (node(i)%diffe<1.0d0) return
  end do
  print *,""
  print *," THIS IS THE END all cells are differentiated and then the simulations stop",trim(carg)//trim(noff),trim(nofi)
  print *,""
  call writesnap                                                                            !>>> Is 25-2-14
  if (len_trim(carg)/=0) then
    print *,trim(carg)//trim(noff),"ki"
    print *,trim(carg)//trim(nofi),"kii"
    open(23,file=trim(carg)//"t",iostat=i)
    print *,"making...",trim(carg)//"t"
    write(23,*,ERR=46) trim(carg)//trim(nofi)
    flush(23)
    !ret = fsync(fnum(23))

    ! Handle possible error
    !if (ret /= 0) stop "Error calling FSYNC"
    open(23,file=trim(carg)//"t",iostat=ii)
    print *,trim(carg)//"t",ii,"iostat"
    read(23,*,END=45,ERR=48) cx
print *,""
print *,cx,"cx HERE"
print *,""
    !close(23)
    call flush(23)
    print *,"done with iostat",ii
    cx="ls -alrt "//trim(carg)//"t"
    call system(cx)
  end if
  stop
45 print *,"end of file error"
  stop
46 print *,"error in writing",trim(carg)//"t"
  stop
48 print *,"other end"
  stop
end if
!>>> Is 4-4-14

end subroutine

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

subroutine cellbreak   !>>> 17-1-14
  real*8  ax,ay,az,ida
  integer doapo(nd)
  integer nocon(nd,nd)
  integer i,j,k,ii,jj,kk,iii,jjj,kkk,ik

  ! now we check that the cell is not split in two or more parts
  doapo=0
  nocon=0

  do i=1,ncels

    kkk=0
    do j=1,cels(i)%nunodes
      ii=cels(i)%node(j)
      if (node(ii)%marge==0) then ; kkk=ii ; ida=node(kkk)%da ; exit ; end if
    end do

    if (kkk==0) then !it means that the cell has no nucleus and then IT MUST DIE!!!!   !>>> Is 11-6-14
      do j=1,cels(i)%nunodes  !>>> Is 11-6-14
        ii=cels(i)%node(j)    !>>> Is 29-6-14
!print *,ii,i,j,cels(i)%node(:cels(i)%nunodes),ii,"ii",nd
        doapo(ii)=1           !>>> Is 11-6-14
      end do                  !>>> Is 11-6-14
      kkk=1  ! >>> Is 11-6-14
    else
      do j=1,cels(i)%nunodes
        ii=cels(i)%node(j)
        ax=node(ii)%x ; ay=node(ii)%y ; az=node(ii)%z
        if (sqrt((ax-node(kkk)%x)**2+(ay-node(kkk)%y)**2+(az-node(ii)%z)**2)<node(ii)%da+ida) then
          doapo(ii)=1
        end if
      end do
    end if   ! >>> Is 11-6-14

    doapo(kkk)=1
    do j=1,cels(i)%nunodes
      ii=cels(i)%node(j)
      ax=node(ii)%x ; ay=node(ii)%y ; az=node(ii)%z
      ida=node(ii)%da
      do jj=1,cels(i)%nunodes
        if (j==jj) cycle
        iii=cels(i)%node(jj)
        if (sqrt((ax-node(iii)%x)**2+(ay-node(iii)%y)**2+(az-node(iii)%z)**2)<node(iii)%da+ida) then
          nocon(ii,iii)=1
          nocon(iii,ii)=1
          if (doapo(ii)==1) then
            doapo(iii)=1
          else
            if (doapo(iii)==1) doapo(ii)=1
          end if
        end if
      end do
    end do

    do k=1,cels(i)%nunodes/2
      do j=1,cels(i)%nunodes
        ii=cels(i)%node(j)
        do jj=1,cels(i)%nunodes
          if (j==jj) cycle
          iii=cels(i)%node(jj)
          if (nocon(ii,iii)==1) then
            if (doapo(ii)==1) then
              doapo(iii)=1
            else
              if (doapo(iii)==1) doapo(ii)=1
            end if
          end if
        end do
      end do
    end do
  end do

  ik=1
  do while(ik<=nd)
    if (doapo(ik)==0) then
      if (node(ik)%tipus>2) then  !we only delete an epithelial node if its altre is also lonely
        call apoptosis(ik)
      else
        if (doapo(node(ik)%altre)==0) then
          call apoptosis(ik)
        end if
      end if
      ik=ik+1
    else
      ik=ik+1   ! I know, it is kind of funky but it should be this way, a loop with nd wont do because nd decreases because of apoptosis
    end if
  end do

end subroutine

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

subroutine polarization
integer:: celd,nnod,tipi,ggr,ccen
real*8::a,b,c,d,e,ax,ay,az,bx,by,bz,cx,cy,cz,ix,iy,iz,alfa,s
      do celd=1,ncels
        tipi=cels(celd)%ctipus
        nnod=cels(celd)%nunodes
        if (nnod==0) cycle      ! >>> Is 10-5-14
        iy=1d10 ; cx=0d0 ; cy=0d0 ; cz=0d0
	a=cels(celd)%cex ; b=cels(celd)%cey ; c=cels(celd)%cez

        do i=1,nnod                                                     ! (gen) in the centroid (in the closest node)
          j=cels(celd)%node(i)
          if(node(j)%tipus==1.or.node(j)%tipus==3)then !in epithelial cells, polarity is planar, so we only take one layer of nodes >>>Miquel22-10-13
            d=sqrt((node(j)%x-a)**2+(node(j)%y-b)**2+(node(j)%z-c)**2)
          end if
          if(d.le.iy)then;iy=d;ccen=j;endif
        end do

        alfa=0.0d0                                                 ! concentration in the central node
        do k=1,npag(nparam_per_node+8)
          kk=whonpag(nparam_per_node+8,k)
          if (gex(ccen,kk)>0.0d0) then
            alfa=alfa+gex(ccen,kk)*gen(kk)%wa(nparam_per_node+8)   ! wa in units of probability such that it makes things to go from 0 to 1
          end if
        end do

        ix=0d0 ; iy=0d0 ; iz=0d0                                        ! vector of the gradient within a cell
        do i=1,nnod
            j=cels(celd)%node(i)
            if(node(j)%tipus==1.or.node(j)%tipus==3)then
              d=sqrt((node(j)%x-a)**2+(node(j)%y-b)**2+(node(j)%z-c)**2)
              if (d<epsilod) cycle
              d=1d0/d                                                   ! module of radial vectors to get unitary vectors
              s=0.0d0
              do k=1,npag(nparam_per_node+8)
                kk=whonpag(nparam_per_node+8,k)
                if (gex(j,kk)>0.0d0) then
                  s=s+gex(j,kk)*gen(kk)%wa(nparam_per_node+8)
                end if
              end do
              ix=ix+((node(j)%x-a)*d)*(s-alfa)                   ! and ignore shape/size effects
              iy=iy+((node(j)%y-b)*d)*(s-alfa)
              iz=iz+((node(j)%z-c)*d)*(s-alfa)
            end if
        end do

        if((ix.eq.0).and.(iy.eq.0).and.(iz.eq.0))then            ! if the gene has uniform expresion, the vector is random ! >>>Miguel1-7-14
          call random_number(a)                                  ! >>>Miguel1-7-14
          k=int(a*nvaloq)+1                                      ! >>>Miguel1-7-14
          cels(celd)%polx=particions_esfera(k,1)                 ! >>>Miguel1-7-14
          cels(celd)%poly=particions_esfera(k,2)                 ! >>>Miguel1-7-14
          cels(celd)%polz=particions_esfera(k,3)                 ! >>>Miguel1-7-14
        else                                                     ! >>>Miguel1-7-14
          a=ix**2+iy**2+iz**2
          if(a==0)then
            cels(celd)%polx=0d0 ; cels(celd)%poly=0d0 ; cels(celd)%polz=0d0	! unitary resultant vector (gradient polarization)
          else
            d=1d0/sqrt(a)
            cels(celd)%polx=ix*d ; cels(celd)%poly=iy*d ; cels(celd)%polz=iz*d	! unitary resultant vector (gradient polarization)
          end if
          if((ix.eq.0d0).and.(iy.eq.0d0).and.(iz.eq.0d0))then                     ! miguel27-11-13
            cels(celd)%polx=0d0 ; cels(celd)%poly=0d0 ; cels(celd)%polz=0d0
          endif   ! miguel27-11-13
        endif                                                    ! >>>Miguel1-7-14
      end do
end subroutine

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

subroutine polarizationisaac
integer i,j,k,ii,kk
real*8 a,b,c,s,sx,sy,sz,d,aa,bb,cc

do i=1,ncels
  a=cels(i)%cex ; b=cels(i)%cey ; c=cels(i)%cez
  sx=0.0d0      ; sy=0.0d0      ; sz=0.0d0
  if (node(cels(i)%node(1))%tipus<3) then
    do j=1,cels(i)%nunodes
      ii=cels(i)%node(j)
      if (node(ii)%tipus==1) then
        s=0.0d0
        do k=1,npag(nparam_per_node+8)
          kk=whonpag(nparam_per_node+8,k)
          if (gex(ii,kk)>0.0d0) then
            s=s+gex(ii,kk)*gen(kk)%wa(nparam_per_node+8)*delta
          end if
        end do
        if (s/=0.0d0) then
          aa=node(ii)%x-a ; bb=node(ii)%y-b ; cc=node(ii)%z-c
          d=s/sqrt(aa**2+bb**2+cc**2)
          sx=sx+d*aa ; sy=sy+d*bb ; sz=sz+d*cc
        end if
      end if
    end do
    aa=sqrt(sx**2+sy**2+sz**2)
    if (aa>0.0d0) then
      d=1d0/aa
      cels(i)%polx=sx*d ; cels(i)%poly=sy*d ; cels(i)%polz=sz*d
    else
      cels(i)%polx=0.0d0 ; cels(i)%poly=0.0d0 ; cels(i)%polz=0.0d0
    end if
  else
    if (node(cels(i)%node(1))%tipus==3) then
      do j=1,cels(i)%nunodes
        ii=cels(i)%node(j)
        s=0.0d0
        do k=1,npag(nparam_per_node+8)
          kk=whonpag(nparam_per_node+8,k)
          if (gex(ii,kk)>0.0d0) then
            s=s+gex(ii,kk)*gen(kk)%wa(nparam_per_node+8)
          end if
        end do
        if (s/=0.0d0) then
          aa=node(ii)%x-a ; bb=node(ii)%y-b ; cc=node(ii)%z-c
          sx=sx+s*aa ; sy=sy+s*bb ; sz=sz+s*cc
        end if
      end do
      aa=sqrt(sx**2+sy**2+sz**2)
      if (aa>0.0d0) then  !epsilod) then
        d=1d0/aa
        cels(i)%polx=sx*d ; cels(i)%poly=sy*d ; cels(i)%polz=sz*d
      else
        cels(i)%polx=0.0d0 ; cels(i)%poly=0.0d0 ; cels(i)%polz=0.0d0
      end if
    end if
  end if
end do
end subroutine

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

subroutine change_minsize_for_div  ! updates the size required for dividint according to gene expression nparam_per_node+10
integer ick,j,k,ii,kk
real*8 a,b,c,s,sx,sy,sz,d

do ick=1,ncels
  s=0.0d0
  do j=1,cels(ick)%nunodes
    ii=cels(ick)%node(j)
    do k=1,npag(nparam_per_node+10)
      kk=whonpag(nparam_per_node+10,k)
      if (gex(ii,kk)>0.0d0) then
        s=s+gex(ii,kk)*gen(kk)%wa(nparam_per_node+10)  !wa in units of number of nodes but it can be roughly understood as space-req units
      end if                                           ! THIS WA CAN BE NEGATIVE
    end do
  end do
  s=s*delta
!print *,delta,"ACHTUNG"
!  if (s>0.0d0) then
    s=s/cels(ick)%nunodes !this way %minsize is independent of cell size !>>>>Miquel2-12-13
    cels(ick)%minsize_for_div=cels(ick)%minsize_for_div+s  !this can make a SUDDEN CHANGE
    if (cels(ick)%minsize_for_div<1) cels(ick)%minsize_for_div=1
!  end if
end do

end subroutine

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

!>>> Is 5-2-14
subroutine change_maxsize_for_div  ! updates the size required for dividint according to gene expression nparam_per_node+10
integer ick,j,k,ii,kk
real*8 a,b,c,s,sx,sy,sz,d

do ick=1,ncels
  s=0.0d0
  do j=1,cels(ick)%nunodes
    ii=cels(ick)%node(j)
    do k=1,npag(nparam_per_node+15)
      kk=whonpag(nparam_per_node+15,k)
      if (gex(ii,kk)>0.0d0) then
        s=s+gex(ii,kk)*gen(kk)%wa(nparam_per_node+15)  ! THIS WA MAY BE NEGATIVE
      end if
    end do
  end do
  s=s*delta
!  if (s>0.0d0) then
    s=s/cels(ick)%nunodes !this way %minsize is independent of cell size !>>>>Miquel2-12-13
    cels(ick)%maxsize_for_div=cels(ick)%maxsize_for_div+s  !wa in units of space-req roughly this can make a SUDDEN CHANGE
    if (cels(ick)%maxsize_for_div<1) cels(ick)%maxsize_for_div=1
!  end if
end do

!>>> Is 5-2-14

end subroutine

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

subroutine emt !epithelial-mensenchymal transitions

               ! this simply makes the tipus equal to 3 but it needs to approach both sides of the originally epithelial cell so that
               ! the two parts do not drift away since they are more far away than da (they are at reqs normally)
               ! this transition has to be sudden otherwise other existing forces may impede the nodes from getting close enough to each other
               ! that wouldnt be realistic since it would be as if the cell would explode.

integer ick,j,k,ii,kk,iii,jjj,kkk,iv
real*8 a,b,c,d,aa,bb,cc,dd,s

e: do ick=1,ncels
  if (cels(ick)%ctipus==3) cycle
  s=0.0
  do j=1,cels(ick)%nunodes
    ii=cels(ick)%node(j)
    do k=1,npag(nparam_per_node+13)
      kk=whonpag(nparam_per_node+13,k)
      a=gex(ii,kk) !don't change this, weird stuff will happen if you do !>>Miquel18-12-14
      s=s+a*gen(kk)%wa(nparam_per_node+13)
   !if(gex(ii,kk)>epsilod) print*,"emt gene",kk,"node",ii,"cell",ick,"conc",gex(ii,kk)
    end do
  end do
  cels(ick)%temt=cels(ick)%temt+s*delta/real(cels(ick)%nunodes)

  if (cels(ick)%temt>=1.0d0) then       ! wa in units of probability: this is an arbitrary value but it is fine since the rest can be re-scaled
    !so this cell goes to emt
 !print*,"EMT HAPPENS CELL",ick
    cels(ick)%ctipus=3  ! >>> Is 4-1-14
    do jjj=1,cels(ick)%nunodes
      iii=cels(ick)%node(jjj)
      if (node(iii)%tipus==1) then
        iv=node(iii)%altre
        a=node(iv)%x-node(iii)%x ; b=node(iv)%y-node(iii)%y ; c=node(iv)%z-node(iii)%z
        d=1d0/sqrt(a**2+b**2+c**2)
        a=a*d ; b=b*d ; c=c*d
        aa=0.5d0*(node(iv)%x+node(iii)%x) ; bb=0.5d0*(node(iv)%y+node(iii)%y) ; cc=0.5d0*(node(iv)%z+node(iii)%z)
        dd=(node(iii)%da+node(iv)%da)*0.1d0
        node(iii)%x=aa-a*dd
        node(iii)%y=bb-b*dd
        node(iii)%z=cc-c*dd
        node(iv)%x=aa+a*dd
        node(iv)%y=bb+b*dd
        node(iv)%z=cc+c*dd
        node(iii)%altre=0
        node(iv)%altre=0
        node(iii)%tipus=3
        node(iv)%tipus=3
        !if only one node per cell both nodes of the epithelial cell get a nucleus
        if (ffu(1)==1) then ; node(iii)%marge=0 ; node(iv)%marge=0 ; end if ! >>> Is 10-10-14
      end if
    end do
  end if
end do e

end subroutine

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

subroutine diffusion_of_reqcr
  real*8 hreqcr(nd),hreqp(nd),hreqc(nd)
  integer i,j,k,ii,jj,kk
  real*8 a,b,c,d

  do i=1,nd
    a=0.0d0 ; b=0.0 ; c=0.0d0
    do ii=1,nneigh(i)
      k=neigh(i,ii)
      if (node(i)%icel/=node(k)%icel) cycle    ! only within the same cell
      if (node(i)%tipus/=node(k)%tipus) cycle  ! only within the same side of the cell
      if (node(i)%tipus>2) cycle               ! only for epithelial cells
      d=dneigh(i,ii)
      a=a+(node(k)%reqcr-node(i)%reqcr) !/(d+1d0)
      !b=b+(node(k)%reqc-node(i)%reqc) !/(d+1d0)
      !c=c+(node(k)%reqp-node(i)%reqp) !/(d+1d0)
    end do
    hreqcr(i)=a*dif_req ! >>> 11-6-14
    !hreqc(i)=b*dif_req  ! >>> 11-6-14
    !hreqp(i)=c*dif_req  ! >>> 11-6-14
  end do
  do i=1,nd
    node(i)%reqcr=node(i)%reqcr+delta*hreqcr(i)
    !node(i)%reqc=node(i)%reqc+delta*hreqc(i)
    !node(i)%reqp=node(i)%reqp+delta*hreqp(i)
  end do
end subroutine

end module nexus

max_stars_repo_name: isalazarciudad/EmbryoMaker

max_issues_repo_name: sebjameswml/EmbryoMaker

max_forks_repo_name: sebjameswml/EmbryoMaker


text:
!  compile using "gfortran -xf95-cpp-input"
!    to invoke the preprocessor

PROGRAM main

use emombm
IMPLICIT NONE
	INTEGER i, ice
	REAL*8, DIMENSION(J) :: Quold, Qlold, Tsold, Qunew, Qlnew, Tsnew
	REAL*8 :: Taveold, Tavenew
	REAL*8, DIMENSION(J+1) :: Uu,Ul,tau

	call emombm_setup()

	print *, ' '
	print *, 'Running the EMomBM once'
	print *, ' '

	if (.NOT. readInputFile ) then
		!  fixed initial condition
		Quold = dble(2.0)*f0/sqrt(2.)*sin(phiq)
		Qlold = Quold
		Taveold = 11.
		Tsold = 40.-40.*2./pi*phiq
	else
		!  get initial condition from a file
		print *, 'Loading initial conditions from file ', inputfilename
		open(unit=10,file=inputfilename)
		read(unit=10,fmt=*) Quold,Qlold
		read(unit=10,fmt=*) Taveold
		read(unit=10,fmt=*) Tsold
		close(unit=10)
	end if

	!  main time loop
time:	DO i=1,numsteps
			call timestep_emombm(Quold,Qlold,Taveold,Tsold,Qunew,Qlnew,Tavenew,Tsnew,ice)
			Quold = Qunew
			Qlold = Qlnew
			Taveold = Tavenew
			Tsold = Tsnew
			!PRINT *, ice
		END DO time

	! final winds
	CALL solveU(Quold,Qlold,Uu,Ul,J)
	tau = eps*(3*Ul - Uu)/2

	100 format(F14.8)
	101 format(E20.10)
	open(unit=10,file='output.dat')
	write(unit=10,fmt=101) Qunew,Qlnew
	write(unit=10,fmt=100) Tavenew
	write(unit=10,fmt=100) Tsnew
	write(unit=10,fmt=101) Uu,Ul,tau
	close(unit=10)
	print *, 'Results saved to file output.dat'
END PROGRAM main

max_stars_repo_name: brian-rose/EMomBM

max_issues_repo_name: brian-rose/EMomBM

max_forks_repo_name: brian-rose/EMomBM


text:
	SUBROUTINE GLTOU  ( np, xl, yl, xm, ym, iret )
C************************************************************************
C* GLTOU								*
C* 									*
C* This subroutine converts points in rotated linear intermediate       *
C* coordinates into non-rotated linear intermediate coordinates.  	*
C* 									*
C* GLTOU  ( NP, XL, YL, XM, YM, IRET )					*
C*									*
C* Input parameters:							*
C*	NP		INTEGER		Number of points		*
C* 	XL (NP)		REAL		X coordinates 			*
C* 	YL (NP)		REAL		Y coordinates 			*
C* 									*
C* Output parameters: 							*
C* 	XM (NP)		REAL		X coordinates			*
C* 	YM (NP)		REAL		Y coordinates			*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Gilbert/NCEP	 8/06	Modified from GLTOW                     *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'ERROR.PRM'
	INCLUDE		'XYDEF.CMN'
C*
	REAL		xl (*), yl (*), xm (*), ym (*)
C--------------------------------------------------------------------------
C*	Select code for current mode.
C
        iret = NORMAL
C
C*	    Map mode.
C
	IF  ( igmode .eq. 1 )  THEN
	    IF  ( .not. mset )  THEN
		iret = NIPROJ
C*
	    ELSE IF  ( mclass .eq. MCCYL )  THEN
C
C*		Unrotate L coordinates.
C
		IF ( rotaml ) CALL GROTTL ( 'M', -1, np, xl, yl, ier )
C*
	    END IF
C
	ELSE IF ( igmode .ne. 2 ) THEN
	    iret = NIMODE
	END IF

        DO I = 1, np
            xm(i) = xl(i)
            ym(i) = yl(i)
        ENDDO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
        SUBROUTINE BRAINFORT(PROG,N,INF,OUF,F)	!Stand strong!
Converts the Brain*uck in PROG into the equivalent furrytran source...
         CHARACTER*(*) PROG	!The Brain*uck source.
         INTEGER N		!A size for the STORE.
         INTEGER INF,OUF,F	!I/O unit numbers.
         INTEGER L		!A stepper.
         INTEGER LABEL,NLABEL,INDEEP,STACK(66)	!Labels cause difficulty.
         CHARACTER*1 C		!The operation of the moment.
         CHARACTER*36 SOURCE	!A scratchpad.
          WRITE (F,1) PROG,N	!The programme heading.
    1     FORMAT (6X,"PROGRAM BRAINFORT",/,	!Name it.
     1     "Code: ",A,/				!Show the provenance.
     2     6X,"CHARACTER*1 STORE(",I0,")",/	!Declare the working memory.
     3     6X,"INTEGER D",/			!The finger to the cell of the moment.
     4     6X,"STORE = CHAR(0)",/		!Clear to nulls, not spaces.
     5     6X,"D = 1",/)			!Start the data finger at the first cell.
          NLABEL = 0		!No labels seen.
          INDEEP = 0		!So, the stack is empty.
          LABEL = 0		!And the current label is absent.
          L = 1			!Start at the start.
Chug through the PROG.
          DO WHILE(L.LE.LEN(PROG))	!And step through to the end.
            C = PROG(L:L)		!The code of the moment.
            SELECT CASE(C)		!What to do?
             CASE(">")			!Move the data finger forwards one.
              WRITE (SOURCE,2) "D = D + ",RATTLE(">")	!But, catch multiple steps.
             CASE("<")			!Move the data finger back one.
              WRITE (SOURCE,2) "D = D - ",RATTLE("<")	!Rather than a sequence of one steps.
             CASE("+")			!Increment the fingered datum by one.
              WRITE (SOURCE,2) "STORE(D) = CHAR(ICHAR(STORE(D)) + ",	!Catching multiple increments.
     1         RATTLE("+"),")"						!And being careful over the placement of brackets.
             CASE("-")			!Decrement the fingered datum by one.
              WRITE (SOURCE,2) "STORE(D) = CHAR(ICHAR(STORE(D)) - ",	!Catching multiple decrements.
     1         RATTLE("-"),")"						!And closing brackets.
             CASE(".")			!Write a character.
              WRITE (SOURCE,2) "WRITE (",OUF,",'(A1,$)') STORE(D)"	!Using the given output unit.
             CASE(",")			!Read a charactger.
              WRITE (SOURCE,2) "READ (",INF,",'(A1)') STORE(D)"		!And the input unit.
             CASE("[")			!A label!
              NLABEL = NLABEL + 1		!Labels come in pairs due to [...]
              LABEL = 2*NLABEL - 1		!So this belongs to the [.
              INDEEP = INDEEP + 1		!I need to remember when later the ] is encountered.
              STACK(INDEEP) = LABEL + 1		!This will be the other label.
              WRITE (SOURCE,2) "IF (ICHAR(STORE(D)).EQ.0) GO TO ",	!So, go thee, therefore.
     1         STACK(INDEEP)			!Its placement will come, all going well.
             CASE("]")			!The end of a [...] pair.
              LABEL = STACK(INDEEP)		!This was the value of the label to be, now to be placed.
              WRITE (SOURCE,2) "IF (ICHAR(STORE(D)).NE.0) GO TO ",	!The conditional part
     1         LABEL - 1			!The branch back destination is known by construction.
              INDEEP = INDEEP - 1		!And we're out of the [...] sequence's consequences.
             CASE DEFAULT		!All others are ignored.
              SOURCE = "CONTINUE"		!So, just carry on.
            END SELECT			!Enough of all that.
    2       FORMAT (A,I0,A)	!Text, an integer, text.
Cast forth the statement.
            IF (LABEL.LE.0) THEN	!Is a label waiting?
              WRITE (F,3) SOURCE		!No. Just roll the source.
    3         FORMAT (<6 + 2*MIN(12,INDEEP)>X,A)!With indentation.
             ELSE			!But if there is a label,
              WRITE (F,4) LABEL,SOURCE		!Slightly more complicated.
    4         FORMAT (I5,<1 + 2*MIN(12,INDEEP)>X,A)	!I align my labels rightwards...
              LABEL = 0				!It is used.
            END IF			!So much for that statement.
            L = L + 1		!Advance to the next command.
          END DO		!And perhaps we're finished.

Closedown.
          WRITE (F,100)		!No more source.
  100     FORMAT (6X,"END")	!So, this is the end.
         CONTAINS	!A function with odd effects.
          INTEGER FUNCTION RATTLE(C)	!Advances thrugh multiple C, counting them.
           CHARACTER*1 C	!The symbol.
            RATTLE = 1		!We have one to start with.
    1       IF (L.LT.LEN(PROG)) THEN	!Further text to look at?
              IF (PROG(L + 1:L + 1).EQ.C) THEN	!Yes. The same again?
              	L = L + 1		!Yes. Advance the finger to it.
                RATTLE = RATTLE + 1	!Count another.
                GO TO 1			!And try again.
              END IF			!Rather than just one at a time.
            END IF			!Curse the double evaluation of WHILE(L < LEN(PROG) & ...)
          END FUNCTION RATTLE	!Computers excel at counting.
        END SUBROUTINE BRAINFORT!Their only need be direction as to what to count...
      END MODULE BRAIN	!Simple in itself.

      PROGRAM POKE	!A tester.
      USE BRAIN		!In a rather bad way.
      CHARACTER*1 STORE(30000)	!Probably rather more than is needed.
      CHARACTER*(*) HELLOWORLD	!Believe it or not...
      PARAMETER (HELLOWORLD = "++++++++[>++++[>++>+++>+++>+<<<<-]"
     1 //" >+>+>->>+[<]<-]>>.>---.+++++++..+++.>>.<-.<.+++.------"
     2 //".--------.>>+.>++.")
      INTEGER F
      KBD = 5		!Standard input.
      MSG = 6		!Standard output.
      F = 10

      STORE = CHAR(0)	!Scrub.

c      CALL RUN(HELLOWORLD,STORE)	!Have a go.

      OPEN (F,FILE="BrainFort.for",STATUS="REPLACE",ACTION="WRITE")
      CALL BRAINFORT(HELLOWORLD,30000,KBD,MSG,F)
      END	!Enough.

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData