text:
!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_sf_sfclay

 REAL    , PARAMETER ::  VCONVC=1.
 REAL    , PARAMETER ::  CZO=0.0185
 REAL    , PARAMETER ::  OZO=1.59E-5

 REAL,   DIMENSION(0:1000 ),SAVE          :: PSIMTB,PSIHTB

CONTAINS

!-------------------------------------------------------------------
   SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w,                    &
                     CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM,      &
                     ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
                     XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
                     U10,V10,TH2,T2,Q2,                            &
                     GZ1OZ0,WSPD,BR,ISFFLX,DX,                     &
                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
                     KARMAN,EOMEG,STBOLT,                          &
                     P1000mb,                                      &
                     ids,ide, jds,jde, kds,kde,                    &
                     ims,ime, jms,jme, kms,kme,                    &
                     its,ite, jts,jte, kts,kte,                    &
                     ustm,ck,cka,cd,cda,isftcflx,iz0tlnd           )
!-------------------------------------------------------------------
      IMPLICIT NONE
!-------------------------------------------------------------------
!-- U3D         3D u-velocity interpolated to theta points (m/s)
!-- V3D         3D v-velocity interpolated to theta points (m/s)
!-- T3D         temperature (K)
!-- QV3D        3D water vapor mixing ratio (Kg/Kg)
!-- P3D         3D pressure (Pa)
!-- dz8w        dz between full levels (m)
!-- CP          heat capacity at constant pressure for dry air (J/kg/K)
!-- G           acceleration due to gravity (m/s^2)
!-- ROVCP       R/CP
!-- R           gas constant for dry air (J/kg/K)
!-- XLV         latent heat of vaporization for water (J/kg)
!-- PSFC        surface pressure (Pa)
!-- ZNT         roughness length (m)
!-- UST         u* in similarity theory (m/s)
!-- USTM        u* in similarity theory (m/s) without vconv correction
!               used to couple with TKE scheme
!-- PBLH        PBL height from previous time (m)
!-- MAVAIL      surface moisture availability (between 0 and 1)
!-- ZOL         z/L height over Monin-Obukhov length
!-- MOL         T* (similarity theory) (K)
!-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
!-- PSIM        similarity stability function for momentum
!-- PSIH        similarity stability function for heat
!-- XLAND       land mask (1 for land, 2 for water)
!-- HFX         upward heat flux at the surface (W/m^2)
!-- QFX         upward moisture flux at the surface (kg/m^2/s)
!-- LH          net upward latent heat flux at surface (W/m^2)
!-- TSK         surface temperature (K)
!-- FLHC        exchange coefficient for heat (W/m^2/K)
!-- FLQC        exchange coefficient for moisture (kg/m^2/s)
!-- CHS         heat/moisture exchange coefficient for LSM (m/s)
!-- QGH         lowest-level saturated mixing ratio
!-- QSFC        ground saturated mixing ratio
!-- U10         diagnostic 10m u wind
!-- V10         diagnostic 10m v wind
!-- TH2         diagnostic 2m theta (K)
!-- T2          diagnostic 2m temperature (K)
!-- Q2          diagnostic 2m mixing ratio (kg/kg)
!-- GZ1OZ0      log(z/z0) where z0 is roughness length
!-- WSPD        wind speed at lowest model level (m/s)
!-- BR          bulk Richardson number in surface layer
!-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
!-- DX          horizontal grid size (m)
!-- SVP1        constant for saturation vapor pressure (kPa)
!-- SVP2        constant for saturation vapor pressure (dimensionless)
!-- SVP3        constant for saturation vapor pressure (K)
!-- SVPT0       constant for saturation vapor pressure (K)
!-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
!-- EP2         constant for specific humidity calculation
!               (R_d/R_v) (dimensionless)
!-- KARMAN      Von Karman constant
!-- EOMEG       angular velocity of earth's rotation (rad/s)
!-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
!-- ck          enthalpy exchange coeff at 10 meters
!-- cd          momentum exchange coeff at 10 meters
!-- cka         enthalpy exchange coeff at the lowest model level
!-- cda         momentum exchange coeff at the lowest model level
!-- isftcflx    =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd, =2 Garratt
!-- iz0tlnd     =0 Carlson-Boland, =1 Czil_new
!-- ids         start index for i in domain
!-- ide         end index for i in domain
!-- jds         start index for j in domain
!-- jde         end index for j in domain
!-- kds         start index for k in domain
!-- kde         end index for k in domain
!-- ims         start index for i in memory
!-- ime         end index for i in memory
!-- jms         start index for j in memory
!-- jme         end index for j in memory
!-- kms         start index for k in memory
!-- kme         end index for k in memory
!-- its         start index for i in tile
!-- ite         end index for i in tile
!-- jts         start index for j in tile
!-- jte         end index for j in tile
!-- kts         start index for k in tile
!-- kte         end index for k in tile
!-------------------------------------------------------------------
      INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte
!
      INTEGER,  INTENT(IN )   ::        ISFFLX
      REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
      REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
      REAL,     INTENT(IN )   ::        P1000mb
!
      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(IN   )   ::                           dz8w

      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(IN   )   ::                           QV3D, &
                                                              P3D, &
                                                              T3D

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(IN   )               ::             MAVAIL, &
                                                             PBLH, &
                                                            XLAND, &
                                                              TSK
      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(OUT  )               ::                U10, &
                                                              V10, &
                                                              TH2, &
                                                               T2, &
                                                               Q2, &
                                                             QSFC

!
      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)               ::             REGIME, &
                                                              HFX, &
                                                              QFX, &
                                                               LH, &
                                                          MOL,RMOL
!m the following 5 are change to memory size
!
      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
                                                        PSIM,PSIH

      REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
                INTENT(IN   )   ::                            U3D, &
                                                              V3D

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(IN   )               ::               PSFC

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)   ::                            ZNT, &
                                                              ZOL, &
                                                              UST, &
                                                              CPM, &
                                                             CHS2, &
                                                             CQS2, &
                                                              CHS

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)   ::                      FLHC,FLQC

      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
                INTENT(INOUT)   ::                                 &
                                                              QGH

      REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX

      REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
                INTENT(OUT)     ::              ck,cka,cd,cda,ustm

      INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND

! LOCAL VARS

      REAL,     DIMENSION( its:ite ) ::                       U1D, &
                                                              V1D, &
                                                             QV1D, &
                                                              P1D, &
                                                              T1D

      REAL,     DIMENSION( its:ite ) ::                    dz8w1d

      INTEGER ::  I,J

      DO J=jts,jte
        DO i=its,ite
          dz8w1d(I) = dz8w(i,1,j)
        ENDDO

        DO i=its,ite
           U1D(i) =U3D(i,1,j)
           V1D(i) =V3D(i,1,j)
           QV1D(i)=QV3D(i,1,j)
           P1D(i) =P3D(i,1,j)
           T1D(i) =T3D(i,1,j)
        ENDDO

        !  Sending array starting locations of optional variables may cause
        !  troubles, so we explicitly change the call.

        CALL SFCLAY1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,               &
                CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
                CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j),   &
                ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j),    &
                MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j),  &
                XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j),     &
                U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),        &
                Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j),      &
                QSFC(ims,j),LH(ims,j),                             &
                GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX,     &
                SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,  &
                P1000mb,                                           &
                ids,ide, jds,jde, kds,kde,                         &
                ims,ime, jms,jme, kms,kme,                         &
                its,ite, jts,jte, kts,kte                          &
#if ( EM_CORE == 1 )
                ,isftcflx,iz0tlnd,                                 &
                USTM(ims,j),CK(ims,j),CKA(ims,j),                  &
                CD(ims,j),CDA(ims,j)                               &
#endif
                                                                   )
      ENDDO

   END SUBROUTINE SFCLAY

!-------------------------------------------------------------------
   SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
                     CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
                     ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,      &
                     XLAND,HFX,QFX,TSK,                            &
                     U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH,              &
                     QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX,             &
                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
                     KARMAN,EOMEG,STBOLT,                          &
                     P1000mb,                                      &
                     ids,ide, jds,jde, kds,kde,                    &
                     ims,ime, jms,jme, kms,kme,                    &
                     its,ite, jts,jte, kts,kte,                    &
                     isftcflx, iz0tlnd,                            &
                     ustm,ck,cka,cd,cda                            )
!-------------------------------------------------------------------
      IMPLICIT NONE
!-------------------------------------------------------------------
      REAL,     PARAMETER     ::        XKA=2.4E-5
      REAL,     PARAMETER     ::        PRT=1.

      INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        J
!
      INTEGER,  INTENT(IN )   ::        ISFFLX
      REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
      REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
      REAL,     INTENT(IN )   ::        P1000mb

!
      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(IN   )               ::             MAVAIL, &
                                                             PBLH, &
                                                            XLAND, &
                                                              TSK
!
      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(IN   )               ::             PSFCPA

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)               ::             REGIME, &
                                                              HFX, &
                                                              QFX, &
                                                         MOL,RMOL
!m the following 5 are changed to memory size---
!
      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
                                                        PSIM,PSIH

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)   ::                            ZNT, &
                                                              ZOL, &
                                                              UST, &
                                                              CPM, &
                                                             CHS2, &
                                                             CQS2, &
                                                              CHS

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)   ::                      FLHC,FLQC

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(INOUT)   ::                                 &
                                                              QGH

      REAL,     DIMENSION( ims:ime )                             , &
                INTENT(OUT)     ::                        U10,V10, &
                                                TH2,T2,Q2,QSFC,LH

      REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX

! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
      REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::  dz8w1d

      REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      UX, &
                                                               VX, &
                                                             QV1D, &
                                                              P1D, &
                                                              T1D

      REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
                INTENT(OUT)     ::              ck,cka,cd,cda,ustm

      INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND

! LOCAL VARS

      REAL,     DIMENSION( its:ite )        ::                 ZA, &
                                                        THVX,ZQKL, &
                                                           ZQKLP1, &
                                                           THX,QX, &
                                                            PSIH2, &
                                                            PSIM2, &
                                                           PSIH10, &
                                                           PSIM10, &
                                                           DENOMQ, &
                                                          DENOMQ2, &
                                                          DENOMT2, &
                                                            WSPDI, &
                                                           GZ2OZ0, &
                                                           GZ10OZ0
!
      REAL,     DIMENSION( its:ite )        ::                     &
                                                      RHOX,GOVRTH, &
                                                            TGDSA
!
      REAL,     DIMENSION( its:ite)         ::          SCR3,SCR4
      REAL,     DIMENSION( its:ite )        ::         THGB, PSFC
!
      INTEGER                               ::                 KL

      INTEGER ::  N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10

      REAL    ::  PL,THCON,TVCON,E1
      REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
      REAL    ::  DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
      REAL    ::  FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,RESTAR2
!-------------------------------------------------------------------
      KL=kte

      DO i=its,ite
! PSFC cb
         PSFC(I)=PSFCPA(I)/1000.
      ENDDO
!
!----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
!
      DO 5 I=its,ite
        TGDSA(I)=TSK(I)
! PSFC cb
!        THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP
        THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
    5 CONTINUE
!
!-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
!     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
!
!     *** NOTE ***
!         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
!         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
!         TENDENCIES.
!
   10 CONTINUE

!     DO 24 I=its,ite
!        UX(I)=U1D(I)
!        VX(I)=V1D(I)
!  24 CONTINUE

   26 CONTINUE

!.....SCR3(I,K) STORE TEMPERATURE,
!     SCR4(I,K) STORE VIRTUAL TEMPERATURE.

      DO 30 I=its,ite
! PL cb
         PL=P1D(I)/1000.
         SCR3(I)=T1D(I)
