text:
      SUBROUTINE ZGET36( RMAX, LMAX, NINFO, KNT, NIN )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            KNT, LMAX, NIN, NINFO
      DOUBLE PRECISION   RMAX
*     ..
*
*  Purpose
*  =======
*
*  ZGET36 tests ZTREXC, a routine for reordering diagonal entries of a
*  matrix in complex Schur form. Thus, ZLAEXC computes a unitary matrix
*  Q such that
*
*     Q' * T1 * Q  = T2
*
*  and where one of the diagonal blocks of T1 (the one at row IFST) has
*  been moved to position ILST.
*
*  The test code verifies that the residual Q'*T1*Q-T2 is small, that T2
*  is in Schur form, and that the final position of the IFST block is
*  ILST.
*
*  The test matrices are read from a file with logical unit number NIN.
*
*  Arguments
*  ==========
*
*  RMAX    (output) DOUBLE PRECISION
*          Value of the largest test ratio.
*
*  LMAX    (output) INTEGER
*          Example number where largest test ratio achieved.
*
*  NINFO   (output) INTEGER
*          Number of examples where INFO is nonzero.
*
*  KNT     (output) INTEGER
*          Total number of examples tested.
*
*  NIN     (input) INTEGER
*          Input logical unit number.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      COMPLEX*16         CZERO, CONE
      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
     $                   CONE = ( 1.0D+0, 0.0D+0 ) )
      INTEGER            LDT, LWORK
      PARAMETER          ( LDT = 10, LWORK = 2*LDT*LDT )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IFST, ILST, INFO1, INFO2, J, N
      DOUBLE PRECISION   EPS, RES
      COMPLEX*16         CTEMP
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   RESULT( 2 ), RWORK( LDT )
      COMPLEX*16         DIAG( LDT ), Q( LDT, LDT ), T1( LDT, LDT ),
     $                   T2( LDT, LDT ), TMP( LDT, LDT ), WORK( LWORK )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           ZCOPY, ZHST01, ZLACPY, ZLASET, ZTREXC
*     ..
*     .. Executable Statements ..
*
      EPS = DLAMCH( 'P' )
      RMAX = ZERO
      LMAX = 0
      KNT = 0
      NINFO = 0
*
*     Read input data until N=0
*
   10 CONTINUE
      READ( NIN, FMT = * )N, IFST, ILST
      IF( N.EQ.0 )
     $   RETURN
      KNT = KNT + 1
      DO 20 I = 1, N
         READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
   20 CONTINUE
      CALL ZLACPY( 'F', N, N, TMP, LDT, T1, LDT )
      CALL ZLACPY( 'F', N, N, TMP, LDT, T2, LDT )
      RES = ZERO
*
*     Test without accumulating Q
*
      CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDT )
      CALL ZTREXC( 'N', N, T1, LDT, Q, LDT, IFST, ILST, INFO1 )
      DO 40 I = 1, N
         DO 30 J = 1, N
            IF( I.EQ.J .AND. Q( I, J ).NE.CONE )
     $         RES = RES + ONE / EPS
            IF( I.NE.J .AND. Q( I, J ).NE.CZERO )
     $         RES = RES + ONE / EPS
   30    CONTINUE
   40 CONTINUE
*
*     Test with accumulating Q
*
      CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDT )
      CALL ZTREXC( 'V', N, T2, LDT, Q, LDT, IFST, ILST, INFO2 )
*
*     Compare T1 with T2
*
      DO 60 I = 1, N
         DO 50 J = 1, N
            IF( T1( I, J ).NE.T2( I, J ) )
     $         RES = RES + ONE / EPS
   50    CONTINUE
   60 CONTINUE
      IF( INFO1.NE.0 .OR. INFO2.NE.0 )
     $   NINFO = NINFO + 1
      IF( INFO1.NE.INFO2 )
     $   RES = RES + ONE / EPS
*
*     Test for successful reordering of T2
*
      CALL ZCOPY( N, TMP, LDT+1, DIAG, 1 )
      IF( IFST.LT.ILST ) THEN
         DO 70 I = IFST + 1, ILST
            CTEMP = DIAG( I )
            DIAG( I ) = DIAG( I-1 )
            DIAG( I-1 ) = CTEMP
   70    CONTINUE
      ELSE IF( IFST.GT.ILST ) THEN
         DO 80 I = IFST - 1, ILST, -1
            CTEMP = DIAG( I+1 )
            DIAG( I+1 ) = DIAG( I )
            DIAG( I ) = CTEMP
   80    CONTINUE
      END IF
      DO 90 I = 1, N
         IF( T2( I, I ).NE.DIAG( I ) )
     $      RES = RES + ONE / EPS
   90 CONTINUE
*
*     Test for small residual, and orthogonality of Q
*
      CALL ZHST01( N, 1, N, TMP, LDT, T2, LDT, Q, LDT, WORK, LWORK,
     $             RWORK, RESULT )
      RES = RES + RESULT( 1 ) + RESULT( 2 )
*
*     Test for T2 being in Schur form
*
      DO 110 J = 1, N - 1
         DO 100 I = J + 1, N
            IF( T2( I, J ).NE.CZERO )
     $         RES = RES + ONE / EPS
  100    CONTINUE
  110 CONTINUE
      IF( RES.GT.RMAX ) THEN
         RMAX = RES
         LMAX = KNT
      END IF
      GO TO 10
*
*     End of ZGET36
*
      END

max_stars_repo_name: nashmit/lapack-release

max_issues_repo_name: nashmit/lapack-release

max_forks_repo_name: nashmit/lapack-release


text:
      FUNCTION ArcMe(Latitude)
c,<901226.1130>
C     ==============
C
C     Calcul de la longueur de l'Arc de Meridien
C     Latitude en Radians , Resultat en metres
C
      Implicit none
C
      REAL*8 Alpha,ArcMe,Beta,Constante,Delta,Gamma,J2,J4,J6,J8,
     *    Latitude,U,V
C
      DATA Alpha/0.0050761277D0/,Beta/0.0000429451D0/,
     *    Constante/6399936.608D0/,Delta/0.0000000013D0/,
     *    Gamma/0.0000001696D0/
C
      V=DCOS(Latitude)
      U=V*DSIN(Latitude)
      V=V**2
      J2=Latitude+U
      J4=(3.*J2+2.*U*V)/4.
      J6=(5.*J4+2.*U*V**2)/3.
      J8=(7.*J6+4.*U*V**3)/8.
C
      ArcMe=Constante*(Latitude-Alpha*J2+Beta*J4-Gamma*J6+Delta*J8)
C
      END
      FUNCTION ArcMer(Ellips,Latitude1,Latitude2)
c,<901226.1130>
C     ===============
C
C     Calcul de l'arc de meridien entre 2 latitudes donnees
C
C->   Ellips(2) - a0,e2 parametres de l'ellipsoide
C->   ArcMer= arc de Latitude1 a Latitude2
C
      Implicit None
C
      INTEGER No
C
      REAL*8 ArcMer,b(8),Beta0,Beta2,Beta4,Beta6,Beta8,e1,Ellips(2),
     *    Latitude1,Latitude2,p1,p2,p3,p4,s1
C
      s1=1.
      e1=1.-Ellips(2)
      b(1)=Ellips(1)*e1
      DO No=1,4
          s1=s1*Ellips(2)*(2*No+1)/(2*No)
          b(2*No)=Ellips(1)*e1*s1
      ENDDO
C
      Beta0=b(1)+.5*b(2)+3./8.*b(4)+5./16.*b(6)+35./128.*b(8)
      Beta2=.5*(b(2)+b(4))+15./32.*b(6)+28./64.*b(8)
      Beta4=1./8.*b(4)+3./16.*b(6)+28./128.*b(8)
      Beta6=1./32.*b(6)+4./64.*b(8)
      Beta8=1./128.*b(8)
C
      p1=Beta0*(Latitude2-Latitude1)
      p2=.5*Beta2*(DSIN(2.*Latitude2)-DSIN(2.*Latitude1))
      p3=.25*Beta4*(DSIN(4.*Latitude2)-DSIN(4.*Latitude1))
      p4=1./6.*Beta6*(DSIN(6.*Latitude2)-DSIN(6.*Latitude1))
C
      ArcMer=p1-p2+p3-p4
C
      END
      FUNCTION DdDmd(Degres_Decimaux)
c,<901226.1130>
C     ==============
C
C     Conversion Deg. decimaux --> Deg. Min. Decimales (DDDMM.MM)
C
      Implicit None
C
      INTEGER Nb_Degres
C
      REAL*8 DdDmd,Degres_Decimaux
C
      Nb_Degres=Degres_Decimaux
C
      DdDmd=Nb_Degres*100.+(Degres_Decimaux-Nb_Degres)*60.
C
      END
      FUNCTION DdDms(Degres_Decimaux)
c,<901226.1130>
C     ==============
C
C     Conversion des Degres Decimaux en Degres Minutes Secondes sous la
C         forme DDDMM.SEC
C
      Implicit None
C
      INTEGER Minute,Nb_Degres
C
      REAL*8 DdDms,Degres_Decimaux,Reste,Seconde
C
      Nb_Degres=Degres_Decimaux
C
      Reste=(Degres_Decimaux-Nb_Degres)*60.
      Minute=Reste
C
      Seconde=Reste-Minute
C
      DdDms=Nb_Degres*100.+Minute+Seconde/100.*60.
C
      END
      FUNCTION DdRad(Degres_Decimaux)
c,<901226.1130>
C     ==============
C
C     Conversion de Degres_Decimaux en Radian
C
      Implicit none
C
      REAL*8 DdRad,Degres_Decimaux
C
      REAL*8 Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      COMMON /Constante/Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
C     Conv_Deg_Rad=DATAN(1.D0)/45.
C
      DdRad=Conv_Deg_Rad*Degres_Decimaux
C
      END
      FUNCTION DmdDd(Deg_Mn_Decimale)
c,<901226.1130>
C     =============
C
C     Conversion Deg. Min. Decimales (DDDMM.MM) --> Deg. Decimaux
C
      Implicit None
C
      INTEGER Nb_Degres
C
      REAL*8 Deg_Mn_Decimale,DmdDd,Mn_Decimale
C
      Nb_Degres=Deg_Mn_Decimale/100.
      DmdDd=Nb_Degres
C
      Mn_Decimale=Deg_Mn_Decimale-Nb_Degres*100.
      DmdDd=DmdDd+Mn_Decimale/60.
C
      END
      FUNCTION DmdRad(Deg_Mn_Decimale)
c,<901226.1130>
C     ===============
C
C     Conversion Deg. Min. Decimales (DDDMM.MM) --> Radian
C
      Implicit None
C
      REAL*8 Deg_Mn_Decimale,DmdDd,DmdRad
C
      REAL*8 Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      COMMON /Constante/Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
C     Conv_Deg_Rad=DATAN(1.D0)/45.
C
      DmdRad=DmdDd(Deg_Mn_Decimale)*Conv_Deg_Rad
C
      END
      FUNCTION DmsDd(Deg_Mn_Sec)
c,<901226.1130>
C     ==============
C
C     Conversion Deg Min Sec (DDDNb_Minutes.Seconde) --> Deg Decimaux
C
      Implicit None
C
      INTEGER Nb_Minutes,Nb_Degres
C
      REAL*8 Deg_Mn_Sec,DmsDd,Mn_Decimale,Seconde
C
      Nb_Degres=Deg_Mn_Sec/100.
      DmsDd=Nb_Degres
C
      Mn_Decimale=Deg_Mn_Sec-Nb_Degres*100.
      Nb_Minutes=Mn_Decimale
C
      Seconde=Mn_Decimale-Nb_Minutes
C
      Mn_Decimale=Nb_Minutes+Seconde*100./60.
      DmsDd=DmsDd+Mn_Decimale/60.
C
      END
      SUBROUTINE GeLam(Longitude,Latitude,X_Lambert,Y_Lambert,
C     ================
     *    No_Systeme)
c,<901226.1130>
C
C     Conversion Geographique (Radians)  -  Lambert (Metres)
C
C     CALL InitG(No_Systeme) initialise le systeme
C
      Implicit None
C
      INTEGER No_Systeme,No_Systeme_Ref
C
      REAL*8 Constante_Lamb,Delta_Long_Ref,Lat_Croissante,
     *    Lat_Croiss_Ref,Lat_Conforme,Latitude,Latitude_Ref,Longitude,
     *    Longitude_Ref,R,R0,Sinus_Lat_Ref,X_Lambert_Ref,X_Lambert,
     *    Y_Lambert_Ref,Y_Lambert
C
      COMMON /GLAM/Sinus_Lat_Ref,X_Lambert_Ref,Y_Lambert_Ref,R0,
     *    Longitude_Ref,Lat_Croiss_Ref,Delta_Long_Ref,No_Systeme_Ref
C
      IF(No_Systeme.NE.No_Systeme_Ref)CALL InitG(No_Systeme)
C
      Constante_Lamb=(Longitude-Longitude_Ref-Delta_Long_Ref)*
     *    Sinus_Lat_Ref
C
      CALL LACRE(Latitude,Lat_Conforme,Lat_Croissante)
C
      R=R0*DEXP((Lat_Croiss_Ref-Lat_Croissante)*Sinus_Lat_Ref)
C
      X_Lambert=X_Lambert_Ref+R*DSIN(Constante_Lamb)
      Y_Lambert=Y_Lambert_Ref+R0-R*DCOS(Constante_Lamb)
C
      END
      SUBROUTINE geotm(Ellips,Param_Projection,Longitude_Rad,
C     ================
     *    Latitude_Rad,X_Cart,Y_Cart)
c,<901226.1130>
C
C->   geographiques -> transverse Mercator
C->   Ellips(2) - a0,e2   : parametres de l'Ellipsoide
C->   Param_Projection(5) : cml,fi0,rk0 : long,lat et facteur d'echelle
C                                         point central
C->                         x0,y0       : constantes point central
C->   Longitude1,Latitude_Rad --> X_Cart,Y_Cart
C->   source j. pagarete
C
      Implicit None
C->
      REAL*8 a0,ArcMer,c,c2,c4,c6,e2,Ellips(2),Latitude_Rad,Longitude1,
     *    Longitude2,Longitude_Rad,Param_Projection(5),qr,qr2,qr3,rk1,
     *    rk2,rk3,rk4,rk5,rk6,rn,rp,s,sigma,t,t2,t4,t6,w,X_Cart,Y_Cart
C
      a0=Ellips(1)
      e2=Ellips(2)
      Longitude1=Longitude_Rad-Param_Projection(1)
      sigma=ArcMer(Ellips,Param_Projection(2),Latitude_Rad)
      s=DSIN(Latitude_Rad)
      w=DSQRT(1.-e2*s*s)
      rn=a0/w
      qr=(w*w)/(1.-e2)
      qr2=qr*qr
      qr3=qr2*qr
      c=DCOS(Latitude_Rad)
      rp=rn*c
      t=s/c
      Longitude2=Longitude1*Longitude1
      c2=c*c*Longitude2
      c4=c2*c2
      c6=c4*c2
      t2=t*t
      t4=t2*t2
      t6=t4*t2
      rk1=qr-t2
      rk2=rk1+4.*qr2
      rk3=qr*(14.-58.*t2)+40.*t2+t4-9.
      rk4=(8.*qr3*(11.-24.*t2)-28.*qr2*(1.-6.*t2)+qr*(1.-32.*t2)-2.*t2)*
     *    qr+t4
      rk5=61.-479.*t2+179.*t4-t6
      rk6=1385.-3111.*t2+543.*t4-t6
      Y_Cart=sigma+.5*rp*Longitude2*s*(1.+c2*rk2/12.+c4*rk4/360.+c6*rk6/
     *    20160.)
      Y_Cart=Param_Projection(3)*Y_Cart+Param_Projection(5)
      X_Cart=Param_Projection(3)*rp*Longitude1*(1.+c2*rk1/6.+c4*rk3/120.
     *    +c6*rk5/5040.)+Param_Projection(4)
C
      end
      SUBROUTINE Geo_Glr(Longitude_Rad,Latitude_Rad,X_Glr,Y_Glr)
c,
C     ==================
c    *    <901226.1130>
C
C-    geographiques (Radians) -> Gauss Laborde Reunion (Metres)
C
      Implicit None
C
      LOGICAL Init
C
      REAL*8 DmdRad,Ellips(2),Latitude_Rad,Longitude_Rad,
     *    Param_Projection(5),RadDmd,X_Glr,Y_Glr
C
C-    Ellipsoide international a0,e2
C
      DATA Ellips/6378200.d0,0.00669342162d0/