!         THCON=(100./PL)**ROVCP
         THCON=(P1000mb*0.001/PL)**ROVCP
         THX(I)=SCR3(I)*THCON
         SCR4(I)=SCR3(I)
         THVX(I)=THX(I)
         QX(I)=0.
   30 CONTINUE
!
      DO I=its,ite
         QGH(I)=0.
         FLHC(I)=0.
         FLQC(I)=0.
         CPM(I)=CP
      ENDDO
!
!     IF(IDRY.EQ.1)GOTO 80
      DO 50 I=its,ite
         QX(I)=QV1D(I)
         TVCON=(1.+EP1*QX(I))
         THVX(I)=THX(I)*TVCON
         SCR4(I)=SCR3(I)*TVCON
   50 CONTINUE
!
      DO 60 I=its,ite
        E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
!  for land points QSFC can come from previous time step
        if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1)
! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
! Q2SAT = QGH IN LSM
        E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
        PL=P1D(I)/1000.
        QGH(I)=EP2*E1/(PL-E1)
        CPM(I)=CP*(1.+0.8*QX(I))
   60 CONTINUE
   80 CONTINUE

!-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
!     LEVEL, AND THE LAYER THICKNESSES.

      DO 90 I=its,ite
        ZQKLP1(I)=0.
        RHOX(I)=PSFC(I)*1000./(R*SCR4(I))
   90 CONTINUE
!
      DO 110 I=its,ite
           ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
  110 CONTINUE
!
      DO 120 I=its,ite
         ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
  120 CONTINUE
!
      DO 160 I=its,ite
        GOVRTH(I)=G/THX(I)
  160 CONTINUE

!-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
!     AKB(1976), EQ(12).

      DO 260 I=its,ite
        GZ1OZ0(I)=ALOG(ZA(I)/ZNT(I))
        GZ2OZ0(I)=ALOG(2./ZNT(I))
        GZ10OZ0(I)=ALOG(10./ZNT(I))
        IF((XLAND(I)-1.5).GE.0)THEN
          ZL=ZNT(I)
        ELSE
          ZL=0.01
        ENDIF
        WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))

        TSKV=THGB(I)*(1.+EP1*QSFC(I))
        DTHVDZ=(THVX(I)-TSKV)
!  Convective velocity scale Vc and subgrid-scale velocity Vsg
!  following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
!                                ... HONG Aug. 2001
!
!       VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
!      Use Beljaars over land, old MM5 (Wyngaard) formula over water
        if (xland(i).lt.1.5) then
        fluxc = max(hfx(i)/rhox(i)/cp                    &
              + ep1*tskv*qfx(i)/rhox(i),0.)
        VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
        else
        IF(-DTHVDZ.GE.0)THEN
          DTHVM=-DTHVDZ
        ELSE
          DTHVM=0.
        ENDIF
        VCONV = 2.*SQRT(DTHVM)
        endif
! Mahrt and Sun low-res correction
        VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
        WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
        WSPD(I)=AMAX1(WSPD(I),0.1)
        BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
!  IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
        IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
!jdf
        RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
!jdf

  260 CONTINUE

!
!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
!
!
!     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
!     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
!
!     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
!
!        1. BR .GE. 0.2;
!               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
!
!        2. BR .LT. 0.2 .AND. BR .GT. 0.0;
!               REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
!               (REGIME=2),
!
!        3. BR .EQ. 0.0
!               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
!
!        4. BR .LT. 0.0
!               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
!
!CCCCC

      DO 320 I=its,ite
!CCCCC
!CC     REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
!CC          IF(BR(I).LT.0..AND.HOL(I,J).GT.1.5)GOTO 310
        IF(BR(I).LT.0.)GOTO 310
!
!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
!
        IF(BR(I).LT.0.2)GOTO 270
        REGIME(I)=1.
        PSIM(I)=-10.*GZ1OZ0(I)
!    LOWER LIMIT ON PSI IN STABLE CONDITIONS
        PSIM(I)=AMAX1(PSIM(I),-10.)
        PSIH(I)=PSIM(I)
        PSIM10(I)=10./ZA(I)*PSIM(I)
        PSIM10(I)=AMAX1(PSIM10(I),-10.)
        PSIH10(I)=PSIM10(I)
        PSIM2(I)=2./ZA(I)*PSIM(I)
        PSIM2(I)=AMAX1(PSIM2(I),-10.)
        PSIH2(I)=PSIM2(I)

!       1.0 over Monin-Obukhov length
        IF(UST(I).LT.0.01)THEN
           RMOL(I)=BR(I)*GZ1OZ0(I) !ZA/L
        ELSE
           RMOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I)) !ZA/L
        ENDIF
        RMOL(I)=AMIN1(RMOL(I),9.999) ! ZA/L
        RMOL(I) = RMOL(I)/ZA(I) !1.0/L

        GOTO 320
!
!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
!
  270   IF(BR(I).EQ.0.0)GOTO 280
        REGIME(I)=2.
        PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))
!    LOWER LIMIT ON PSI IN STABLE CONDITIONS
        PSIM(I)=AMAX1(PSIM(I),-10.)
!.....AKB(1976), EQ(16).
        PSIH(I)=PSIM(I)
        PSIM10(I)=10./ZA(I)*PSIM(I)
        PSIM10(I)=AMAX1(PSIM10(I),-10.)
        PSIH10(I)=PSIM10(I)
        PSIM2(I)=2./ZA(I)*PSIM(I)
        PSIM2(I)=AMAX1(PSIM2(I),-10.)
        PSIH2(I)=PSIM2(I)

        ! Linear form: PSIM = -0.5*ZA/L; e.g, see eqn 16 of
        ! Blackadar, Modeling the nocturnal boundary layer, Preprints,
        ! Third Symposium on Atmospheric Turbulence Diffusion and Air Quality,
        ! Raleigh, NC, 1976
        ZOL(I) = BR(I)*GZ1OZ0(I)/(1.00001-5.0*BR(I))

        if ( ZOL(I) .GT. 0.5 ) then ! linear form ok
           ! Holtslag and de Bruin, J. App. Meteor 27, 689-704, 1988;
           ! see also, Launiainen, Boundary-Layer Meteor 76,165-179, 1995
           ! Eqn (8) of Launiainen, 1995
           ZOL(I) = ( 1.89*GZ1OZ0(I) + 44.2 ) * BR(I)*BR(I)    &
                + ( 1.18*GZ1OZ0(I) - 1.37 ) * BR(I)
           ZOL(I)=AMIN1(ZOL(I),9.999)
        end if

        ! 1.0 over Monin-Obukhov length
        RMOL(I)= ZOL(I)/ZA(I)

        GOTO 320
!
!-----CLASS 3; FORCED CONVECTION:
!
  280   REGIME(I)=3.
        PSIM(I)=0.0
        PSIH(I)=PSIM(I)
        PSIM10(I)=0.
        PSIH10(I)=PSIM10(I)
        PSIM2(I)=0.
        PSIH2(I)=PSIM2(I)

        IF(UST(I).LT.0.01)THEN
          ZOL(I)=BR(I)*GZ1OZ0(I)
        ELSE
          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
        ENDIF

        RMOL(I) = ZOL(I)/ZA(I)

        GOTO 320
!
!-----CLASS 4; FREE CONVECTION:
!
  310   CONTINUE
        REGIME(I)=4.
        IF(UST(I).LT.0.01)THEN
          ZOL(I)=BR(I)*GZ1OZ0(I)
        ELSE
          ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
        ENDIF
        ZOL10=10./ZA(I)*ZOL(I)
        ZOL2=2./ZA(I)*ZOL(I)
        ZOL(I)=AMIN1(ZOL(I),0.)
        ZOL(I)=AMAX1(ZOL(I),-9.9999)
        ZOL10=AMIN1(ZOL10,0.)
        ZOL10=AMAX1(ZOL10,-9.9999)
        ZOL2=AMIN1(ZOL2,0.)
        ZOL2=AMAX1(ZOL2,-9.9999)
        NZOL=INT(-ZOL(I)*100.)
        RZOL=-ZOL(I)*100.-NZOL
        NZOL10=INT(-ZOL10*100.)
        RZOL10=-ZOL10*100.-NZOL10
        NZOL2=INT(-ZOL2*100.)
        RZOL2=-ZOL2*100.-NZOL2
        PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
        PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
        PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))
        PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
        PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))
        PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))

!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
!---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
!       PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
!       PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
        PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
        PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
        PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
        PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
! AHW: mods to compute ck, cd
        PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))

        RMOL(I) = ZOL(I)/ZA(I)

  320 CONTINUE
!
!-----COMPUTE THE FRICTIONAL VELOCITY:
!     ZA(1982) EQS(2.60),(2.61).
!
      DO 330 I=its,ite
        DTG=THX(I)-THGB(I)
        PSIX=GZ1OZ0(I)-PSIM(I)
        PSIX10=GZ10OZ0(I)-PSIM10(I)
!     LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
!     ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
        PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)

        IF((XLAND(I)-1.5).GE.0)THEN
          ZL=ZNT(I)
        ELSE
          ZL=0.01
        ENDIF
        PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)
        PSIT2=GZ2OZ0(I)-PSIH2(I)
        PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-PSIH2(I)
! AHW: mods to compute ck, cd
        PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-PSIH10(I)
        IF ( PRESENT(ISFTCFLX) ) THEN
           IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
! v3.1
!             Z0Q = 1.e-4 + 1.e-3*(MAX(0.,UST(I)-1.))**2
! hfip1
!             Z0Q = 0.62*2.0E-5/UST(I) + 1.E-3*(MAX(0.,UST(I)-1.5))**2
! v3.2
              Z0Q = 1.e-4
              PSIQ=ALOG(ZA(I)/Z0Q)-PSIH(I)
              PSIT=PSIQ
              PSIQ2=ALOG(2./Z0Q)-PSIH2(I)
              PSIQ10=ALOG(10./Z0Q)-PSIH10(I)
              PSIT2=PSIQ2
           ENDIF
           IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN
! AHW: Garratt formula: Calculate roughness Reynolds number
!        Kinematic viscosity of air (linear approc to
!                 temp dependence at sea levle)
              VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
!!            VISC=1.5E-5
              RESTAR=UST(I)*ZNT(I)/VISC
              RESTAR2=2.48*SQRT(SQRT(RESTAR))-2.
              PSIT=GZ1OZ0(I)-PSIH(I)+RESTAR2
              PSIQ=GZ1OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
              PSIT2=GZ2OZ0(I)-PSIH2(I)+RESTAR2
              PSIQ2=GZ2OZ0(I)-PSIH2(I)+2.28*SQRT(SQRT(RESTAR))-2.
              PSIQ10=GZ10OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
           ENDIF
        ENDIF
        IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
           Ck(I)=(karman/psix10)*(karman/psiq10)
           Cd(I)=(karman/psix10)*(karman/psix10)
           Cka(I)=(karman/psix)*(karman/psiq)
           Cda(I)=(karman/psix)*(karman/psix)
        ENDIF
        IF ( PRESENT(IZ0TLND) ) THEN
           IF ( IZ0TLND.EQ.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN
              ZL=ZNT(I)
!             CZIL RELATED CHANGES FOR LAND
              VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
              RESTAR=UST(I)*ZL/VISC
!             Modify CZIL according to Chen & Zhang, 2009

              CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) )

              PSIT=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
              PSIQ=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
              PSIT2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
              PSIQ2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)

           ENDIF
        ENDIF
! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
        UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
! TKE coupling: compute ust without vconv for use in tke scheme
        WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
        IF ( PRESENT(USTM) ) THEN
        USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
        ENDIF
        U10(I)=UX(I)*PSIX10/PSIX
        V10(I)=VX(I)*PSIX10/PSIX
        TH2(I)=THGB(I)+DTG*PSIT2/PSIT
        Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