C
C-    Projection Gauss laborde Reunion
C     long - lat du centre de Projection facteur d'echelle + x0,y0
C
      DATA Param_Projection/5532.d0,-2107.d0,1.d0,1.6d5,5.d4/
C
      Param_Projection(1)=DmdRad(Param_Projection(1))
      Param_Projection(2)=DmdRad(Param_Projection(2))
C
      CALL geotm(Ellips,Param_Projection,Longitude_Rad,Latitude_Rad,
     *    X_Glr,Y_Glr)
C
      Param_Projection(1)=RadDmd(Param_Projection(1))
      Param_Projection(2)=RadDmd(Param_Projection(2))
C
      END
      SUBROUTINE Glr_Geo(X_Glr,Y_Glr,Longitude_Rad,Latitude_Rad)
c,
C     ==================
c    *    <901226.1130>
C
C-    Gauss Laborde Reunion (Metres) -> geographiques (Radians)
C
      Implicit None
C
      REAL*8 DmdRad,Ellips(2),Latitude_Rad,Longitude_Rad,
     *    Param_Projection(5),RadDmd,X_Glr,Y_Glr
C
C-    Ellipsoide international a0,e2
C
      DATA Ellips/6378200.d0,0.00669342162d0/
C
C-    Projection Gauss laborde Reunion
C     long - lat du centre de Projection facteur d'echelle + x0,y0
C
      DATA Param_Projection/5532.d0,-2107.d0,1.d0,1.6d5,5.d4/
C
      Param_Projection(1)=DmdRad(Param_Projection(1))
      Param_Projection(2)=DmdRad(Param_Projection(2))
C
      CALL tmgeo(Ellips,Param_Projection,X_Glr,Y_Glr,Longitude_Rad,
     *    Latitude_Rad)
C
      Param_Projection(1)=RadDmd(Param_Projection(1))
      Param_Projection(2)=RadDmd(Param_Projection(2))
C
      END
      FUNCTION GradRad(Grade)
c,<901226.1130>
C     ================
C
C     Conversion de Grades en Radian
C
      Implicit none
C
      REAL*8 Grade,GradRad
C
      REAL*8 Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      COMMON /Constante/Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
C     Conv_Deg_Rad=DATAN(1.D0)/50.
C
      GradRad=Conv_Grad_Rad*Grade
C
      END
      SUBROUTINE GUTM(Longitude,Latitude,X_UTM,Y_UTM,No_Fuseau_UTM)
c,
C     ===============
c    *    <901226.1130>
C
C     Conversion Geographique (Radian)  -  UTM (metres)
C
C     Longitude et Latitude en Radians.
C     No_Fuseau_UTM=(Longitude+186.)/6.
C     Utilise le Sspg ArcMe
C
      Implicit None
C
      INTEGER No_Fuseau_UTM
C
      REAL*8 ArcMe,Constante_UTM,Cosinus,DdRad,ECOXI,Eta,Eprim,GN,
     *    Latitude,Longitude,Longitud_Relativ,X_UTM,XI,Y_UTM
C
      DATA Constante_UTM/6399936.608D0/,Eprim/0.0822688896073373D0/
C
      Longitud_Relativ=Longitude-DdRad(No_Fuseau_UTM*6.D0-183.D0)
      Cosinus=DCOS(Latitude)
C
      XI=Cosinus*DSIN(Longitud_Relativ)
      XI=0.5*DLOG((1.+XI)/(1.-XI))
C
      Eta=DATAN(DSIN(Latitude)/(Cosinus*DCOS(Longitud_Relativ)))-
     *    Latitude
      GN=Constante_UTM/DSQRT(1.+(Eprim*Cosinus)**2)
      ECOXI=(Eprim*Cosinus*XI)**2
C
      X_UTM=500000.+0.9996*GN*XI*(1.+ECOXI/6.)
C
      Y_UTM=0.9996*(GN*Eta*(1.+ECOXI/2.)+ArcMe(Latitude))
C
      END
      SUBROUTINE InitG(No_Systeme)
c,<901226.1130>
C     ================
C
C     Initialise les constantes de conversion Geographiques - Lambert
C         No_Systeme de 1 a 10
C
      Implicit None
C
      INTEGER No_Systeme,No_Systeme_Ref

      REAL*8 DdRad,Delta_Long_Ref,GradRad,Lat_Croiss_Ref,
     *    Lat_Conforme_Ref,Latitude_Ref,Longitude_Ref,Sinus_Lat_Ref,R0,
     *    Param_Projection(6,10),X_Lambert_Ref,Y_Lambert_Ref
C
      COMMON /GLAM/Sinus_Lat_Ref,X_Lambert_Ref,Y_Lambert_Ref,R0,
     *    Longitude_Ref,Lat_Croiss_Ref,Delta_Long_Ref,No_Systeme_Ref
C
      DATA Param_Projection /
C Syst  Latitude_Ref  R0              Longitude_Ref   X    Y   DELT
C                                                   Lambert_Ref
C
C    1                  France   Nord
C
     *  55.D0, 5457616.68D0, 0.D0,     0.6D6,        0.2D6, 2.3372292D0,
C
C    2                  France   Centre
C
     *  52.D0, 5999695.77D0, 0.D0,     0.6D6,        0.2D6, 2.3372292D0,
C
C    3                  France   Sud
C
     *  49.D0, 6591905.08D0, 0.D0,     0.6D6,        0.2D6, 2.3372292D0,
C
C    4                  France   Corse
C
     *46.85D0, 7053300.18D0, 0.D0, 234.538D0, 185861.369D0, 2.3372292D0,

C
C    5                  Nord     Maroc
C
     *  37.D0,9716290.594D0,-6.D0,     0.5D6,        0.3D6,        0.D0,
C
C    6                  Sud      Maroc
C
     *  33.D0,11187308.682D0,-6.D0,    0.5D6,        0.3D6,        0.D0,
C
C    7                  Nord     Algerie
C
     *  40.D0,8785951.541D0, 3.D0,      5.D5,         3.D5,        0.D0,
C
C    8                  Sud      Algerie
C
     *  37.D0,9716290.594D0, 3.D0,      5.D5,         3.D5,        0.D0,
C
C    9                  Nord     Tunisie
C
     *  40.D0,8785951.541D0,11.D0,      5.D5,         3.D5,        0.D0,
C
C   10                  Sud      Tunisie
C
     *  37.D0,9716290.594D0,11.D0,      5.D5,         3.D5,        0.D0/
C
      No_Systeme_Ref=No_Systeme
      X_Lambert_Ref=Param_Projection(4,No_Systeme)
      Y_Lambert_Ref=Param_Projection(5,No_Systeme)
      R0=Param_Projection(2,No_Systeme)
C
      Latitude_Ref=GradRad(Param_Projection(1,No_Systeme))
      Longitude_Ref=GradRad(Param_Projection(3,No_Systeme))
      Delta_Long_Ref=DdRad(Param_Projection(6,No_Systeme))
      Sinus_Lat_Ref=DSIN(Latitude_Ref)
C
      CALL LACRE(Latitude_Ref,Lat_Conforme_Ref,Lat_Croiss_Ref)
C
      END
      SUBROUTINE LACRE(Latitude,Lat_Conforme,Lat_Croissante)
c,
C     ================
c     *    <901226.1130>
C
C     Sspg de calcul de Latitude Conforme et de Latitude Croissante
C
C     Calcul de la Latitude conforme Lat_Conforme
C     Calcul de la Latitude Croissante Ellipsoidique Lat_Croissante
C     Donnees et Resultats en Radians
C
      Implicit None
C
      REAL*8 Constante,Lat_Conforme,Lat_Croissante,Latitude,Teta,Valeur
C
      DATA Constante/0.0824832568D0/
C
      Valeur=0.5*Latitude
      Teta=DTAN(Valeur)
      Valeur=2.*Constante*Teta/(1.+Teta**2)
      Valeur=0.5*Constante*DLOG((1.+Valeur)/(1.-Valeur))
      Lat_Croissante=DLOG((1.+Teta)/(1.-Teta))-Valeur
      Valeur=DTANH(0.5*Valeur)
      Lat_Conforme=2.*DATAN((Teta-Valeur)/(1.-Valeur*Teta))
C
      END
      SUBROUTINE LamGe(X_Lambert,Y_Lambert,Longitude,Latitude,
C     ================
     *    No_Systeme)