!        T2(I) = TH2(I)*(PSFC(I)/100.)**ROVCP
        T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP
!       LATER Q2 WILL BE OVERWRITTEN FOR LAND POINTS IN SURFCE
!       QA2(I,J) = Q2(I)
!       UA10(I,J) = U10(I)
!       VA10(I,J) = V10(I)
!       write(*,1002)UST(I),KARMAN*WSPD(I),PSIX,KARMAN*WSPD(I)/PSIX
!
        IF((XLAND(I)-1.5).LT.0.)THEN
          UST(I)=AMAX1(UST(I),0.1)
        ENDIF
        MOL(I)=KARMAN*DTG/PSIT/PRT
        DENOMQ(I)=PSIQ
        DENOMQ2(I)=PSIQ2
        DENOMT2(I)=PSIT2
  330 CONTINUE
!
  335 CONTINUE

!-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:

      DO i=its,ite
        QFX(i)=0.
        HFX(i)=0.
      ENDDO

      IF (ISFFLX.EQ.0) GOTO 410

!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).

      DO 360 I=its,ite
        IF((XLAND(I)-1.5).GE.0)THEN
          ZNT(I)=CZO*UST(I)*UST(I)/G+OZO
! AHW: change roughness length, and hence the drag coefficients Ck and Cd
          IF ( PRESENT(ISFTCFLX) ) THEN
             IF ( ISFTCFLX.NE.0 ) THEN
!               ZNT(I)=10.*exp(-9.*UST(I)**(-.3333))
                ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333))
                ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
                ZNT(I)=MIN(ZNT(I),2.85e-3)
                ZNT(I)=MAX(ZNT(I),1.27e-7)
             ENDIF
          ENDIF
          ZL = ZNT(I)
        ELSE
          ZL = 0.01
        ENDIF
        FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I)
!       FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
!               ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
        DTTHX=ABS(THX(I)-THGB(I))
        IF(DTTHX.GT.1.E-5)THEN
          FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))
!         write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
 1001   format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
        ELSE
          FLHC(I)=0.
        ENDIF
  360 CONTINUE

!
!-----COMPUTE SURFACE MOIST FLUX:
!
!     IF(IDRY.EQ.1)GOTO 390
!
      DO 370 I=its,ite
        QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
        QFX(I)=AMAX1(QFX(I),0.)
        LH(I)=XLV*QFX(I)
  370 CONTINUE

!-----COMPUTE SURFACE HEAT FLUX:
!
  390 CONTINUE
      DO 400 I=its,ite
        IF(XLAND(I)-1.5.GT.0.)THEN
          HFX(I)=FLHC(I)*(THGB(I)-THX(I))
          IF ( PRESENT(ISFTCFLX) ) THEN
             IF ( ISFTCFLX.NE.0 ) THEN
! AHW: add dissipative heating term
                HFX(I)=HFX(I)+RHOX(I)*USTM(I)*USTM(I)*WSPDI(I)
             ENDIF
          ENDIF
        ELSEIF(XLAND(I)-1.5.LT.0.)THEN
          HFX(I)=FLHC(I)*(THGB(I)-THX(I))
          HFX(I)=AMAX1(HFX(I),-250.)
        ENDIF
  400 CONTINUE

      DO I=its,ite
         IF((XLAND(I)-1.5).GE.0)THEN
           ZL=ZNT(I)
         ELSE
           ZL=0.01
         ENDIF
         CHS(I)=UST(I)*KARMAN/DENOMQ(I)
!        GZ2OZ0(I)=ALOG(2./ZNT(I))
!        PSIM2(I)=-10.*GZ2OZ0(I)
!        PSIM2(I)=AMAX1(PSIM2(I),-10.)
!        PSIH2(I)=PSIM2(I)
         CQS2(I)=UST(I)*KARMAN/DENOMQ2(I)
         CHS2(I)=UST(I)*KARMAN/DENOMT2(I)
      ENDDO

  410 CONTINUE
!jdf
!     DO I=its,ite
!       IF(UST(I).GE.0.1) THEN
!         RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I))
!       ELSE
!         RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1)
!       ENDIF
!     ENDDO
!jdf

!
   END SUBROUTINE SFCLAY1D

!====================================================================
   SUBROUTINE sfclayinit( allowed_to_read )

   LOGICAL , INTENT(IN)      ::      allowed_to_read
   INTEGER                   ::      N
   REAL                      ::      ZOLN,X,Y

   DO N=0,1000
      ZOLN=-FLOAT(N)*0.01
      X=(1-16.*ZOLN)**0.25
      PSIMTB(N)=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))- &
                2.*ATAN(X)+2.*ATAN(1.)
      Y=(1-16*ZOLN)**0.5
      PSIHTB(N)=2*ALOG(0.5*(1+Y))
   ENDDO

   END SUBROUTINE sfclayinit

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

END MODULE module_sf_sfclay

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
On a cool, autumn morning in 1992, Steve Walter and Lucas Taber were riding their bikes to school along the Old Davis Highway bike path.  During their journey, they discovered that during the night someone had severely vandalized the Murder Burger billboard that stood along Interstate 80.  The evidence suggested that an animal rights activist group had taken a chain saw to the structure, leaving little more than a wooden pile of debris and graffiti that read: Meat is Murder.  Beneath the carnage , Taber noticed what looked like lettuce and tomato, barely visible in the rubble.  Taber and Walter probed further, eventually uncovering the wooden hamburger that had previously adorned the billboard.  In effort to preserve this valuable piece of Davis history, the two boys took the burger with them to Davis High School, where they proudly hung the relic on the wall of the The Hub office, to be enjoyed by future generations.  This is a True South Davis Tales true story.

Hey, does anyone know if this is still there? what happened to it? updates to historical oddities is necessary for awesomeness ~Users/DavidPoole

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module utils_mod

use constants_mod

interface kronproduct
    procedure kronproduct_matmat
    procedure kronproduct_diagMat
end interface

interface pp
    procedure vp
    procedure mp
    procedure bp
    procedure tp
    procedure tp5
end interface

contains

!! pure elemental function logmean(uL,uR) result(r)
!!     real(dp), intent(in)    :: uL,uR
!!     real(dp)                :: r
!!
!!     r = (uL-uR+TOL)/(log(uL)-log(uR)+TOL)
!! end function

pure function calc_dist(a,b) result(dist)

    real(dp), intent(in)    :: a(N_DIMS),b(N_DIMS)
    real(dp)                :: dist
    real(dp)                :: v(N_DIMS)

    v = b-a
    dist = SQRT(SUM(v*v))

end function

pure elemental function logmean(aL,aR) result(r)

    real, intent(in)         :: aL, aR
    real(dp)                     :: x,u,r
    real(dp),parameter           :: eps = 1.0e-02_dp

    x = al/ar
    u = (x*(x-2.0_dp)+1.0_dp)/(x*(x+2.0_dp)+1.0_dp)
    r = merge((aL+aR)*52.50_dp/(105.0_dp + u*(35.0_dp + u*(21.0_dp +u*15.0_dp))), (aL-aR)/log(x), u < eps)

end function