c,<901226.1130>
C
C     Conversion Lambert (Metres)  -  Geographiques (Radians)
C     CALL InitG(No_Systeme)  Initialise le systeme
C     Ellipsoide de CLARKE IGN
C
      Implicit None
C
      INTEGER No_Systeme,No_Systeme_Ref
C
      REAL*8 Delta_Long_Ref,Latitude,Lat_Crois_Ref,Lat_Croissante,
     *    Longitude,Longitude_Ref,R,R0,RLAC,Sinus_Lat_Ref,X_Lambert,
     *    X_Lambert_Ref,Y_Lambert,Y_Lambert_Ref
C
      COMMON /GLAM/Sinus_Lat_Ref,X_Lambert_Ref,Y_Lambert_Ref,R0,
     *    Longitude_Ref,Lat_Crois_Ref,Delta_Long_Ref,No_Systeme_Ref
C
      IF(No_Systeme.NE.No_Systeme_Ref)CALL InitG(No_Systeme)
C
      Longitude=X_Lambert-X_Lambert_Ref
      Latitude=R0-(Y_Lambert-Y_Lambert_Ref)
      R=Latitude*DSQRT(1.+(Longitude*Longitude)/(Latitude*Latitude))
      Longitude=Longitude_Ref+DATAN(Longitude/Latitude)/Sinus_Lat_Ref
      Lat_Croissante=Lat_Crois_Ref+DLOG(R0/R)/Sinus_Lat_Ref
      Lat_Croissante=DTANH(0.5D0*Lat_Croissante)
C
      Latitude=RLAC(Lat_Croissante)
      Longitude=Longitude+Delta_Long_Ref
C
      END
      FUNCTION RadDd(Radian)
c,<901226.1130>
C     ==============
C
C     Conversion de radians en degres decimaux
C
      Implicit None
C
      REAL*8 RadDd,Radian
C
      REAL*8 Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      COMMON /Constante/Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
C     Conv_Deg_Rad=DATAN(1.D0)/45.
C
      RadDd=Radian/Conv_Deg_Rad
C
      END
      FUNCTION RadDmd(Radian)
c,<901226.1130>
C     ===============
C
C     Conversion de Radians en Degres Minutes decimales
C
      Implicit None
C
      REAL*8 DdDmd,RadDmd,Radian
C
      REAL*8 Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      COMMON /Constante/Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
C     Conv_Deg_Rad=datan(1.d0)/45.
C
      RadDmd=DdDmd(Radian/Conv_Deg_Rad)
C
      END
      FUNCTION RLAC(Tangente_Lat_xxx)
c,<901226.1130>
C     =============
C
C     Fonction inverse du sspg LACRE
C     Calcul retour de latitude croissante ou de latitude conforme
C
C         Tangente_Lat_xxx=TAN(Lat_Conforme/2.)=TANH(Lat_Croissante/2.)
C
C     Resultat : Latitude en Radians
C
      Implicit None
C
      INTEGER Nb
C
      REAL*8 Constante,Q,RLAC,Tangente_Lat_xxx,Teta
C
      DATA Constante/0.0824832568D0/
C
      Teta=Tangente_Lat_xxx
      DO Nb=1,3
          Q=2.*Constante*Teta/(1.+Teta**2)
          Q=DTANH(0.25*Constante*DLOG((1.+Q)/(1.-Q)))
          Teta=(Tangente_Lat_xxx+Q)/(1.+Q*Tangente_Lat_xxx)
      ENDDO
C
      RLAC=2.*DATAN(Teta)
C
      END
      SUBROUTINE tmgeo(Ellips,Param_Projection,X_Cart,Y_Cart,
C     ================
     *    Longitude_Rad,Latitude_Rad)
c,<901226.1130>
C
C->   transverse Mercator -> geogrphiques
C->   Ellips - a0,e2 - parametres de l'Ellipsoide
C->   Param_Projection - cm0,fi0,rk0,x0,y0
C-    coordonnes du point central et facteur d'echelle
C->   X_Cart,Y_Cart --> Longitude_Rad,Latitude_Rad
C->          metres --> radians
C
      Implicit None
C
      INTEGER No
C     INTEGER LU_1
C
      REAL*8 a0,c,dLatitude,dLongitude,dist,dist2,dx,dy,e2,Ellips(2),
     *    Latitude_Rad,Longitude_Rad,prec2,Param_Projection(5),qr,rn,rp,
     *    s,w,X_Cart,xe,Y_Cart,ye
C
C     COMMON /LUecran/LU_1
C
      a0=Ellips(1)
      e2=Ellips(2)
      prec2=0.0001
      No=0
C
C->   on obtient une approximation de Longitude_Rad et Latitude_Rad par
C         une formule standard
C
      CALL tmgeo1(Ellips,Param_Projection,X_Cart,Y_Cart,Longitude_Rad,
     *    Latitude_Rad)
   10 CALL geotm(Ellips,Param_Projection,Longitude_Rad,Latitude_Rad,xe,
     *    ye)
      dx=X_Cart-xe
      dy=Y_Cart-ye
C     IF(No.EQ.0)WRITE(LU_1,*)'dx dy --> ',dx,dy
      s=DSIN(Latitude_Rad)
      w=DSQRT(1.-e2*s*s)
      rn=a0/w
      qr=(w*w)/(1.-e2)
      c=DCOS(Latitude_Rad)
      rp=rn*c
      dLatitude=dy*qr/rn
      dLongitude=dx/rp
      Latitude_Rad=Latitude_Rad+dLatitude
      Longitude_Rad=Longitude_Rad+dLongitude
      No=No+1
      IF(No.LE.20)then
          dist2=dx*dx+dy*dy
          if(dist2.GT.prec2)GO TO 10
C     ELSE
C         WRITE(LU_1,*)' tmgeo : pas de convergence '
      ENDIF
C
      END
      SUBROUTINE tmgeo1(Ellips,Param_Projection,X_Cart,Y_Cart,
C     =================
     *    Longitude_Rad,Latitude_Rad)
c,<901226.1130>
C
C->   transverse Mercator --> geographiques
C->   formule directe approchee (utilisee dans tmgeo)
C->   l'erreur n'excede pas 2 metres dans un fuseau de +- 10 deg
C->   adapte de utmgr pour des Ellipsoides variables
C->   a0,e2 - parametres de l'Ellipsoide
C->   cm0,fi0,rk0 - long,lat point central et facteur d'echelle
C->   X_Cart,Y_Cart --> Longitude_Rad,Latitude_Rad
C
      Implicit None
C
      real*8 a0,a1,a2,ArcMer,cofip,cophi,e2,ecofi,ecxn2,Ellips(2),eprim,
     *    eprim2,eta,expxi,gn,gphi,Latitude_Rad,Longitude_Rad,phi,phipr,
     *    Param_Projection(5),siphi,vm2,x,X_Cart,xi,y,Y_Cart
c
      EQUIVALENCE (cofip,ecxn2),(gphi,expxi),(vm2,ecofi)
C
C-    gain de place ?
C
      a0=Ellips(1)
      e2=Ellips(2)
      eprim2=e2/(1.-e2)
      eprim=DSQRT(eprim2)
      a1=a0/DSQRT(1.-e2)
      a2=a0*DSQRT(1.+eprim2)
      y=Y_Cart-Param_Projection(5)
      y=y/Param_Projection(3)
      x=(X_Cart-Param_Projection(4))/Param_Projection(3)
      phi=y/a2+Param_Projection(2)
      siphi=DSIN(phi)
      cophi=DCOS(phi)
      ecofi=eprim*cophi
      gn=a1/DSQRT(1.+ecofi**2)
      y=y-ArcMer(Ellips,Param_Projection(2),phi)
      xi=x/gn
      ecxn2=(ecofi*xi)**2
      xi=xi*(1.-ecxn2/6.)
      eta=(y/gn)*(1.-ecxn2/2.)
      phipr=phi+eta
      cofip=DCOS(phipr)
      expxi=DEXP(xi)
      Longitude_Rad=DATAN((expxi**2-1.)/(2.*expxi*cofip))
      gphi=DATAN(DCOS(Longitude_Rad)*DSIN(phipr)/cofip)
      vm2=1.+ecofi**2-1.5*eprim*siphi*ecofi*(gphi-phi)
      Longitude_Rad=Longitude_Rad+Param_Projection(1)
      Latitude_Rad=phi+vm2*(gphi-phi)