!! subroutine write_headers(outputfp, buflen, buffer)
!!
!!     character(len=*), intent(in)    :: outputfp
!!     integer,intent(in)                 :: buflen
!!     character(len=*), intent(in)    :: buffer(:)
!!
!!     integer :: ioUnit
!!     integer :: openStat
!!
!!     character(range(buflen)+2) :: countStr
!!
!!     OPEN(NEWUNIT  = ioUnit             , &
!!          FILE     = TRIM(outputfp)     , &
!!          FORM     = 'FORMATTED'        , &
!!          STATUS   = 'UNKNOWN'          , &
!!          POSITION = 'REWIND'           , &
!!          RECL     = 50000              , &
!!          IOSTAT   = openStat             )
!!
!!     IF(openStat.NE.0) THEN
!!         CALL abort('ERROR: cannot open '// TRIM(outputfp))
!!     END IF
!!
!!     !! write nr. of columns to string
!!     write(countStr,'(i0)') buflen+1
!!     WRITE(ioUnit, '('// trim(countStr) //'(A))') "# ", buffer(1:buflen)
!!     CLOSE(ioUnit)
!!
!! end subroutine
!!

! This is a simple function to search for an available unit.
! LUN_MIN and LUN_MAX define the range of possible LUNs to check.
! The UNIT value is returned by the function, and also by the optional
! argument. This allows the function to be used directly in an OPEN
! statement, and optionally save the result in a local variable.
! If no units are available, -1 is returned.
function utils_get_newunit(unit) result(newunit)

    integer, intent(out), optional :: unit

    integer, parameter :: LUN_MIN=10, LUN_MAX=1000
    logical :: opened
    integer :: lun

    integer :: newunit

    newunit=-1

    do lun=LUN_MIN,LUN_MAX
        inquire(unit=lun,opened=opened)
        if (.not. opened) then
            newunit=lun
            exit
        end if
    end do

    if (present(unit)) unit = newunit

end function

subroutine write2file(outputfp, buflen, buffer)

    character(len=*), intent(in)    :: outputfp
    integer,intent(in) :: buflen
    real(dp),intent(in) :: buffer(:)

    integer :: ioUnit
    integer :: openStat

    character(range(buflen)+2) :: countStr

    !OPEN(NEWUNIT  = ioUnit             , &
    OPEN(UNIT     = utils_get_newunit(ioUnit), &
         FILE     = TRIM(outputfp)     , &
         FORM     = 'FORMATTED'        , &
         STATUS   = 'UNKNOWN'          , &
         POSITION = 'APPEND'           , &
         RECL     = 50000              , &
         IOSTAT   = openStat             )

    IF(openStat.NE.0) THEN
        !CALL abort('ERROR: cannot open '// TRIM(outputfp))
        write (*,*) 'ERROR: cannot open '// TRIM(outputfp)
        return
    END IF

    !! write nr. of columns to string
    write(countStr,'(i0)') buflen
    WRITE(ioUnit, '('// trim(countStr) //'(ES32.23,1x))') buffer(1:buflen)
    CLOSE(ioUnit)

end subroutine write2file

function incr(n)

    implicit none

    integer, intent(inout) :: n
    integer :: incr

    n = n + 1
    incr = n

end function

pure function krondelta(i,j) result(k)

    integer, intent(in) :: i,j
    integer :: k

    if (i == j) then
        k = 1
    else
        k = 0
    end if

end function

pure function isintegerdivisible(a,b) result(c)

    integer, intent(in) :: a,b
    integer :: c

    if (mod(a,b) == 0) then
        c = 1
    else
        c = 0
    end if

end function

function kronproduct_matmat(A,B) result(C)

    real(dp), intent(in) :: A(:,:)
    real(dp), intent(in) :: B(:,:)

    real(dp) :: C(size(A,1)*size(B,1),size(A,2)*size(B,2))

    integer     :: i,j
    integer     :: bn(2),cn(2)

    bn = shape(B)
    cn = shape(A)*shape(B)

    do i = 1,cn(1); do j = 1,cn(2)
        C(i,j) = A((i-1)/bn(1)+1,(j-1)/bn(2)+1) * B(mod(i-1,bn(1))+1, mod(j-1,bn(2))+1)
    end do; end do

end function

function kronproduct_diagMat(A,B) result(C)

    real(dp), intent(in) :: A(:)
    real(dp), intent(in) :: B(:,:)

    real(dp) :: C(size(A)*size(B,1),size(A)*size(B,2))

    integer     :: i,j, iA, jA
    integer     :: bn(2),cn(2)

    bn = shape(B)
    cn = size(A)*shape(B)

    do i = 1,cn(1); do j = 1,cn(2)
        iA = (i-1)/bn(1)+1
        jA = (j-1)/bn(2)+1

        if (iA == jA) then
            C(i,j) = A(iA) * B(mod(i-1,bn(1))+1, mod(j-1,bn(2))+1)
        else
            C(i,j) = 0.0_dp
        end if
    end do; end do

end function

function eyeproduct(N,B) result(C)

    integer, intent(in)     :: N
    real(dp), intent(in)    :: B(:,:)

    real(dp) :: C(N*size(B,1),N*size(B,2))

    integer     :: i,j, iA, jA
    integer     :: bn(2),cn(2)

    bn = shape(B)
    cn = N*shape(B)

    C = 0.0_dp

    do i = 1,cn(1); do j = 1,cn(2)
        iA = (i-1)/bn(1)+1
        jA = (j-1)/bn(2)+1

        if (iA == jA) then
            C(i,j) = B(mod(i-1,bn(1))+1, mod(j-1,bn(2))+1)
        end if
    end do; end do

end function

pure function calc_index_window(nelems,nranks,rankid) result(window)

    integer(kind=8), intent(in) :: nelems
    integer(kind=4), intent(in) :: nranks,rankid
    integer(kind=8)             :: window(2)

    integer(kind=8) :: a,b,NelemsPerRank,NelemsResidual

    NelemsPerRank = nelems / nranks
    NelemsResidual = mod(nelems,INT(nranks,kind=8))

    if (rankid < NelemsResidual) then
        a = rankid * NelemsPerRank + mod(INT(rankid,kind=8),NelemsResidual) + 1
        b = a + NelemsPerRank
    else
        a = rankid * NelemsPerRank + NelemsResidual + 1
        b = a + NelemsPerRank - 1
    end if

    window = (/a,b/)

end function

pure elemental function vmin(a,b) result(c)

    real(dp), intent(in)    :: a,b
    real(dp)                :: c

    c = min(a,b)

end function

pure elemental function vmax(a,b) result(c)

    real(dp), intent(in)    :: a,b
    real(dp)                :: c

    c = max(a,b)

end function

function tostr(num) result(str)

    integer, intent(in) :: num
    character(len=16)   :: str

    write(str,'(i0)') num

end function

subroutine tp5(title,blk)

    character(len=*), intent(in)    :: title
    real(dp), intent(in)            :: blk(:,:,:,:,:)

    integer :: i,k,l,q

    do q = 1,size(blk,5)
        do l = 1,size(blk,4)
            do k = 1,size(blk,3)
                write (*,'((a12),(a6),(i3),(a6),(i3),(a6),(i3))') trim(title),' | q = ', q, ' | w = ', l, ' | z = ', k
                do i = 1,size(blk,1)
                    write (*,'(99(ES12.4))') blk(i,:,k,l,q)
                end do
                write (*,*)
            end do
            write (*,*)
        end do
        write (*,*)
    end do
    write (*,*)

end subroutine

subroutine tp(title,blk)

    character(len=*), intent(in)    :: title
    real(dp), intent(in)            :: blk(:,:,:,:)

    integer :: i,k,l

    do l = 1,size(blk,4)
        do k = 1,size(blk,3)
            write (*,'((a12),(a6),(i3),(a6),(i3))') trim(title),' | w = ', l, ' | z = ', k
            do i = 1,size(blk,1)
                write (*,'(99(ES12.4))') blk(i,:,k,l)
            end do
            write (*,*)
        end do
        write (*,*)
    end do
    write (*,*)

end subroutine

subroutine bp(title,blk)

    character(len=*), intent(in)    :: title
    real(dp), intent(in)            :: blk(:,:,:)

    integer :: i,k

    do k = 1,size(blk,3)
        write (*,*) trim(title) // ' | z = ' // tostr(k)
        do i = 1,size(blk,1)
            write (*,'(99(ES12.4))') blk(i,:,k)
        end do
        write (*,*)
    end do
    write (*,*)

end subroutine

subroutine mp(title,mat)

    character(len=*), intent(in)    :: title
    real(dp), intent(in)            :: mat(:,:)

    integer :: i

    write (*,*) trim(title)
    do i = 1,size(mat,1)
        write (*,'(99(ES12.4))') mat(i,:)
    end do
    write (*,*)

end subroutine

subroutine vp(title,vec)

    character(len=*), intent(in)    :: title
    real(dp), intent(in)            :: vec(:)

    integer :: i

    write (*,*) trim(title)
    write (*,'(99(ES12.4))') vec
    write (*,*)

end subroutine

end module

max_stars_repo_name: jmark/nemo

max_issues_repo_name: jmark/nemo

max_forks_repo_name: jmark/nemo


text:
      SUBROUTINE ZERO(LX,X)
      DIMENSION X(LX)
      IF( LX .LE. 0 ) RETURN
      DO I = 1, LX
         X(I) = 0.
      ENDDO
      RETURN
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
      PROGRAM STOCH
      USE MOD_ZA  ! HYCOM array I/O interface
      IMPLICIT NONE
C
C     STOCHASTIC ARRAYS.
C
      INTEGER, ALLOCATABLE :: MSK(:,:)
      REAL*4,  ALLOCATABLE :: PLON(:,:),PLAT(:,:)
      REAL*4,  ALLOCATABLE :: VPM(:,:,:),TPM(:,:,:)
C
      LOGICAL   LARCTIC
      CHARACTER PREAMBL(5)*79
C
C     NAMELIST.
C
      CHARACTER*40     CTITLE
      NAMELIST/AFTITL/ CTITLE
C
      REAL*8           FSTART,DTI,SINC,SSTART,WSTART,TSTART,TMAX
      REAL*4           STFR_SPC(3),STFR_TIM(3),STFR_RMS(3)
      NAMELIST/AFTIME/ FSTART,DTI,SINC,SSTART,WSTART,TSTART,TMAX,
     +                 STFR_SPC,   STFR_TIM,   STFR_RMS
C
C**********
C*
C 1)  CREATE T,S,U,V MODEL GRID STOCHASTIC SURFACE FORCING
C      FILES SUITABLE FOR INPUT TO HYCOM OVER THE GIVEN REGION.
C
C 2)  NAMELIST INPUT:
C
C     /AFTITL/
C        CTITLE - ONE (40-CHARACTER) LINE TITLE.
C
C     /AFTIME/
C        FSTART - TIME OF FLUX START, IGNORED HERE         (DAYS)
C        WSTART - TIME OF WIND START, IGNORED HERE         (DAYS)
C        DTI    - NOMINAL TIME STEP FOR STOCHASTIC FORCING (SECS)
C        SINC   - TIME INCREMENT BETWEEN OUTPUT FIELDS     (DAYS)
C        SSTART - TIME OF STOCHASTIC FORCING               (DAYS)
C        TMAX   - TIME OF END   OF CURRENT INTEGRATION     (DAYS)
C        TSTART - TIME OF START OF CURRENT INTEGRATION     (DAYS)
C
C        STFR_SPC - HOR SPATIAL SCALE FOR STOCHASTIC FORCING (M)
C                    (1:T;2:S;3:UV)
C        STFR_TIM - TEMPORAL    SCALE FOR STOCHASTIC FORCING (HR)
C                    (1:T;2:S;3:UV)
C        STFR_RMS - RMS AMPLITUDE     FOR STOCHASTIC FORCING (UNITS/HR)
C                    (1:T;2:S;3:UV)
C
C     NAMELIST /AFTIME/ IS PATTERNED AFTER /XXTIME/ SO THAT THE
C      MODEL''S STANDARD AWK-BASED RUN SCRIPT CUSTOMIZER CAN ALSO
C      WORK FOR THE FLUX GENERATION SCRIPT.  IN PARTICULAR, 'WSTART'
C      AND 'FSTART' ARE READ IN, BUT NOT USED.
C
C 4)  INPUT:
C        ON UNIT  5:    NAMELIST /AFTITL/, /AFTIME/
C     OUTPUT:
C        ON UNIT 10:    UNFORMATTED MODEL T STOCASTIC FILE, SEE (6).
C        ON UNIT 11:    UNFORMATTED MODEL S STOCASTIC FILE, SEE (6).
C        ON UNIT 12:    UNFORMATTED MODEL U STOCASTIC FILE, SEE (6).
C        ON UNIT 13:    UNFORMATTED MODEL V STOCASTIC FILE, SEE (6).
C
C 6)  THE OUTPUT FIELDS ARE AT EVERY GRID POINT OF THE MODEL'S
C     'P' GRID.  ARRAY SIZE IS 'IDM' BY 'JDM'.
C
C 7)  SEVERAL STATISTICS ARE WRITTEN OUT IN ORDER TO CHECK THE
C      INTERPOLATION BETWEEN VARIOUS MACHINES.  MIN, MAX, MEAN AND RMS
C      OF THE ENTIRE BASIN ARE OUTPUT FOR kPAR FOR EACH
C      RECORD.  NOTE HOWEVER THAT THESE VALUES MAY NOT REPRESENT THE
C      STATISTICS OF THE FLUXS AS SEEN BY THE MODEL, IF THE INPUT FLUX
C      DATA HAS NON-REALISTIC VALUES OVER LAND.  IT IS UP TO THE USER
C      TO CHECK THE LOG FILES FOR CONSISTENCY BETWEEN MACHINES.
C
C 8)  ALAN J. WALLCRAFT,  NRL, SEPTEMBER 2017.
C      BASED ON NCOM VERSIONS BY PAUL MARTIN AND CLARK ROWLEY.
C*
C**********
C
      EXTERNAL AVERMS,MINMAX
C
      REAL*4     ZERO,RADIAN
      PARAMETER (ZERO=0.0, RADIAN=57.2957795)
C
      INTEGER NREC,MREC,NSTEP
C
      CHARACTER*80 CLINE
      REAL*8  WDAY8
      REAL*4  HMINA,HMINB,HMAXA,HMAXB
C
      INTEGER I,INIT,J,KREC,KSTEP
      REAL*4  FDY,WYR,DXREF,DTI4,
     +        XMIN,XMAX,XAVE,XRMS
C
C --- MODEL ARRAYS.
C
      CALL XCSPMD  !define idm,jdm
      ALLOCATE(  MSK(IDM,JDM),
     +          PLON(IDM,JDM),
     +          PLAT(IDM,JDM),
     +           TPM(IDM,JDM,2),
     +           VPM(IDM,JDM,2) )
C
C     NAMELIST INPUT.
C
      CALL ZHOPEN(6, 'FORMATTED', 'UNKNOWN', 0)
C
      CTITLE = ' '
      WRITE(6,*) 'READING /AFTITL/'
      CALL ZHFLSH(6)
      READ( 5,AFTITL)
      WRITE(6,AFTITL)
C
      FSTART = 0.0
      DTI    = 0.0
      SINC   = 0.0
      SSTART = 0.0
      WSTART = 0.0
      TSTART = 0.0
      TMAX   = 0.0
      WRITE(6,*) 'READING /AFTIME/'
      CALL ZHFLSH(6)
      READ( 5,AFTIME)
      IF     (SINC.GT.0.03D0) THEN
        SINC = NINT(SINC*24.D0)/24.D0  !assume multiple of 1 hour
      ENDIF
      WRITE(6,AFTIME)
      WRITE(6,*)
      CALL ZHFLSH(6)
C
      IF     (SSTART.LT.TSTART) THEN
        WRITE(6,*) 'ERROR - SSTART MUST BE AT LEAST TSTART'
        WRITE(6,*)
        CALL ZHFLSH(6)
        STOP
      ENDIF
      IF     (SINC.LE.0.0) THEN
        WRITE(6,*) 'ERROR - SINC MUST BE POSITIVE'
        WRITE(6,*)
        CALL ZHFLSH(6)
        STOP
      ENDIF
      IF     (DTI.LE.0.0) THEN
        WRITE(6,*) 'ERROR - DTI MUST BE POSITIVE'
        WRITE(6,*)
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
C     GRID INPUT.
C
      CALL ZAIOST