C
      END
      SUBROUTINE TRANSF(X_Init,Y_Init,X_Final,Y_Final,No_Systeme_Init,
C     =================
     *    No_Systeme_Final,No_Ref_Init,No_Ref_Final)
c,<1.1>
C
C     Conversion de coordonnees
C              X_Init,Y_Init      en       X_Final,Y_Final
C         No_Systeme_Init  -  No_Systeme_Final
C              1      2       3       4       5       6
C             UTM    DMS     DD      DMD     LAM     GLR
C     No_Ref_Init,No_Ref_Final Parametres supplementaires (No Lambert
C         ou Fuseau UTM)
C     Dans le cas des Geographiques : Deplacement du meridien origine
C              1      2       3       4       5
C             GRNW   PARIS   ROME    .....    ......
C
      IMPLICIT None
C
      INTEGER LU_1,No_Ref_Final,No_Ref_Init,No_Reference,
     *    No_Systeme_Final,No_Systeme_Init
C
      REAL*8 DdDmd,DdDms,DdRad,DmdDd,DmsDd,Longitude_Ref(5),RadDd,
     *    X_Init,X_Final,X_Inter,Y_Init,Y_Final,Y_Inter
C
      REAL*8 Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      COMMON /Constante/Conv_Deg_Rad,Pi,Conv_Grad_Rad
C
      DATA No_Reference/0/,Longitude_Ref/0.D0,2.3372292D0,12.452222D0,
     *    0.D0,0.D0/
C
C     Initialisation unique via le Common Constante sinon la mettre dans
C         TRANSG ou TRANSF suivant le programme :
C             GEODIS et PARNU (GEO) : TRANSG
C             COORD  et HYPO        : TRANSF
C
      Pi=4.*DATAN(1.D0)
      Conv_Deg_Rad=Pi/180.
      Conv_Grad_Rad=Pi/200.
C
      IF(No_Systeme_Init.EQ.No_Systeme_Final)GO TO 90
   10 GO TO (20,30,40,50,60,70),No_Systeme_Init
C
      WRITE(*,'("Transf : Systeme initial indefini")')
      RETURN
C
   20 CALL UTMGR(X_Init*1000.,Y_Init*1000.,X_Inter,Y_Inter,No_Ref_Init)
C
      GO TO 80
C
   30 X_Final=DmsDd(X_Init)+Longitude_Ref(No_Ref_Init)
      Y_Final=DmsDd(Y_Init)
      GO TO 100
C
   40 X_Final=X_Init+Longitude_Ref(No_Ref_Init)
      Y_Final=Y_Init
      GO TO 100
C
   50 X_Final=DmdDd(X_Init)+Longitude_Ref(No_Ref_Init)
      Y_Final=DmdDd(Y_Init)
      GO TO 100
C
   60 CALL LamGe(1000.*X_Init,1000.*Y_Init,X_Inter,Y_Inter,No_Ref_Init)
      GO TO 80
C
   70 CALL Glr_Geo(1000.*X_Init,1000.*Y_Init,X_Inter,Y_Inter)
C
   80 X_Final=RadDd(X_Inter)
      Y_Final=RadDd(Y_Inter)
      GO TO 100
C
   90 IF(No_Ref_Init.NE.No_Ref_Final)GO TO 10
      X_Final=X_Init
      Y_Final=Y_Init
      RETURN
C
  100 GO TO (110,120,130,140,150,160),No_Systeme_Final
C
      WRITE(*,'("Transf :  Systeme final indefini")')
      RETURN
C
  110 IF(No_Ref_Final.NE.0)No_Reference=No_Ref_Final
      IF(No_Reference.EQ.0)No_Reference=(X_Final+186.)/6.
      No_Ref_Final=No_Reference
C
      CALL GUTM(DdRad(X_Final),DdRad(Y_Final),X_Inter,Y_Inter,
     *    No_Reference)
      GO TO 170
C
  120 X_Final=DdDms(X_Final-Longitude_Ref(No_Ref_Final))
      Y_Final=DdDms(Y_Final)
      RETURN
C
  130 X_Final=X_Final-Longitude_Ref(No_Ref_Final)
      RETURN
C
  140 X_Final=DdDmd(X_Final-Longitude_Ref(No_Ref_Final))
      Y_Final=DdDmd(Y_Final)
      RETURN
C
  150 CALL GeLam(DdRad(X_Final),DdRad(Y_Final),X_Inter,Y_Inter,
     *    No_Ref_Final)
      GO TO 170
C
  160 CALL Geo_Glr(DdRad(X_Final),DdRad(Y_Final),X_Inter,Y_Inter)
C
  170 X_Final=X_Inter/1000.
      Y_Final=Y_Inter/1000.
C
      END
      SUBROUTINE UTMGR(X_UTM,Y_UTM,Longitude,Latitude,No_Fuseau_Ref)
c,
C     ================
c    *    <901226.1130>
C
C     Conversion Coordonnees UTM en Geographiques
C
C     Donnees en Metres , Resultats en Radians
C     Origine des longitudes : le meridien international
C     Definition normale de la projection mais Y_UTM est negatif dans
C         l'hemisphere sud , prendre ABS(Y_UTM) pour avoir le NORTHING
C         habituel
C
C     Longitude_Ref est la longitude du meridien central du fuseau
C
      Implicit None
C
      INTEGER No_Fuseau_Ref
C
      REAL*8 ArcMe,C,COFIP,CoPhi,DdRad,ECOFI,ECXN2,Eprim,Eta,EXPXI,GN,
     *    GPhi,Latitude,Longitude,Longitude_Ref,Phi,PhiPr,SiPhi,VM2,X,
     *    X_UTM,XI,Y,Y_UTM
C
      EQUIVALENCE(COFIP,ECXN2),(GPhi,EXPXI),(VM2,ECOFI)
C
      DATA C/6399936.608D0/,Eprim/0.0822688896073373D0/
C
      Longitude_Ref=DdRad(No_Fuseau_Ref*6.D0-183.)
      Y=DABS(Y_UTM)
C
C     IF(Y_UTM.LT.0.0)Y=-Y_UTM-10000000.
C
      Y=Y/0.9996
      Phi=Y/6366197.724
      SIPhi=DSIN(Phi)
      COPhi=DCOS(Phi)
      ECOFI=Eprim*COPhi
      GN=C/DSQRT(1.+ECOFI**2)
      X=(X_UTM-500000.)/0.9996
      Y=Y-ArcMe(Phi)
      XI=X/GN
      ECXN2=(ECOFI*XI)**2
      XI=XI*(1.-ECXN2/6.)
      Eta=(Y/GN)*(1.-ECXN2/2.)
      PhiPR=Phi+Eta
      COFIP=DCOS(PhiPR)
      EXPXI=DEXP(XI)
      Longitude=DATAN((EXPXI**2-1.)/(2.*EXPXI*COFIP))
      GPhi=DATAN(DCOS(Longitude)*DSIN(PhiPR)/COFIP)
      VM2=1.+ECOFI**2-1.5*Eprim*SIPhi*ECOFI*(GPhi-Phi)
      Longitude=Longitude+Longitude_Ref
      Latitude=Phi+VM2*(GPhi-Phi)
      IF(Y_UTM.LT.0.)Latitude=-Latitude
C
      END