C
      CALL ZHOPNC(21, 'regional.grid.b', 'FORMATTED', 'OLD', 0)
      CALL ZAIOPF(    'regional.grid.a', 'OLD', 21)
C
      READ(21,*) ! skip idm
      READ(21,*) ! skip jdm
      READ(21,*) ! skip mapflg
      READ(21,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ (CLINE(I+1:),*)   HMINB,HMAXB
      CALL ZAIORD(PLON,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (plon):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      READ(21,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ (CLINE(I+1:),*)   HMINB,HMAXB
      CALL ZAIORD(PLAT,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (plat):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      CLOSE(UNIT=21)
      CALL ZAIOCL(21)
C
C     CHECK FOR AN ARCTIC BIPOLAR PATCH
C
      LARCTIC = PLON(3*IDM/4+1,JDM-1).EQ.PLON(IDM/4,JDM) .AND.
     &          PLAT(3*IDM/4+1,JDM-1).EQ.PLAT(IDM/4,JDM)
C
      WRITE(6,*)
      WRITE(6,*) 'LARCTIC = ',LARCTIC
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     INITIALIZE OUTPUT.
C
      WRITE(6,6000) 'OUTPUT:',CTITLE
      CALL ZHFLSH(6)
C
      CALL ZAIOPN('NEW', 10)
      CALL ZAIOPN('NEW', 11)
      CALL ZAIOPN('NEW', 12)
      CALL ZAIOPN('NEW', 13)
C
      CALL ZHOPEN(10, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(11, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(12, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(13, 'FORMATTED', 'NEW', 0)
C
      PREAMBL(1) = CTITLE
      PREAMBL(3) = ' '
      PREAMBL(4) = ' '
      WRITE(PREAMBL(5),'(A,2I5)')
     +        'i/jdm =',
     +       IDM,JDM
      WRITE(PREAMBL(2),'(A,F10.2,A,A,F10.2,A,A,F10.6,A)')
     +        'stfr_spc =',STFR_SPC(1),' m   ',
     +        'stfr_tim =',STFR_TIM(1),' hr   ',
     +        'stfr_rms =',STFR_RMS(1),' degC/hr'
      WRITE(10,4101) PREAMBL
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      WRITE(PREAMBL(2),'(A,F10.2,A,A,F10.2,A,A,F10.6,A)')
     +        'stfr_spc =',STFR_SPC(2),' m   ',
     +        'stfr_tim =',STFR_TIM(2),' hr   ',
     +        'stfr_rms =',STFR_RMS(2),' psu/hr'
      WRITE(11,4101) PREAMBL
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      WRITE(PREAMBL(2),'(A,F10.2,A,A,F10.2,A,A,F10.6,A)')
     +        'stfr_spc =',STFR_SPC(3),' m   ',
     +        'stfr_tim =',STFR_TIM(3),' hr   ',
     +        'stfr_rms =',STFR_RMS(3),' m/s/hr'
      PREAMBL(3) = 'on u-grid'
      WRITE(12,4101) PREAMBL
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      PREAMBL(3) = 'on v-grid'
      WRITE(13,4101) PREAMBL
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      WRITE(6,*)
C
C     PROCESS ALL THE STOCHASTIC RECORDS.
C
      NREC  = NINT( (TMAX   - TSTART) / SINC ) + 1
      MREC  = NINT( (SSTART - TSTART) / SINC ) + 1
      NSTEP = NINT((SINC*86400.0d0)/DTI)
      DXREF = (PLAT(IDM/2,JDM/2+1)-PLAT(IDM/2+1,JDM/2))*111200.0
      DTI4  = DTI
C
      WRITE(6,'(a,3i4,f12.3)')
     &  'nrec,mrec,nstep,dxref = ',nrec,mrec,nstep,dxref
      WRITE(6,*)
C
      TPM(:,:,1) = 0.0
      TPM(:,:,2) = 0.0
      VPM(:,:,1) = 0.0
      VPM(:,:,2) = 0.0
C
      DO 810 KREC= 1,NREC
C
        WDAY8 = TSTART + (KREC-1)*SINC
C
        IF     (KREC.EQ.MREC) THEN
C
C         FORM INITIAL STOCHASTIC FIELDS.
C
          init = 0
          call stoch_r(init,idm,jdm,2,dxref,dti4,
     &                 stfr_spc(1),stfr_tim(1),stfr_rms(1),tpm)
          init = 0
          call stoch_v(init,idm,jdm,  dxref,dti4,
     &                 stfr_spc(3),stfr_tim(3),stfr_rms(3),vpm(1,1,1))
          init = 1
          call stoch_v(init,idm,jdm,  dxref,dti4,
     &                 stfr_spc(3),stfr_tim(3),stfr_rms(3),vpm(1,1,2))
        ELSEIF (KREC.GT.MREC) THEN
C
C         TIME ADVANCE STOCHASTIC FIELDS.
C
          init = 1
          do kstep= 1,nstep
            call stoch_r(init,idm,jdm,2,dxref,dti4,
     &                   stfr_spc(1),stfr_tim(1),stfr_rms(1),tpm)
            call stoch_v(init,idm,jdm,  dxref,dti4,
     &                   stfr_spc(3),stfr_tim(3),stfr_rms(3),vpm(1,1,1))
            call stoch_v(init,idm,jdm,  dxref,dti4,
     &                   stfr_spc(3),stfr_tim(3),stfr_rms(3),vpm(1,1,2))
          enddo !istep
        ENDIF !krec
C
C       WRITE OUT STATISTICS.
C
        CALL ARCUPD(TPM(1,1,1),IDM,JDM, LARCTIC,.TRUE.)
        CALL MINMAX(TPM(1,1,1),IDM,JDM, XMIN,XMAX)
        CALL AVERMS(TPM(1,1,1),IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'stocht', XMIN,XMAX,XAVE,XRMS
C
        CALL ARCUPD(TPM(1,1,2),IDM,JDM, LARCTIC,.TRUE.)
        CALL MINMAX(TPM(1,1,2),IDM,JDM, XMIN,XMAX)
        CALL AVERMS(TPM(1,1,2),IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'stochs', XMIN,XMAX,XAVE,XRMS
C
        CALL ARCUPD(VPM(1,1,1),IDM,JDM, LARCTIC,.TRUE.)
        CALL MINMAX(VPM(1,1,1),IDM,JDM, XMIN,XMAX)
        CALL AVERMS(VPM(1,1,1),IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'stochu', XMIN,XMAX,XAVE,XRMS
C
        CALL ARCUPD(VPM(1,1,2),IDM,JDM, LARCTIC,.TRUE.)
        CALL MINMAX(VPM(1,1,2),IDM,JDM, XMIN,XMAX)
        CALL AVERMS(VPM(1,1,2),IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'stochv', XMIN,XMAX,XAVE,XRMS
C
C       WRITE OUT HYCOM FIELDS.
C
        CALL ZAIOWR(TPM(1,1,1),MSK,.FALSE., XMIN,XMAX, 10, .FALSE.)
        WRITE(10,4122) 'stocht',WDAY8,SINC,XMIN,XMAX
        CALL ZHFLSH(10)
C
        CALL ZAIOWR(TPM(1,1,2),MSK,.FALSE., XMIN,XMAX, 11, .FALSE.)
        WRITE(11,4122) 'stochs',WDAY8,SINC,XMIN,XMAX
        CALL ZHFLSH(11)
C
        CALL ZAIOWR(VPM(1,1,1),MSK,.FALSE., XMIN,XMAX, 12, .FALSE.)
        WRITE(12,4122) 'stochu',WDAY8,SINC,XMIN,XMAX
        CALL ZHFLSH(12)
C
        CALL ZAIOWR(VPM(1,1,2),MSK,.FALSE., XMIN,XMAX, 13, .FALSE.)
        WRITE(13,4122) 'stochv',WDAY8,SINC,XMIN,XMAX
        CALL ZHFLSH(13)
C
        CALL WNDAY(WDAY8, WYR,FDY)
        WRITE(6,6350) KREC,WDAY8,FDY,NINT(WYR)
        CALL ZHFLSH(6)
  810 CONTINUE
C
      CALL ZAIOCL(10)
      CALL ZAIOCL(11)
      CALL ZAIOCL(12)
      CALL ZAIOCL(13)
      CLOSE( UNIT=10)
      CLOSE( UNIT=11)
      CLOSE( UNIT=12)
      CLOSE( UNIT=13)
C
C     SUMMARY.
C
      CALL WNDAY(TSTART, WYR,FDY)
      WRITE(6,6450) NREC,FDY,NINT(WYR),WDAY8-TSTART
      CALL ZHFLSH(6)
      STOP
C
 4101 FORMAT(A79)
 4112 FORMAT(2X,A,': range = ',1P2E16.7)
 4122 FORMAT(2X,A,': day,span,range =',F12.5,F10.6,1P2E16.7)
 6000 FORMAT(1X,A,2X,A40 //)
 6200 FORMAT(/ 1X,'MIN,MAX I COORDS = ',F8.2,',',F8.2
     +       / 1X,'MIN,MAX J COORDS = ',F8.2,',',F8.2 /)
 6350 FORMAT(10X,'WRITING STOC RECORD',I5,
     +           '    SDAY =',F10.3,
     +            '  SDATE =',F8.3,'/',I4 /)
 6450 FORMAT(I5,' STOC RECORDS STARTING ON',F8.3,'/',I4,
     +   ' COVERING',F12.5,' DAYS')
 8100 FORMAT(1X,A,': MIN=',F12.8,' MAX=',F12.8,
     +             ' AVE=',F12.8,' RMS=',F12.8)
C     END OF PROGRAM STOCH
      END
      SUBROUTINE WNDAY(WDAY, YEAR,DAY)
      IMPLICIT NONE
      REAL*8 WDAY
      REAL*4 YEAR,DAY
C
C**********
C*
C  1) CONVERT 'FLUX DAY' INTO JULIAN DAY AND YEAR.
C
C  2) THE 'FLUX DAY' IS THE NUMBER OF DAYS SINCE 001/1901 (WHICH IS
C      FLUX DAY 1.0).
C     FOR EXAMPLE:
C      A) YEAR=1901.0 AND DAY=1.0, REPRESENTS 0000Z HRS ON 001/1901
C         SO WDAY WOULD BE 1.0.
C      B) YEAR=1901.0 AND DAY=2.5, REPRESENTS 1200Z HRS ON 002/1901
C         SO WDAY WOULD BE 2.5.
C     YEAR MUST BE NO LESS THAN 1901.0, AND NO GREATER THAN 2099.0.
C     NOTE THAT YEAR 2000 IS A LEAP YEAR (BUT 1900 AND 2100 ARE NOT).
C
C  3) ALAN J. WALLCRAFT, PLANNING SYSTEMS INC., FEBRUARY 1993.
C*
C**********
C
      INTEGER IYR,NLEAP
      REAL*8  WDAY1
C
C     FIND THE RIGHT YEAR.
C
      IYR   = (WDAY-1.0)/365.25
      NLEAP = IYR/4
      WDAY1 = 365.0*IYR + NLEAP + 1.0
      DAY   = WDAY - WDAY1 + 1.0
      IF     (WDAY1.GT.WDAY) THEN
        IYR   = IYR - 1
      ELSEIF (DAY.GE.367.0) THEN
        IYR   = IYR + 1
      ELSEIF (DAY.GE.366.0 .AND. MOD(IYR,4).NE.3) THEN
        IYR   = IYR + 1
      ENDIF
      NLEAP = IYR/4
      WDAY1 = 365.0*IYR + NLEAP + 1.0
C
C     RETURN YEAR AND JULIAN DAY.
C
      YEAR = 1901 + IYR
      DAY  = WDAY - WDAY1 + 1.0
      RETURN
C     END OF WNDAY.
      END

      subroutine stoch_v(init,n,m,delx,dti,
     &                   stfv_ss,stfv_ts,stfv_rms, stfv)
c
c  created 2014-04-07, Paul J Martin, NRL.
      implicit none
c
c  declare passed variables.
      integer init,n,m
      real*4  delx,dti
      real*4  stfv_ss,stfv_ts,stfv_rms
      real*4  stfv(n,m)
c
c  compute/update stochastic forcing fields on 1st call for new time step.
*     write(6,*) 'v: init,n,m    =',init,n,m
*     write(6,*) 'v: delx,dti    =',delx,dti
*     call zhflsh(6)
*
      call stoch_field(init,delx,dti,n,m,1,
     &  stfv_ss,stfv_ts,stfv_rms,stfv)
c
      return
      end

      subroutine stoch_r(init,n,m,nr,delx,dti,
     &                   stfr_ss,stfr_ts,stfr_rms, stfr)
c
c  created 2014-04-07, Paul J Martin, NRL.
      implicit none
c
c  declare passed variables.
      integer init,n,m,nr
      real*4  delx,dti
      real*4  stfr_ss(nr),stfr_ts(nr),stfr_rms(nr)
      real*4  stfr(n,m,nr)
c
c  compute/update stochastic forcing fields on 1st call for new time step.
*     write(6,*) 'r: init,n,m,nr =',init,n,m,nr
*     write(6,*) 'r: delx,dti    =',delx,dti
*     call zhflsh(6)
*
      call stoch_field(init,delx,dti,n,m,nr,
     &  stfr_ss,stfr_ts,stfr_rms,stfr)
c
      return
      end

      subroutine stoch_field(init,delx,dti,n,m,nstf,
     &  stf_ss,stf_ts,stf_rms,stf)
c  subroutine to compute stochastic forcing field.
c      init   = integer flag to denote initial calculation of stochastic
c               forcing field (when =0).
c      delx   = approximate grid spacing in m.
c      dti    = time step
c      n,m    = horizontal dimensions of local tile.
c      l      = number of vertical layers +1.
c      nstf   = number of stochastic forcing fields being computed.
c      stf_ss = horizontal spatial scale for stochastic forcing (+m).
c      stf_ts = temporal scale for stochastic forcing (hr).
c      stf_rms= rms amplitude for stochastic forcing (units/hr).
c      stf    = returned stochastic forcing field.
c  created 2014-04-04, Paul J Martin, NRL.
      implicit none
c
c  declare passed variables.
      integer init,n,m,nstf
      real*4  delx,dti
      real*4  stf_ss(nstf),stf_ts(nstf),stf_rms(nstf)
      real*4  stf(n,m,nstf)
c
c  declare local temporary variables.
      integer i,j,ir,iflt,nflt
      real*4  a,b,c,d
      real*4, allocatable :: r1(:,:),r2(:,:)
      real*8  a8,b8,c8,d8
      integer indgrd
      real*4  timed,dum,amult,cint,vscale
c
c  diagnostic printout.
      if (init .eq. 0) then
        call init_random_seed
      endif
c
c  allocate needed local arrays.
      allocate( r1(n,m) )
      allocate( r2(n,m) )
c
c  loop through all the stochastic forcing fields being computed.
      do ir=1,nstf
*       write(6,*) 'ir,delx,dti =',ir,delx,dti
*       write(6,*) 'ir,stf      =',ir,stf_ss(ir),stf_ts(ir),stf_rms(ir)
*       call zhflsh(6)
c
c  calc random 2D field (note:  random_number is a system routine).
c  this routine returns uniformly-distributed random numbers from 0 to 1.
      call random_number(r1)
      r2 = 0.
c
c  average over 3x3 boxes to give an approximate normal distribution.
c  this has small tails because of small number of values (9) averaged.
c  the value "a" is subtracted to provide an ~ mean value of zero.
c  note:  value of "a" used here may need slight tweek to give mean =0.!!
      a=4.5
      do j=1,m
        do i=2,n-1
          r2(i,j)=r1(i-1,j)+r1(i,j)+r1(i+1,j)
        enddo
        r2(1,j)=r2(2,j)
        r2(n,j)=r2(n-1,j)
      enddo
      do i=1,n
        do j=2,m-1
          r1(i,j)=(r2(i,j-1)+r2(i,j)+r2(i,j+1)) - a
        enddo
        r1(i,1)=r1(i,2)
        r1(i,m)=r1(i,m-1)
      enddo
c
c  use spatial filtering to achieve desired spatial correlation scale.
c
c  estimate number of passes of spatial filter needed.
c  note:  this formula may need some refinement.!!
      nflt=int(0.300*(stf_ss(ir)/delx)**2)
c
c  apply spatial filter.
      do iflt=1,nflt
c  apply filter in x-direction.
        do j=1,m
          do i=2,n-1
            r2(i,j)=(r1(i-1,j)+r1(i+1,j)) + 2.*r1(i,j)
          enddo
          r2(1,j)=r2(2,j)
          r2(n,j)=r2(n-1,j)
        enddo
c  apply filter in y-direction.
        a=1./16.0
        do i=1,n
          do j=2,m-1
            r1(i,j)=( (r2(i,j-1)+r2(i,j+1)) + 2.*r2(i,j) )*a
          enddo
          r1(i,1)=r1(i,2)
          r1(i,m)=r1(i,m-1)
        enddo
      enddo
c
c  scale the amplitude of the stochastic forcing field.
c  note:  this formula may need some tuning.!!
      a=stf_rms(ir)*(nflt+1)**0.465/0.86366
      do j=1,m
        do i=1,n
          r1(i,j)=r1(i,j)*a
        enddo
      enddo
c
c  apply temporal correlation.
c  note that the particular temporal weightings (a,b) used here may seem
c  a bit strange, but these will maintain a fairly constant amplitude
c  for the computed stochastic forcing field (stf) with time (which is
c  desireable).
      a=dti/(3600.0*stf_ts(ir))
      c=exp(-a)
      d=sqrt(1. - exp(-2.*a))
      if (init .eq. 0) then
        do j=1,m
          do i=1,n
            stf(i,j,ir)=r1(i,j)
          enddo
        enddo
      else
*       write(6,*) 'ir,a,c,d    =',ir,a,c,d
*       call zhflsh(6)
        do j=1,m
          do i=1,n
            stf(i,j,ir)=c*stf(i,j,ir)+d*r1(i,j)
          enddo
        enddo
      endif
c
      enddo
c
c  deallocate local temporary arrays.
      deallocate(r1,r2)
c
      return
      end

      subroutine init_random_seed
c  subroutine to get random seed values for random number generator
c  from file = /dev/urandom .
c  note that the intel workstations return n=34 and 10-digit integer
c  seed values with both + and - signs.
c  created 2014-06-27, Paul J Martin, NRL.
      implicit none
c
      integer n,istat,i,datetime(8)
      integer, allocatable :: iseed(:)
c
      call random_seed(size = n)
*     write(6,*) 'rs: n =',n
*     call zhflsh(6)
      allocate( iseed(n) )
c
      call date_and_time(values=datetime)
      do i= 1,min(8,n)
        iseed(i) = datetime(9-i)
      enddo
      do i= 9,n
        iseed(i) = i
      enddo
*     write(6,*) 'dtime =',datetime
      write(6,*) 'iseed =',iseed(1:n)
      call zhflsh(6)
      call random_seed(put = iseed)
c
      deallocate( iseed )
      return
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
      PROGRAM TESTPCG
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C   SOLVER : PCG
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     -- Déclarations ---------------------------------------
      IMPLICIT REAL*8(A-H,O-Z)

      INTEGER N, IA, JA, NELEM, ierr,it_max,iwk,iflag
      parameter(N = 100, NELEM = 3500,iwk=10000)

      DIMENSION A(NELEM), IA(N+1), JA(NELEM),
     #          ALU(iwk),JLU(iwk),JU(N+1),
     #          b(N), x(N), z(N),
     #          r(N), w(N), p(N),DR(N),JI(N),xref(N)
      LOGICAL JT(N)

C     -- Chargement de la matrice A et du vecteur b ---------
      CALL loadmat(A,IA,JA,b,xref,ierr)

C     -- Initialisation du vecteur x à 0.0d0 ----------------
      CALL vectinit(N,x)

C     -- Paramètres du GMRES --------------------------------
      epsilon = 1.0D-8
      it_max  = 1000

C     -- Création du préconditionneur -----------------------
      WRITE(*,*)' 1 : IC0'
      WRITE(*,*)' 2 : IC0  + M réduc. = RIC0'
      WRITE(*,*)' 3 : MIC0'
      WRITE(*,*)' 4 : MIC0 + M réduc. = RMIC0'
      WRITE(*,*)'precond ? (0 = non)'
      READ(*,*)iflag
      ierr= 0
C
      if (iflag.eq.1) then
         CALL IC0 (N,A,IA,JA,ALU,JU,JLU,JT,JI,DR,iwk,
     $             ierr)
      else if (iflag.eq.2) then
         CALL RIC0 (N,A,IA,JA,ALU,JU,JLU,JT,JI,DR,iwk,
     $             ierr)
      else if (iflag.eq.3) then
         WRITE(*,*)'alpha ?'
         READ(*,*)alpha
         CALL MIC0 (N,A,IA,JA,ALU,JU,JLU,JT,JI,DR,iwk,
     $             ierr,alpha)
      else if (iflag.eq.4) then
         WRITE(*,*)'alpha ?'
         READ(*,*)alpha
         CALL RMIC0 (N,A,IA,JA,ALU,JU,JLU,JT,JI,DR,iwk,
     $               ierr,alpha)
      endif
C     -- Appel de la sous-routine CG ---------------------

      CALL pcg(N, x, b, A, IA, JA, ALU,JU,JLU, epsilon,
     #         it_max,r, w, p, z, iflag)

C     -- Sauvegarde du résultat vers MATLAB ----------------

      OPEN (UNIT = 1, FILE = 's_gmres.m', STATUS='unknown')
      DO 100 i = 1, N
         WRITE(1,*)'x(',i,')=',x(i),';'
100   CONTINUE
      CLOSE (UNIT = 1)

C     -- Visualisation de la solution ----------------------

      WRITE(*,*)
      WRITE(*,*)' <ENTER> pour voir la solution'
      READ(*,*)
      CALL VectAff2(N, x, xref)

      END

ext: for

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
!    GNOMO 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 single_node
use general
use neighboring
use growth
use mitosis

contains

subroutine should_I_divide_single
integer ick,j,k,ii,jj,kk
real*8 a,b,c,s,sx,sy,sz,d

do ick=1,ncels
  s=0.0d0
  ii=cels(ick)%node(1)
  if(node(ii)%hold==1) cycle !hold nodes doesn't perform cell behaviors
  if(cels(ick)%ctipus<3)then
    jj=node(ii)%altre
    c=1-0.5*(node(ii)%diffe+node(jj)%diffe)
    do k=1,npag(nparam_per_node+2)
      kk=whonpag(nparam_per_node+2,k)
      if (gex(ii,kk)>0.0d0.or.gex(jj,kk)>0.0d0) then
        s=s+(gex(ii,kk)+gex(jj,kk))*gen(kk)%wa(nparam_per_node+2)*c
      end if
    end do
    s=s*0.5 !this way %fase is independent of cell size !>>>>Miquel2-12-13
  else
    c=1-node(ii)%diffe
    do k=1,npag(nparam_per_node+2)
      kk=whonpag(nparam_per_node+2,k)
      if (gex(ii,kk)>0.0d0) then
        s=s+gex(ii,kk)*gen(kk)%wa(nparam_per_node+2)*c
      end if
    end do
  end if
  cels(ick)%fase=cels(ick)%fase+s*delta
  if (cels(ick)%fase>=1.0d0) then
    if (ffu(25)==1) then  !>>> Is 18-4-15
      if (node(cels(ick)%node(1))%req>=mmae) then !>>> Is 18-4-15 !>> Miquel20-4-15
        call division_single_node(ick)!>>> Is 18-4-15
        cels(ick)%fase=0d0  !>>> Is 18-4-15
        cels(ncels)%fase=0d0  !>>> Is 18-4-15
      end if
    else
      call division_single_node(ick)    !>>> Is 18-4-15
      cels(ick)%fase=0d0    !>>> Is 18-4-15
      cels(ncels)%fase=0d0  !>>> Is 18-4-15
    end if   !>>> Is 18-4-15
  end if  !>>> Is 18-4-15
end do

end subroutine should_I_divide_single

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

subroutine division_single_node(celd)   ! IMPORTANT, THIS IMPLIES AN INELUDIBLE BIAS IF THE INITIAL CONDITIONS ARE REGULAR SINCE THE NEW CELL HAS TO ARISE IN A SPECIFIC DIRECTION
integer::celd,nnod,tipi,nnoda,nnodb,i,j,k,ii,jj,kk
real*8::a,b,c,d,ax,ay,az,bx,by,bz,cx,cy,cz,dotp,ix,iy,iz,nodda
integer,dimension(:)::nodea(cels(celd)%nunodes),nodeb(cels(celd)%nunodes)

  ii=cels(celd)%node(1) !the node
  tipi=node(ii)%tipus
  nodda=node(ii)%da
  !determining the physicial division plane (in relation with the neighboring cells)
  !just taking the longest vector connecting two of its neighbors
  if(tipi<3)then !epithelial : only looks for neighbors of the same %tipus
    aa=0
    bx=0;by=0;bz=0
    if(nneigh(ii)>2)then !>>Miquel20-3-15
      do i=1,nneigh(ii)-1
        jj=neigh(ii,i)
        if(node(jj)%tipus==tipi.and.(dneigh(ii,i)<nodda+node(jj)%da))then
          do j=i+1,nneigh(ii)
            kk=neigh(ii,j)
            if(node(kk)%tipus==tipi.and.dneigh(ii,j)<nodda+node(kk)%da)then
              a=node(kk)%x-node(jj)%x ; b=node(kk)%y-node(jj)%y ; c=node(kk)%z-node(jj)%z
              d=sqrt(a**2+b**2+c**2)
             if(d>aa)then
                aa=d
                bx=a ; by=b ; bz=c
              end if
            end if
          end do
        end if
      end do
    else
      call random_number(bx)
      call random_number(by)
      call random_number(bz)
      aa=sqrt(bx**2+by**2+bz**2)
    end if

    d=1/aa
    bx=bx*d ; by=by*d ; bz=bz*d
    jj=node(ii)%altre
    ax=node(jj)%x-node(ii)%x ; ay=node(jj)%y-node(ii)%y ; az=node(jj)%z-node(ii)%z
    d=1/sqrt(ax**2+ay**2+az**2)
    ax=ax*d;ay=ay*d;az=az*d
    dotp=ax*bx+ay*by+az*bz
    cels(celd)%hpolx=bx-ax*dotp ; cels(celd)%hpoly=by-ay*dotp ; cels(celd)%hpolz=bz-az*dotp    !physical vector
  else    !mesenchymal : looks for any neighbor
    aa=0
    bx=0;by=0;bz=0
    if(nneigh(ii)>1)then !>>Miquel20-3-15
      do i=1,nneigh(ii)-1
        jj=neigh(ii,i)
        if(dneigh(ii,i)<nodda+node(jj)%da)then
          do j=i+1,nneigh(ii)
            kk=neigh(ii,j)
            if(dneigh(ii,j)<nodda+node(kk)%da)then
              a=node(kk)%x-node(jj)%x ; b=node(kk)%y-node(jj)%y ; c=node(kk)%z-node(jj)%z
              d=sqrt(a**2+b**2+c**2)
              if(d>aa)then
                aa=d
                bx=a ; by=b ; bz=c
              end if
            end if
          end do
        end if
      end do
    else
      call random_number(bx)
      call random_number(by)
      call random_number(bz)
      aa=sqrt(bx**2+by**2+bz**2)
    end if
    cels(celd)%hpolx=bx ; cels(celd)%hpoly=by ; cels(celd)%hpolz=bz    !physical vector
  end if

  d=0.0d0
  do j=1,npag(nparam_per_node+11)    !number of genes affecting growth
    k=whonpag(nparam_per_node+11,j)  !which are those genes
    if(tipi<3)then
      d=d+(gex(ii,k)+gex(jj,k))*gen(k)%wa(nparam_per_node+11) !this is the differential of growth for the node
    else
      d=d+gex(ii,k)*gen(k)%wa(nparam_per_node+11) !this is the differential of growth for the node
    end if
  end do
  d=1d0/(1d0+d)  !ponderacion entre vector fisico i quimico
 !print*,"d ponderation",d,node(nd)%tipus
  !d=0d0
  !d    ; dependence of the gradient vector (how many it affects to polarization vector)
  !d=0  ; polarization vector comes only from its shape  (default mode)
  !d=1  ; polarization vector comes only from the gradient
  !d=0.5; polarization vector comes equally from both the gradient and the shape
  cx=((1-d)*cels(celd)%polx)+(d*cels(celd)%hpolx)   !vector resultante para division
  cy=((1-d)*cels(celd)%poly)+(d*cels(celd)%hpoly)
  cz=((1-d)*cels(celd)%polz)+(d*cels(celd)%hpolz)
  d=1/sqrt(cx**2+cy**2+cz**2)
  cx=cx*d ; cy=cy*d ; cz=cz*d
  ncels=ncels+1
  if(tipi<3)then
    ndepi=ndepi+2
    nd=nd+2
    nnodb=2
    nnoda=2
    !i=jj
    !jj=ii
    !ii=i
  else
    ndmes=ndmes+1
    nd=nd+1
    nnodb=1
    nnoda=1
  end if

  call addanode(ii)
  node(nd)%icel=ncels
!  node(nd)%icel=ncels!pfh17-3-15

  node(nd)%x=node(ii)%x+desmax*cx ; node(nd)%y=node(ii)%y+desmax*cy ; node(nd)%z=node(ii)%z+desmax*cz
 !print *,cx,cy,cz,cels(celd)%polx,cels(celd)%poly,cels(celd)%polz,"koop",desmax
  node(nd)%marge=node(ii)%marge !>>> Is 1-15
  nodeo(nd)%x=node(nd)%x ; nodeo(nd)%y=node(nd)%y ; nodeo(nd)%z=node(nd)%z !pfh18-3-15

  if (ffu(25)==0) then !>>> Is 18-4-15
    node(nd)%reqcr=node(ii)%reqcr !>>> Is 18-4-15
  else                            !>>> Is 18-4-15
    node(nd)%reqcr=reqmin         !>>> Is 18-4-15
  end if                          !>>> Is 18-4-15

  if (prop_noise==0d0) node(nd)%e=0d0 !pfh 20-3-15
  if(tipi<3)then
    node(nd-1)%icel=ncels
    node(nd-1)%x=node(jj)%x+desmax*cx ; node(nd-1)%y=node(jj)%y+desmax*cy ; node(nd-1)%z=node(jj)%z+desmax*cz
    if (ffu(25)==0) then !>>> Is 18-4-15
      node(nd-1)%reqcr=node(ii)%reqcr !>>> Is 18-4-15
    else                              !>>> Is 18-4-15
      node(nd-1)%reqcr=reqmin         !>>> Is 18-4-15
    end if                            !>>> Is 18-4-15
    ! node(nd)%marge=0 >>> Is 1-15
    node(nd-1)%marge=1-node(nd)%marge !>>> Is 1-15
    nodeo(nd-1)%x=node(nd-1)%x ; nodeo(nd-1)%y=node(nd-1)%y ; nodeo(nd-1)%z=node(nd-1)%z !pfh18-3-15
    if (prop_noise==0d0) node(nd-1)%e=0d0  !pfh 20-3-15
  else
    node(nd)%marge=0
  end if

  !cels(celd)%nunodes=nnoda
  cels(ncels)%nunodes=nnodb
  !deallocate(cels(celd)%node)
  !cels(celd)%nodela=nnoda+10
  cels(ncels)%nodela=nnodb
  allocate(cels(ncels)%node(nnodb))
  cels(celd)%nunodes=nnoda
  cels(celd)%node(:)=0
  cels(celd)%node(1)=ii

  cels(ncels)%node(:)=0
  cels(ncels)%node(1)=nd

  cels(ncels)%hpolx=0 ; cels(ncels)%hpoly=0 ; cels(ncels)%hpolz=0  !>>Miquel12-9-14
  cels(ncels)%polx=0 ; cels(ncels)%poly=0 ; cels(ncels)%polz=0     !>>Miquel12-9-14

  if(tipi<3)then
    cels(celd)%node(2)=node(ii)%altre
    cels(ncels)%node(2)=nd-1
    cels(ncels)%cex=node(nd-1)%x ; cels(ncels)%cey=node(nd-1)%y ; cels(ncels)%cez=node(nd-1)%z
  else
    cels(ncels)%cex=node(nd)%x ; cels(ncels)%cey=node(nd)%y ; cels(ncels)%cez=node(nd)%z
  end if
  cels(ncels)%ctipus=cels(celd)%ctipus
  cels(ncels)%maxsize_for_div=cels(celd)%maxsize_for_div!pfh-17-3-15
  cels(ncels)%minsize_for_div=cels(celd)%minsize_for_div!pfh-17-3-15
  cels(ncels)%temt=cels(celd)%temt			!pfh-17-3-15

  if (ncels>=ncals) call recels

end subroutine division_single_node

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

subroutine polarization_single  !it only works for diffusible molecules, should be implemented for membrane-tethered molecules
real*8::vx,vy,vz,ix,iy,iz,ivx,ivy,ivz,ax,ay,az,d,alfa
integer::iv,celi,i,j,k,jj,kk

  !print*,"POLARISATION_SINGLE"
  do i=1,nd
    if(node(i)%tipus==2)then !epithelial
      vx=0 ; vy=0 ; vz=0
      iv=node(i)%altre
      ix=node(i)%x ; iy=node(i)%y ; iz=node(i)%z
      ivx=node(iv)%x ; ivy=node(iv)%y ; ivz=node(iv)%z
      celi=node(i)%icel
      do j=1,nneigh(i)
        jj=neigh(i,j)
        if(node(jj)%tipus/=node(i)%tipus) cycle !only compare nodes of the same layer
        ax=node(jj)%x-ix ; ay=node(jj)%y-iy ; az=node(jj)%z-iz
        d=1/dneigh(i,j)
        do k=1,npag(nparam_per_node+8)
          kk=whonpag(nparam_per_node+8,k)
          if (gex(jj,kk)>0.0d0) then
            alfa=gex(jj,kk)*gen(kk)%wa(nparam_per_node+8)*d!*gen(kk)%diffu !proportional to the diffusion differential
            vx=vx+ax*alfa ; vy=vy+ay*alfa ; vz=vz+az*alfa
          else
            cycle
          end if                                                                         !just a bit different (for simplification)
        end do
      end do
      do j=1,nneigh(iv)
        jj=neigh(iv,j)
        if(node(jj)%tipus/=node(iv)%tipus) cycle !only compare nodes of the same layer
        ax=node(jj)%x-ivx ; ay=node(jj)%y-ivy ; az=node(jj)%z-ivz
        d=1/dneigh(iv,j)
        do k=1,npag(nparam_per_node+8)
          kk=whonpag(nparam_per_node+8,k)
          if (gex(jj,kk)>0.0d0) then
            alfa=gex(jj,kk)*gen(kk)%wa(nparam_per_node+8)*d!*gen(kk)%diffu !proportional to the diffusion differential
            vx=vx+ax*alfa ; vy=vy+ay*alfa ; vz=vz+az*alfa
          else
            cycle
          end if                                                                         !just a bit different (for simplification)
        end do
      end do
      if(vx==0.and.vy==0.and.vz==0)then
        d=0
      else
        d=1/sqrt(vx**2+vy**2+vz**2)
      end if
      cels(celi)%polx=vx*d ; cels(celi)%poly=vy*d ; cels(celi)%polz=vz*d
!      cels(celi)%polx=0.0d0 ; cels(celi)%poly=0.0d0 ; cels(celi)%polz=1.0d0 !ACHTUNG  ACHTUNG  ACHTUNG  ACHTUNG  ACHTUNG  ACHTUNG  ACHTUNG  ACHTUNG
    elseif(node(i)%tipus==3)then  !mesenchymal
      vx=0 ; vy=0 ; vz=0
      ix=node(i)%x ; iy=node(i)%y ; iz=node(i)%z
      celi=node(i)%icel
      do j=1,nneigh(i)
        jj=neigh(i,j)
        ax=node(jj)%x-ix ; ay=node(jj)%y-iy ; az=node(jj)%z-iz
        d=1/dneigh(i,j)
        do k=1,npag(nparam_per_node+8)
          kk=whonpag(nparam_per_node+8,k)
          if (gex(jj,kk)>0.0d0) then
            alfa=gex(jj,kk)*gen(kk)%wa(nparam_per_node+8)*d!*gen(kk)%diffu !proportional to the diffusion differential
            vx=vx+ax*alfa ; vy=vy+ay*alfa ; vz=vz+az*alfa
          else
            cycle
          end if                                                                         !just a bit different (for simplification)
        end do
      end do
      if(vx==0.and.vy==0.and.vz==0)then
        d=0
      else
        d=1/sqrt(vx**2+vy**2+vz**2)
      end if
      cels(celi)%polx=vx*d ; cels(celi)%poly=vy*d ; cels(celi)%polz=vz*d
      !print*,i,"pol",cels(celi)%polx,cels(celi)%poly,cels(celi)%polz
    end if
  end do

end subroutine polarization_single

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

subroutine emt_singleold  !AIXO NO VA MIQUEL >>> Is 4-10-14
integer:: i,j,k,ii,jj,kk,tipi
real*8:: a,b,c,aa,bb,cc,aaa,bbb,ccc,dotp,d,dd,mean,an

  an=-pi/6
  do i=1,nd
    tipi=node(i)%tipus
    if(tipi==2)then
      ii=0 !tipus 1 counter
      jj=0 !tipus 2 counter
      mean=0
      k=node(i)%altre
      a=node(i)%x ; b=node(i)%y ; c=node(i)%z
      aaa=node(k)%x ; bbb=node(k)%y ; ccc=node(k)%z
      do j=1,nneigh(i)
        kk=neigh(i,j)
        !if(node(kk)%tipus==1)then
        !  ii=ii+1
        !  a=a+dneigh(i,j)
        !elseif(node(kk)%tipus==2)then
        !  jj=jj+1
        !  b=b+dneigh(i,j)
        !end if
        if(node(kk)%tipus==2)then
          jj=jj+1
          aa=node(kk)%x ; bb=node(kk)%y ; cc=node(kk)%z
          d=sqrt((aa-a)**2+(bb-b)**2+(cc-c)**2)
          dd=sqrt((aaa-a)**2+(bbb-b)**2+(ccc-c)**2)
          dotp=((aa-a)*(aaa-a)+(bb-b)*(bbb-b)+(cc-c)*(ccc-c))/(dd*d)
          mean=mean+dotp
        end if
      end do
      !a=a/real(ii)
      !b=b/real(jj)
      mean=mean/real(jj) ; print*,i,"TRANS mean",mean

      if(mean<sin(an))then
      !print*,i,"TRANS","a",a,"b",b
      !if(a<b)then !the basal node is getting out of the basal epithelial lining
      node(i)%tipus=3 ; node(k)%tipus=3
      node(i)%altre=0 ; node(k)%altre=0
      cels(node(i)%icel)%ctipus=3
        !gex(i,3)=0  !THIS IS PROVISIONAL!!!!!!!!!!!!!!!!!!!!!!!!!!
        print*,"TRANSITIOOOOOOOOOOON",i,k,"a",a,"b",b
      end if
    end if
  end do

end subroutine emt_singleold

subroutine emt_single  !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)
      s=s+gex(ii,kk)*gen(kk)%wa(nparam_per_node+13)
    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
    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.5d0
        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
      end if
    end do
  end if
end do e

end subroutine

end module single_node

max_stars_repo_name: isalazarciudad/EmbryoMaker

max_issues_repo_name: sebjameswml/EmbryoMaker

max_forks_repo_name: sebjameswml/EmbryoMaker


text:
        IMPLICIT REAL*8(A-H,O-Z)
        DIMENSION A(6),AE(6),AO(6),p(2),po(2),pe(2),ap(4),apo(4),ape(4)
	 dimension xab(36),aogal(6),agal(6)
        DIMENSION AGO(6)
	character*150 arxi, arxo
       character afact*1
        INCLUDE 'const_math.h'
        INCLUDE 'const_ast.h'
        INCLUDE 'const_pop_color.h'
        INCLUDE 'const_pot.h'

c Initializations:
        idum=-15350

        WRITE(*,*)'INPUT FILE WITH GALACTOCENTRIC COORDINATES'
        READ(*,*)ARXI
        WRITE(*,*)ARXI

        WRITE(*,*)'OUTPUT FILE TO WRITE OBSERVED GALACTOCENTRIC COORD.'
        READ(*,*)ARXO
        WRITE(*,*)ARXO

        NCI=23
        OPEN(NCI,FILE=ARXI)

        NCO=4
        OPEN(NCO,FILE=ARXO)

	iff=0
	ies=0
	ia=0
	ifa=0

9	continue

c Reading the galactocentric coordinates of the star
        ip=0
10	continue
	read(nci,*,err=10,end=111)X,Y,Z,VX,VY,VZ,aMv,vi_o
c        write(*,*)X,Y,Z,VX,VY,VZ,aMv,vi_o
        ip=ip+1

c Transform the position of the star from galactocientric to heliocentric
c coordinates
         dist=dsqrt((x+R0)*(x+R0)+y*y+z*z)   !in kpc
         distpc=dist*1000.d0                 !in pc
         xpi=1.d0/dist                       !parallax in mas
         CALL carte_to_equatorial(X,Y,Z,ALPHA,DELTA)
         a(1)=alpha
         a(2)=delta
         a(3)=xpi

C Transforming velocities from galactocentric to galactic helicocentric
c coordinates and to equatorial heliocentric coordinates
         call Carte_to_UVWH(VX,VY,VZ,UH,VH,WH)
         CALL UVWH_to_equatorials(UH,VH,WH,A)

c galactocentric coordinates to galactic coordinates
c         call carte_to_galactic(X,Y,Z,gl,gb)
	 call lb_cal(ALPHA,DELTA,gl,gb)
         call gal_mu(dist,gl,gb,UH,VH,WH,xmub,xmuls,vr)

c---------------------------------------------------------------
c Drimmel extinction law, input itype=2 (equatorial cood) and
c afact='y', using scaling factors
         itype=3
         afact='y'
         call  DrimmelAv(dist,gl/deg,gb/deg,itype,afact,Av,ifg,ifac)
c         write(*,*)Av
c------------------------------------------------------------------
c do not consider stars with Av greater than 40 (fittings of G not considered)
         if(Av.gt.40.d0)goto 10
c  K apparent magnitude of the source
c         xK=aMk+5.d0*log10(distpc)-5.d0+0.112d0*Av
c  V apparent magnitude of the source
         V=aMv+5.d0*log10(distpc)-5.d0 + Av
c Compute the observed (V-I) of the source
         vi=factE*Av+vi_o

c Compute the observed (J-K) of the source
c         xJK=factE2*Av+jk_o

c Assign atmospheric parameters to the source
         ap(1)=Teff
         ap(2)=xlogg
         ap(3)=FeH    !No radial metallicity gradient is assumed here
         ap(4)=Av     !we hare assuming A0=Av

c Introduce Gaia errors
         jflag=-1                          !weighted errors
c       jflag= 1                          !mean errors
c         write(*,*)V,VI,Av
         call Gaia_errors(jflag,V,VI,a,ao,ae,p,po,pe,ap,apo,
     &                      ape,GRVS,idum)
c         write(*,*)p(1)
c	if (iflag.eq.0) then
c	 iff=iff+1
c	go to 10
c	end if
c---------------------------------------------------------
C transform velocities from equatorial heliocentric to
c galactocentric
         CALL equatorials_to_UVWH(ao,Uobs,Vobs,Wobs)
         CALL UVWH_TO_Carte(Uobs,Vobs,Wobs,Vxobs,VYobs,VZobs)

c transform heliocentric to galactocentric positions
         distobs=(1.d0/ao(3))
         call equatorial_to_carte(distobs,ao(1),ao(2),xobs,yobs,zobs)
         call lb_cal(ao(1),ao(2),globs,gbobs)
         distobs=1.d0/ao(3) !in kpc

	call gal_mu(distobs,globs,gbobs,Uobs,Vobs,Wobs,xmubo,xmulso
     &              ,vro)
c Computation of (U,V,W) with errors
        call Gal_to_UVW(Xobs,Yobs,Zobs,VXobs,VYobs,VZobs,Uobs,Vobs,Wobs)

c Computation of (U,V,W) without errors
         call Gal_to_UVW(X,Y,Z,VX,VY,VZ,UU1,VV1,WW1)

         Gmag=p(1)
         vrad=a(6)
         vrob=ao(6)
         xpar=a(3)
         xparobs=ao(3)
         relerr=AE(3)/A(3)
         relerr_mu=AE(5)/A(5)
         relerr_vrad=AE(6)/A(6)

         if (Gmag.le.20.7) then

         WRITE(NCO,100)Av,V,Gmag,GRVS,relerr,X,Y,Z,VX,VY,VZ,
     &              xpar,gl/deg,gb/deg,dist,xmuls*1000.d0,
     &              xmub*1000.d0,vrad,Xobs,Yobs,Zobs,VXobs,VYobs,VZobs
     &              ,xparobs,globs/deg,gbobs/deg,distobs
     &              ,xmulso*1000.d0,xmubo*1000.d0,vrob
         endif
 100     FORMAT(31(F18.10,1X))
         GOTO 10
 111      CLOSE(NCI)
        CLOSE(nco)

        END


text:
c---+----1----+----2----+----3----+----4
      program main
      double precision price_pay,price_org,discout,tax
      discount=0.1d0
      tax=0.05d0
      price_org=980.0d0
      price_pay=price_org*(1.0d0-discount)*(1.0d0+tax)
      write(6,*) 'You should pay',price_pay,'Yen'
      stop
      end

max_stars_repo_name: terasakisatoshi/Fortran

max_issues_repo_name: terasakisatoshi/Fortran

max_forks_repo_name: terasakisatoshi/Fortran