ext: f

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
  MODULE mo_salsa_driver
  USE classSection, ONLY : Section
  USE util, ONLY : getMassIndex !!! IS it good to import this here??? The function is anyway handy here too.
  USE mo_salsa_types, ONLY : aero, cloud, precp, ice, allSALSA, rateDiag, iaero, icloud, iprecp, iice
  USE mo_submctl
  USE classFieldArray, ONLY : FieldArray
  USE mo_structured_datatypes
  IMPLICIT NONE

   !---------------------------------------------------------------
   !
   ! MO_SALSA_DRIVER:
   ! Contains the primary SALSA input/output variables as well as
   ! Subroutines used to call the main SALSA routine.
   !
   ! Juha Tonttila, FMI, 2014
   !
   !---------------------------------------------------------------

  SAVE

  ! Derived type for taking copies of the "old" values
  TYPE old
     REAL, ALLOCATABLE :: d(:)
  END TYPE old
  INTERFACE old
     PROCEDURE old_cnstr
  END INTERFACE old

   ! JT: Variables from SALSA
   ! --------------------------------------------
   ! grid points for SALSA
   INTEGER, PARAMETER :: kproma = 1
   INTEGER, PARAMETER :: kbdim = 1
   INTEGER, PARAMETER :: klev = 1
   INTEGER, PARAMETER :: krow = 1

   REAL, PARAMETER    :: init_rh(kbdim,klev) = 0.3

   ! -- Local gas compound tracers [# m-3]
   REAL :: zgso4(kbdim,klev),   &
           zghno3(kbdim,klev),  &
           zgnh3(kbdim,klev),   &
           zgocnv(kbdim,klev),  &
           zgocsv(kbdim,klev)

   ! To help coupling and backing up for tendencies
   TYPE(old)          :: npart(4), mpart(4)  ! To store the old values
   TYPE(FloatArray1d) :: ntend(4), mtend(4)  ! Arranged pointers to tendencies

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

CONTAINS

   !
   !----------------------------------------------------
   ! RUN_SALSA
   ! Performs necessary unit and dimension conversion between
   ! the host model and SALSA module, and calls the main SALSA
   ! routine
   !
   ! Partially adobted form the original SALSA boxmodel version.
   !
   ! Now takes masses in as kg/kg from LES!! Converted to m3/m3 for SALSA

   ! Juha Tonttila, FMI, 2014
   ! Jaakko Ahola, FMI, 2016
   !
  SUBROUTINE run_SALSA(Diag, Prog, nzp, nxp, nyp, ns, wp,    &
                       pa_nactd, pa_vactd, tstep, time,      &
                       level, initialize                     )

      USE mo_salsa, ONLY : salsa
      USE mo_salsa_properties, ONLY  : equilibration
      IMPLICIT NONE

      TYPE(FieldArray), INTENT(inout) :: Diag, Prog

      INTEGER, INTENT(in) :: nzp,nxp,nyp,ns                       ! dimensions: x,y,z,number of chemical species
      REAL, INTENT(in)    :: tstep                                ! Model timestep length
      REAL, INTENT(in)    :: time
      LOGICAL, INTENT(in) :: initialize
      REAL, INTENT(in)    :: wp(nzp,nxp,nyp)
      REAL, INTENT(out)   :: pa_vactd(nzp,nxp,nyp,ns*ncld) ! mass concentrations of newly activated droplets for calculating the
                                                           ! actual tendency due to new droplet formation.
      REAL, INTENT(out)   :: pa_nactd(nzp,nxp,nyp,ncld)   ! Same for number concentration
      INTEGER, INTENT(in) :: level                         ! thermodynamical level

      ! Named pointers to get SALSA variables
      TYPE(FloatArray3d), POINTER :: press => NULL(), tk => NULL(), rv => NULL(), &
                                     rt => NULL(), rs => NULL(), rsi => NULL(),   &
                                     pdn => NULL()
      TYPE(FloatArray4d), POINTER :: naerop => NULL(), ncloudp => NULL(),   &
                                     nprecpp => NULL(), nicep => NULL(),    &
                                     naerot => NULL(), ncloudt => NULL(),   &
                                     nprecpt => NULL(), nicet => NULL(),    &

                                     maerop => NULL(), mcloudp => NULL(),   &
                                     mprecpp => NULL(), micep => NULL(),    &
                                     maerot => NULL(), mcloudt => NULL(),   &
                                     mprecpt => NULL(), micet => NULL(),    &

                                     indefp => NULL(),                      &

                                     gaerop => NULL(), gaerot => NULL()

      REAL :: in_p(kbdim,klev), in_t(kbdim,klev), in_rv(kbdim,klev), in_rs(kbdim,klev),&
              in_w(kbdim,klev), in_rsi(kbdim,klev)
      REAL :: rv_old(kbdim,klev)

      ! For process rate diagnostics
      TYPE(FloatArray3d), POINTER :: autoconversion => NULL(), accretion => NULL(),             &
                                     ACcoll => NULL(), APcoll => NULL(), AIcoll => NULL(),      &
                                     activation => NULL(), Icehom => NULL(), Icedep => NULL(),  &
                                     Iceimm => NULL(), Conda => NULL(), Condc => NULL(),        &
                                     Condp => NULL(), Condi => NULL()

      TYPE(Section) :: actd(kbdim,klev,ncld) ! Activated droplets - for interfacing with SALSA

      INTEGER :: jj,ii,kk,str,end,nc,nb,nbloc,ndry,nwet,iwa,irim,icat

      INTEGER :: bin_numbers(4), bin_starts(4) ! Number of bins and the starting indices in the
                                               ! allSALSA array according to the phase indexing in classSection.f90
      INTEGER :: catnbins

      ndry = spec%getNSpec(type="dry")
      nwet = spec%getNSpec(type="wet")  ! excludes rimed ice
      iwa = spec%getIndex("H2O")    ! water/unrimed ice
      irim = spec%getIndex("rime")  ! rimed ice; returns 0 if level<5 and rime not used

      ! The order here should follow the phase indentifier indexing!!
      bin_numbers(1:4) = [nbins,ncld,nprc,nice]
      bin_starts(1:4) = [iaero,icloud,iprecp,iice]

      ! Get the pointers to variable arrays
      CALL Diag%getData(1,press,name="press")
      CALL Diag%getData(1,tk,name="temp")
      CALL Diag%getData(1,rs,name="rsl")
      IF (level == 5) &
           CALL Diag%getData(1,rsi,name="rsi")
      CALL Diag%getData(1,pdn,name="dn")
      CALL Prog%getData(1,rv,name="rp")
      CALL Prog%getData(2,rt,name="rp")
      ! Note that maybe there could be a "getSeveral" method or something to make this less verbose?
      CALL Prog%getData(1,naerop,name="naero")
      CALL Prog%getData(2,naerot,name="naero")
      CALL Prog%getData(1,ncloudp,name="ncloud")
      CALL Prog%getData(2,ncloudt,name="ncloud")
      CALL Prog%getData(1,nprecpp,name="nprecp")
      CALL Prog%getData(2,nprecpt,name="nprecp")
      IF (level == 5) THEN
           CALL Prog%getData(1,nicep,name="nice")
           CALL Prog%getData(2,nicet,name="nice")
      END IF

      CALL Prog%getData(1,maerop,name="maero")
      CALL Prog%getData(2,maerot,name="maero")
      CALL Prog%getData(1,mcloudp,name="mcloud")
      CALL Prog%getData(2,mcloudt,name="mcloud")
      CALL Prog%getData(1,mprecpp,name="mprecp")
      CALL Prog%getData(2,mprecpt,name="mprecp")
      IF (level == 5) THEN
         CALL Prog%getData(1,micep,name="mice")
         CALL Prog%getData(2,micet,name="mice")
         IF (ice_theta_dist) CALL Prog%getData(1,indefp,name="indef")
      END IF

      CALL Prog%getData(1,gaerop,name="gaero")
      CALL Prog%getData(2,gaerot,name="gaero")

      ! Get pointers to process rate output streams
      CALL Diag%getData(1,autoconversion,name="autoconversion")
      CALL Diag%getData(1,accretion,name="accretion")
      CALL Diag%getData(1,ACcoll,name="ACcoll")
      CALL Diag%getData(1,APcoll,name="APcoll")
      CALL Diag%getData(1,activation,name="activation")
      IF (level == 5) THEN
         CALL Diag%getData(1,AIcoll,name="AIcoll")
         CALL Diag%getData(1,Icehom,name="icehom")
         CALL Diag%getData(1,Icedep,name="icedep")
         CALL Diag%getData(1,Iceimm,name="iceimm")
         CALL Diag%getData(1,condi,name="condi")
      END IF
      CALL Diag%getData(1,conda,name="conda")
      CALL Diag%getData(1,condc,name="condc")
      CALL Diag%getData(1,condp,name="condp")

      ! Initialize the arrays for managing "old" values and tendencies
      IF (initialize) &
           CALL initialize_arrays(nwet,bin_numbers)

      in_p(:,:) = 0.; in_t(:,:) = 0.; in_rs(:,:) = 0.; in_rsi(:,:) = 0.; in_w(:,:) = 0.
      in_rv(:,:) = 0.; rv_old(:,:) = 0.

      ! Set the SALSA runtime config
      CALL set_salsa_runtime(time)

      ! Convert input concentrations for SALSA into #/m3 or m3/m3 instead of kg/kg (multiplied by pdn/divided by substance density)
      DO jj = 3, nyp-2
         DO ii = 3, nxp-2
            DO kk = nzp-1, 2, -1

               ! Set inputs
               in_p(1,1) = press%d(kk,ii,jj)
               in_t(1,1) = tk%d(kk,ii,jj)
               in_rs(1,1) = rs%d(kk,ii,jj)
               IF (level == 5) &
                    in_rsi(1,1) = rsi%d(kk,ii,jj)
               in_w(1,1) = wp(kk,ii,jj)

               ! For initialization and spinup, limit the RH with the parameter rhlim (assign in namelist.salsa)
               IF ( lsfreeRH%state ) THEN
                  in_rv(1,1) = rv%d(kk,ii,jj)
               ELSE
                  in_rv(1,1) = MIN(rv%d(kk,ii,jj), rs%d(kk,ii,jj)*rhlim)
               END IF
               rv_old(1,1) = in_rv(1,1)

               ! Update volume concentrations
               ! ---------------------------------------------------------------------------------------------------

               ! Take a copies for the old values. Note that the order of the categories in these arrays is CRUCIAL.
               ! It follows the phase indexing in classSection. Use the derived type "old"
               ! Also make use of these when preparing SALSA inputs below. this is a bit verbose, but what are you gonna do...
               npart(1)%d = naerop%d(kk,ii,jj,:)
               npart(2)%d = ncloudp%d(kk,ii,jj,:)
               npart(3)%d = nprecpp%d(kk,ii,jj,:)
               IF (level == 5) &
                    npart(4)%d = nicep%d(kk,ii,jj,:)
               mpart(1)%d = maerop%d(kk,ii,jj,:)
               mpart(2)%d = mcloudp%d(kk,ii,jj,:)
               mpart(3)%d = mprecpp%d(kk,ii,jj,:)
               IF (level == 5) &
                    mpart(4)%d = micep%d(kk,ii,jj,:)

               ! Make arranged pointers (not copies!) to tendency arrays in the same manner
               ntend(1)%d => naerot%d(kk,ii,jj,:)
               ntend(2)%d => ncloudt%d(kk,ii,jj,:)
               ntend(3)%d => nprecpt%d(kk,ii,jj,:)
               IF (level == 5) &
                    ntend(4)%d => nicet%d(kk,ii,jj,:)
               mtend(1)%d => maerot%d(kk,ii,jj,:)
               mtend(2)%d => mcloudt%d(kk,ii,jj,:)
               mtend(3)%d => mprecpt%d(kk,ii,jj,:)
               IF (level == 5) &
                    mtend(4)%d => micet%d(kk,ii,jj,:)

               ! Update SALSA input arrays
               DO nb = 1,ntotal
                  icat = allSALSA(1,1,nb)%phase   ! Phase indentifier, this should correspond to index in bin_starts and bin_numbers, mpart and npart, ntend and mtend
                  IF (icat == 4 .AND. level < 5) EXIT  ! skip ice if level < 5; since ice gets the last index in classSection, we can use EXIT

                  nbloc = nb - bin_starts(icat) + 1
                  catnbins = bin_numbers(icat)

                  ! Update volume concentrations. Note, rime handled separately below
                  IF (icat < 4) THEN
                     DO nc = 1,nwet
                        str = getMassIndex(catnbins,nbloc,nc)
                        allSALSA(1,1,nb)%volc(nc) = mpart(icat)%d(str)*pdn%d(kk,ii,jj)/spec%rholiq(nc)
                     END DO

                     IF (level == 5 .AND. ice_theta_dist) THEN
                        allSALSA(1,1,nb)%INdef = indefp%d(kk,ii,jj,nb)
                     END IF

                  ELSE IF (icat == 4) THEN
                     DO nc = 1,nwet
                        str = getMassIndex(catnbins,nbloc,nc)
                        allSALSA(1,1,nb)%volc(nc) = mpart(icat)%d(str)*pdn%d(kk,ii,jj)/spec%rhoice(nc)
                     END DO
                     str = getMassIndex(catnbins,nbloc,irim)
                     allSALSA(1,1,nb)%volc(irim) = mpart(icat)%d(str)*pdn%d(kk,ii,jj)/spec%rhori

                  END IF

                  ! Update number concentrations
                  allSALSA(1,1,nb)%numc = npart(icat)%d(nbloc)*pdn%d(kk,ii,jj)

                  ! Update particle size and densities
                  !CALL allSALSA(1,1,nb)%updateDiameter(.TRUE.,type="all")
                  CALL allSALSA(1,1,nb)%updateRhomean()

                  IF (allSALSA(1,1,nb)%numc > allSALSA(1,1,nb)%nlim) THEN
                     allSALSA(1,1,nb)%core = SUM(allSALSA(1,1,nb)%volc(1:ndry))/allSALSA(1,1,nb)%numc
                  ELSE
                     allSALSA(1,1,nb)%core = pi6*(allSALSA(1,1,nb)%dmid)**3
                  END IF

               END DO

               ! If this is an initialization call, calculate the equilibrium particle
               If (initialize) CALL equilibration(kproma,kbdim,klev,   &
                                                  init_rh,in_t,.TRUE.)

               ! Convert to #/m3
               zgso4(1,1) = gaerop%d(kk,ii,jj,1)*pdn%d(kk,ii,jj)
               zghno3(1,1) = gaerop%d(kk,ii,jj,2)*pdn%d(kk,ii,jj)
               zgnh3(1,1) = gaerop%d(kk,ii,jj,3)*pdn%d(kk,ii,jj)
               zgocnv(1,1) = gaerop%d(kk,ii,jj,4)*pdn%d(kk,ii,jj)
               zgocsv(1,1) = gaerop%d(kk,ii,jj,5)*pdn%d(kk,ii,jj)

               ! Reset process rate diagnostics
               CALL rateDiag%Reset()

               ! ***************************************!
               !                Run SALSA               !
               ! ***************************************!
               CALL salsa(kproma, kbdim,  klev,   krow,     &
                          in_p,   in_rv,  in_rs,  in_rsi,   &
                          in_t,   tstep,  zgso4,  zgocnv,   &
                          zgocsv, zghno3, zgnh3,  actd,     &
                          in_w,   level                     )

               ! Update tendency arrays
               DO nb = 1,ntotal
                  icat = allSALSA(1,1,nb)%phase   ! Phase indentifier, this should correspond to index in npart and mpart
                  IF (icat == 4 .AND. level < 5) EXIT  ! skip ice if level < 5; since ice gets the last index in classSection, we can use EXIT

                  nbloc = nb - bin_starts(icat) + 1
                  catnbins = bin_numbers(icat)

                  ! Update mass tendencies
                  IF (icat < 4) THEN
                     DO nc = 1,nwet
                        str = getMassIndex(catnbins,nbloc,nc)
                        mtend(icat)%d(str) = mtend(icat)%d(str) + &
                             ( (allSALSA(1,1,nb)%volc(nc)*spec%rholiq(nc)/pdn%d(kk,ii,jj)) -   &
                               mpart(icat)%d(str) ) / tstep
                     END DO

                     IF (level == 5 .AND. ice_theta_dist) THEN
                        indefp%d(kk,ii,jj,nb) = allSALSA(1,1,nb)%INdef
                     END IF

                  ELSE IF (icat == 4) THEN
                     DO nc = 1,nwet
                        str = getMassIndex(catnbins,nbloc,nc)
                        mtend(icat)%d(str) = mtend(icat)%d(str) + &
                             ( (allSALSA(1,1,nb)%volc(nc)*spec%rhoice(nc)/pdn%d(kk,ii,jj)) -   &
                               mpart(icat)%d(str) ) / tstep
                     END DO
                     str = getMassIndex(catnbins,nbloc,irim)
                     mtend(icat)%d(str) = mtend(icat)%d(str) + &
                          ( (allSALSA(1,1,nb)%volc(irim)*spec%rhori/pdn%d(kk,ii,jj)) -     &
                            mpart(icat)%d(str) ) / tstep
                  END IF

                  ! Update number tendencies
                  ntend(icat)%d(nbloc) = ntend(icat)%d(nbloc) + &
                       ( (allSALSA(1,1,nb)%numc/pdn%d(kk,ii,jj)) - npart(icat)%d(nbloc) )/tstep

               END DO

               ! Activated droplets
               pa_nactd(kk,ii,jj,1:ncld) = actd(1,1,1:ncld)%numc/pdn%d(kk,ii,jj)

               DO nc = 1,iwa
                  ! Activated droplets (in case of cloud base activation)
                  str = getMassIndex(ncld,1,nc)
                  end = getMassIndex(ncld,ncld,nc)
                  pa_vactd(kk,ii,jj,str:end) = actd(1,1,1:ncld)%volc(nc)*spec%rholiq(nc)/pdn%d(kk,ii,jj)
               END DO

               IF (lscndgas) THEN
                  gaerot%d(kk,ii,jj,1) = gaerot%d(kk,ii,jj,1) + &
                                          ( (zgso4(1,1)/pdn%d(kk,ii,jj)) - gaerop%d(kk,ii,jj,1) )/tstep

                  gaerot%d(kk,ii,jj,2) = gaerot%d(kk,ii,jj,2) + &
                                          ( (zghno3(1,1)/pdn%d(kk,ii,jj)) - gaerop%d(kk,ii,jj,2) )/tstep

                  gaerot%d(kk,ii,jj,3) = gaerot%d(kk,ii,jj,3) + &
                                          ( (zgnh3(1,1)/pdn%d(kk,ii,jj)) - gaerop%d(kk,ii,jj,3) )/tstep

                  gaerot%d(kk,ii,jj,4) = gaerot%d(kk,ii,jj,4) + &
                                          ( (zgocnv(1,1)/pdn%d(kk,ii,jj)) - gaerop%d(kk,ii,jj,4) )/tstep

                  gaerot%d(kk,ii,jj,5) = gaerot%d(kk,ii,jj,5) + &
                                          ( (zgocsv(1,1)/pdn%d(kk,ii,jj)) - gaerop%d(kk,ii,jj,5) )/tstep
               END IF

               ! Tendency of water vapour mixing ratio
               rt%d(kk,ii,jj) = rt%d(kk,ii,jj) + &
                  ( in_rv(1,1) - rv_old(1,1) )/tstep

               ! Store the process rate diagnostics for output
               autoconversion%d(kk,ii,jj) = rateDiag%Autoconversion%volc(iwa)*spec%rhowa/pdn%d(kk,ii,jj)
               accretion%d(kk,ii,jj) = rateDiag%Accretion%volc(iwa)*spec%rhowa/pdn%d(kk,ii,jj)
               ACcoll%d(kk,ii,jj) = SUM(rateDiag%ACcoll%volc(1:ndry)*spec%rholiq(1:ndry))/pdn%d(kk,ii,jj)
               APcoll%d(kk,ii,jj) = SUM(rateDiag%APcoll%volc(1:ndry)*spec%rholiq(1:ndry))/pdn%d(kk,ii,jj)
               activation%d(kk,ii,jj) = rateDiag%Activation%numc/pdn%d(kk,ii,jj)
               conda%d(kk,ii,jj) = rateDiag%Cond_a%volc(iwa)/pdn%d(kk,ii,jj)
               condc%d(kk,ii,jj) = rateDiag%Cond_c%volc(iwa)/pdn%d(kk,ii,jj)
               condp%d(kk,ii,jj) = rateDiag%Cond_p%volc(iwa)/pdn%d(kk,ii,jj)
               IF (level == 5) THEN
                  AIcoll%d(kk,ii,jj) = SUM(rateDiag%AIcoll%volc(1:ndry)*spec%rhoice(1:ndry))/pdn%d(kk,ii,jj)
                  Icehom%d(kk,ii,jj) = rateDiag%Ice_hom%numc/pdn%d(kk,ii,jj)
                  Icedep%d(kk,ii,jj) = rateDiag%Ice_dep%numc/pdn%d(kk,ii,jj)
                  Iceimm%d(kk,ii,jj) = rateDiag%Ice_imm%numc/pdn%d(kk,ii,jj)
                  Condi%d(kk,ii,jj) = ( rateDiag%Cond_i%volc(iwa)*spec%rhoic    +  &
                                        rateDiag%Cond_i%volc(irim)*spec%rhori ) /  &
                                        pdn%d(kk,ii,jj)
               END IF

            END DO !kk
         END DO ! ii
      END DO ! jj

      ! Nullify pointers and clean up
      DO icat = 1,4
         ntend(icat)%d => NULL()
         mtend(icat)%d => NULL()
      END DO

      press => NULL(); tk => NULL(); rv => NULL()
      rt => NULL(); rs => NULL(); rsi => NULL()
      pdn => NULL()
      naerop => NULL(); ncloudp => NULL()
      nprecpp => NULL(); nicep => NULL()
      naerot => NULL(); ncloudt => NULL()
      nprecpt => NULL(); nicet => NULL()

      maerop => NULL(); mcloudp => NULL()
      mprecpp => NULL(); micep => NULL()
      maerot => NULL(); mcloudt => NULL()
      mprecpt => NULL(); micet => NULL()

      gaerop => NULL(); gaerot => NULL()

      indefp => NULL()

   END SUBROUTINE run_SALSA

   !
   !---------------------------------------------------------------
   ! SET_SALSA_RUNTIME
   ! Set the master process %state:s based on the values of %switch and %delay
   !
   ! Juha Tonttila, FMI, 2014
   !
   SUBROUTINE set_SALSA_runtime(time)
     USE mo_submctl, ONLY : Nmaster, lsmaster, lsfreeRH
     IMPLICIT NONE
     REAL, INTENT(in) :: time
     INTEGER :: i

     DO i = 1,Nmaster
        IF( lsmaster(i)%switch .AND. time > lsmaster(i)%delay ) lsmaster(i)%state = .TRUE.
     END DO

     ! Some other switches
     IF ( lsfreeRH%switch .AND. time > lsfreeRH%delay ) lsfreeRH%state = .TRUE.
     IF ( lsFreeTheta%switch .AND. time > lsFreeTheta%delay) lsFreeTheta%state = .TRUE.

   END SUBROUTINE set_SALSA_runtime

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

   SUBROUTINE initialize_arrays(nwet,bins)
     INTEGER, INTENT(in) :: nwet
     INTEGER, INTENT(in) :: bins(4)
     INTEGER :: icat
     DO icat = 1,4
        ntend(icat) = FloatArray1d()
        mtend(icat) = FloatArray1d()
        npart(icat) = old(bins(icat))
        IF (icat < 4) THEN
           mpart(icat) = old(bins(icat)*nwet)
        ELSE
           mpart(icat) = old(bins(icat)*(nwet+1))
        END IF
     END DO
   END SUBROUTINE initialize_arrays

   FUNCTION old_cnstr(n)
     TYPE(old) :: old_cnstr
     INTEGER, INTENT(in) :: n
     ALLOCATE(old_cnstr%d(n))
     old_cnstr%d(:) = 0.
   END FUNCTION old_cnstr

END MODULE mo_salsa_driver

max_stars_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_issues_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_forks_repo_name: UCLALES-SALSA/UCLALES-SALSA


text:
!   FLIPS (Fortran Linear Inverse Problem Solver)
!
!   Copyright 2005-2009 University of Oulu, Finland. All rights reserved.
!   Written by  Mikko Orispaa <mikko.orispaa@oulu.fi>
!
!   Redistribution and use in source and binary forms, with or without modification,
!   are permitted provided that the following conditions are met:
!
!      1. Redistributions of source code must retain the above copyright notice, this
!         list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright notice, this list
!         of conditions and the following disclaimer in the documentation and/or other materials
!         provided with the distribution.
!
!   THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF OULU ``AS IS'' AND ANY EXPRESS OR IMPLIED
!   WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
!   FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF OULU
!   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
!   OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
!   ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
!   ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!   The views and conclusions contained in the software and documentation are those of
!   the authors and should not be interpreted as representing official policies, either
!   expressed or implied, of University of Oulu.

!*****************************************************************
! rind - gives the storage index of i,j-element of matrix R (band version)
!*****************************************************************

function rind(i,j,n,bw) result(res)
  implicit none

  integer, intent(in) :: i,j,n
  integer, intent(in), optional :: bw
  integer :: res

  if (present(bw)) then
     ! If matrix band width is given
     if (i <= n - bw) then
        res = yind(i,j-i+1,bw)
     else
        res = bw*(n-bw) + frind(i-(n-bw),j-(n-bw),bw)
     end if
  else
     ! If the matrix band width is omitted, assume that the matrix is full
     res = j + (i-1)*(2*n - i)/2
  end if

end function rind

max_stars_repo_name: morispaa/flips

max_issues_repo_name: morispaa/flips

max_forks_repo_name: morispaa/flips