text:
	SUBROUTINE IM_LABL ()
C************************************************************************
C* IM_LABL								*
C*									*
C* This subroutine creates the proper label for the color bar based	*
C*	 upon image type.						*
C*									*
C* IM_LABL ()								*
C*									*
C* Input parameters:							*
C*	NONE								*
C* Output parameters:							*
C*	NONE								*
C **									*
C * Log:								*
C * T. Piper/SAIC	07/06	Created					*
C * S. Chiswell/Unidata	07/06	Added check for supplied calibrations	*
C ***********************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'IMGDEF.CMN'
	INCLUDE		'imcmn.cmn'
C
	REAL		PR_TMCK, PR_TMKC
	CHARACTER	dirnam*256, basnam*(MXFLSZ)
C------------------------------------------------------------------------
C
	IF ( imcalbar .eq. 1 ) THEN
C
C*	    We have already created the levels
C
	    imndlv = min ( 256, (immxpx-immnpx+1) )
	    RETURN
	ENDIF
C
C*  For sat images, create the data level labels
C
	IF ( imradf .eq. 0 )  THEN
	    CALL FL_PATH(cmfile, dirnam, basnam, ier)
	    IF ( basnam(1:4) .eq. "CAPE" )  THEN
		imndlv = 64
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "J/kg"
		cmblev(1)  = "0"
		cmblev(5)  = "1000"
		cmblev(12) = "2000"
		cmblev(18) = "3000"
		cmblev(24) = "4000"
		cmblev(30) = "5000"
		cmblev(36) = "6000"
		cmlutf = "gdpicape.tbl"
	    ELSEIF ( basnam(1:4) .eq. "CINH" )  THEN
		imndlv = 64
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "J/kg"
		cmblev(1) = "0"
		cmblev(2) = "30"
		cmblev(4) = "60"
		cmblev(6) = "90"
		cmblev(10) = "120"
		cmlutf = "gdpicinh.tbl"
	    ELSEIF ( basnam(1:2) .eq. "LI" )  THEN
		imndlv = 64
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "Deg C"
		cmblev(25) = "+5"
		cmblev(31) = "0"
		cmblev(38) = "-5"
		cmblev(44) = "-10"
		cmlutf = "gdpili.tbl"
	    ELSEIF ( basnam(1:2) .eq. "PW" )  THEN
		imndlv = 64
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "Inches"
		cmblev(44) = "0.0"
		cmblev(47) = "0.5"
		cmblev(51) = "1.0"
		cmblev(53) = "1.5"
		cmblev(57) = "2.0"
		cmblev(60) = "2.5"
		cmlutf = "micropw.tbl"
	    ELSEIF ( basnam(1:4) .eq. "RAIN" )  THEN
		imndlv = 64
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "In/hour"
		cmblev(27) = "0.00"
		cmblev(32) = "0.25"
		cmblev(37) = "0.50"
		cmblev(41) = "0.75"
		cmblev(44) = "1.00"
		cmblev(47) = "1.25"
		cmblev(50) = "1.50"
		cmblev(53) = "1.75"
		cmblev(56) = "2.00"
		cmlutf = "microrain.tbl"
	    ELSEIF ( basnam(1:3) .eq. "TPW" )  THEN
		imndlv = 64
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "Inches"
		cmblev(1) = "0.0"
		cmblev(10) = "0.5"
		cmblev(19) = "1.0"
		cmblev(29) = "1.5"
		cmblev(38) = "2.0"
		cmblev(49) = "2.5"
		cmlutf = "gdpipw.tbl"
	    ELSEIF ( basnam(1:4) .eq. "WIND" )  THEN
		imndlv = 32
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "Knots"
		cmblev(14) = "0"
		cmblev(19) = "10"
		cmblev(21) = "15"
		cmblev(23) = "19"
		cmblev(25) = "28"
		cmblev(26) = "36"
		cmlutf = "ssmiwind.tbl"
C
C*  For IR images calibrated to brightness, create Deg C levels
C
	    ELSEIF ( cmcalb .eq. "BRIT" .and.
     +		   ( imtype .eq. 8 .or. imtype .eq. 128) )  THEN
		imndlv = min ( 256, (immxpx-immnpx+1) )
		DO ii = 1, imndlv
		    cmblev(ii) = ' '
		ENDDO
		cmbunt = "Deg C"
		CALL IM_BTOT ( 1, immnpx, tmpk, ier)
		tmpc = PR_TMKC(tmpk)
		CALL ST_RLCH ( tmpc, 1, cmblev(1), ier)
		ibrit = (imndlv-1) + immnpx
		CALL IM_BTOT ( 1, ibrit, tmpk, ier )
		tmpc = PR_TMKC(tmpk)
		CALL ST_RLCH ( tmpc, 1, cmblev(imndlv), ier)
		tmpc = -100.0
		DO WHILE ( tmpc .lt. 51.0 )
		    tmpk = PR_TMCK(tmpc)
		    CALL IM_TTOB ( 1, tmpk, brit, ier )
		    ibrit = nint ( brit ) - immnpx
		    IF ( ibrit .gt. 0 )  THEN
			CALL ST_INCH ( INT (tmpc), cmblev(ibrit), ier )
		    ENDIF
		    tmpc = tmpc + 10.0
		ENDDO
C
C*  For anything else, create labels in native units
C
	    ELSE
		cmbunt = cmcalb
		imndlv = min ( 256, (immxpx-immnpx+1) )
		ratio = (immxpx - immnpx)/255.0
		CALL ST_INCH ( immnpx, cmblev(1), ier )
		DO ii = 2, imndlv-1
		    IF ( mod(ii-1,16) .eq. 0 )  THEN
			level = nint((ii-1)*ratio) + immnpx
			CALL ST_INCH ( level, cmblev(ii), ier )
		    ELSE
			cmblev(ii) = ' '
		    ENDIF
		ENDDO
		level = nint( (imndlv - 1) * ratio) + immnpx
		CALL ST_INCH ( level, cmblev(imndlv), ier)
	    ENDIF
	ENDIF
C
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine zgcbd(simul,n,binf,bsup,x,f,g,zero,napmax,itmax,indgc
     &  ,ibloc,nfac,imp,io,epsx,epsf,epsg,dir,df0,diag,x2,
     &izs,rzs,dzs,y,s,z,ys,zs,nt,index,wk1,wk2,alg,ialg,nomf)
c     Copyright INRIA

c
      implicit double precision (a-h,o-z)
      real rzs(*)
      double precision dzs(*)
      dimension x2(n),dir(n),epsx(n)
      dimension binf(n),bsup(n),x(n),g(n),diag(n),ibloc(n),izs(*)
      dimension y(nt,n),s(nt,n),z(nt,n),ys(nt),zs(nt)
      dimension wk1(n),wk2(n),alg(15)
      character*6 nomf
      integer index(nt),ialg(15)
      external simul
c
      if(imp.ge.4) then
      write(io,10000)
10000 format (' dans gcbd. algorithme utilise: ')
      if(ialg(1).eq.1)write(io,10001)
10001 format ('        emploi correction de powell ')
      if(ialg(2).eq.1)write(io,10002)
10002 format ('  mise a jour de diag par la methode bfgs')
      if(ialg(3).eq.1)write(io,10003)
10003 format ('  mise a echelle de diag par methode de shanno-phua')
      if(ialg(3).eq.2)write(io,10004)
10004 format ('  mise a echelle de diag seulement a la 2e iter')
      if(ialg(4).eq.1)write(io,10005)
10005 format ('      memorisation pour choix iteration ')
      if(ialg(5).eq.1)write(io,10006)
10006 format ('      memorisation par variable')
      if(ialg(6).eq.1)write(io,10007)
10007 format ('      relachememt de variables a toutes les iteration')
      if(ialg(6).eq.2)write(io,10008)
10008 format ('      relachement de vars si decroissance g_norme')
      if(ialg(6).eq.10)write(io,10009)
10009 format ('      relachement de vars si dec f % iter_init du cycle')
      if(ialg(6).eq.11)write(io,10010)
10010 format ('      relachement de vars si dec f % dec du cycle')
      if(ialg(7).eq.1)write(io,10011)
10011 format ('      choix de vars a relacher par bertsekas modifiee')
      if(ialg(8).eq.1)write(io,10012)
10012 format ('      choix de dir descente par methode de gradient')
      if(ialg(8).eq.2)write(io,10013)
10013 format ('      choix de dir descente par methode qn')
      if(ialg(8).eq.3)write(io,10014)
10014 format ('      choix de dir descente par qn sans memoire.nt depl')
      if(ialg(8).eq.4)write(io,10015)
10015 format ('      choix de dir descente par qn -mem,redem,sans acc.')
      if(ialg(8).eq.5)write(io,10016)
10016 format ('     choix de dir descente par qn -mem,redem,avec acc.')
      if(ialg(9).eq.2)write(io,10017)
10017 format ('      redem si relachement de vars')
      if(ialg(9).eq.10)write(io,10018)
10018 format ('      redem si dec f % dec iter_init du cycle')
      if(ialg(9).eq.11)write(io,10019)
10019 format ('      redem si dec f % dec totale du cycle.')
      if(ialg(9).eq.12)write(io,10020)alg(9)
10020 format ('    redem si diminution du gradient des var libres d un',
     & 'facteur',d11.4)
      endif
c
c     section 1  initialisations
c     irl nombre de rech lin 'lentes'
c     nred nombre de redemarrage de la direction de descente
c     icycl nombre de cycles de minimisation
c
      epsgcp=1.0d-5
      indsim=4
      indrl=1
      irl=0
      irl=0
      nred=1
      icycl=1
      nap=0
c
      iresul=1
      call proj(n,binf,bsup,x)
      indsim=4
      call simul(indsim,n,x,f,g,izs,rzs,dzs)
      nap=nap+1
      if(indsim.gt.0)go to 99
      indgc=-1
      if(indsim.eq.0)indgc=0
      if(imp.gt.0)write(io,123)indgc
      go to 900
99    continue
      ceps0=20.0d+0
      eps0=0.0d+0
      do 100 i=1,n
100   eps0=eps0+epsx(i)
      eps0=ceps0*eps0/n
c
c     calcul de zng
      znog0=rednor(n,binf,bsup,x,epsx,g)
      zng=znog0
      zngrit=znog0
      zngred=znog0
c
      do 130 i=1,n
130   ibloc(i)=0
      izag=3
      izag1=izag
      nap=0
      iter=0
      scal=1.0d+0
      nfac=n
      np=0
      lb=1
      nb=2
      if(ialg(8).eq.3) nb=1
      do 140 i=1,nt
140   index(i)=i
      tetaq=alg(9)
      condm=alg(2)
      param=alg(1)
      indgc1=indgc
c     si indgc=0 on init diag a k*ident puis scal a it=2
c
      if(indgc.eq.1.or.indgc.ge.100)go to 150
      if(indgc.eq.2)go to 180
      indgc=-13
      if(imp.gt.0) write(io,123) indgc
      go to 900
c
150   continue
c     on initialise diag par approximation quadratique
c     df0 decroissance prevue . si mod quad df0=((dh)-1g,g)/2
c     et on cherche dh diag de la forme cst/(dx)**2
c     donc cst=som((g(i)*(dx))**2))/(2*df0)
      sy=0.0d+0
      do 160 i=1,n
160    sy=sy+(g(i)*epsx(i))**2
      sy=sy/(2.0d+0*df0)
      do 170 i=1,n
170   diag(i)=(sy + zero)/(epsx(i)**2 + zero)
180   continue
c
c
c     bouclage
200   iter=iter +1
      indgc=1
      if(iter.gt.itmax)then
         indgc=5
         go to 900
      endif
201   continue
      if(imp.ge.2)write(io,1210)iter,f
1210  format(/' dans gcbd  iter=',i3,'  f=',d15.7)
      if(iter.eq.1)then
         irit=1
         goto 301
      endif
c
      call majysa(n,nt,np,y,s,ys,lb,g,x,wk2,wk1,index,ialg,nb)
      inp=index(np)
c
c
c     correction powell sur y si (y,s) trop petit
      if(ialg(1).ne.1) go to 290
      param1=1.-param
      bss=0.0d+0
      do 260 i=1,n
260   bss=bss + diag(i)*s(inp,i)**2
      bss2=param*bss
      if(ys(inp).gt.bss2)go to 290
      if(imp.gt.2)write(io,1270)ys(inp)
1270  format(' gcbd. emploi correction powell (y,s)=',d11.4)
      teta=param1*bss/(bss-ys(inp))
      teta1=1.0d+0-teta
      do 274 i=1,n
274   y(inp,i)=teta*y(inp,i)+teta1*diag(i)*s(inp,i)
      ys(inp)=bss2
c     verif correction powell (facultatif; faire go to 300)
      ys1=ddot(n,s(inp,1),1,y(inp,1),1)
      ys1=abs(bss2-ys1)/bss2
      if(imp.gt.2)write(io,1280)ys1
1280  format(' erreur relative correction powell =',d11.4)
c
c mise a jour de diag
290   continue
      if(ialg(2).eq.1)
     &  call bfgsd(diag,n,nt,np,y,s,ys,condm,param,zero,index)
c
      if(ialg(3).eq.1.or.(ialg(3).eq.2.and.iter.eq.2))
     &  call shanph(diag,n,nt,np,y,s,ys,scal,index,io,imp)
c
      call majz(n,np,nt,y,s,z,ys,zs,diag,index)
c
c     section 3 determination des variables libres et bloquees
300   continue
c     -----decision de relachement a l'iteration courante
c          relachement si irit=1 (sinon irit=0)
      irit=0
      if(ialg(6).eq.1) irit=1
      if(ialg(6).eq.2.and.znglib.le.alg(6)*zngrit)irit=1
      if(ialg(6).eq.10.and.diff.le.dfrit1*alg(6))irit=1
      if(ialg(6).eq.11.and.diff.le.difrit*alg(6))irit=1
      if(irit.eq.1) nred=nred+1
c    ----choix des variables a relacher
      imp1=imp
301   if(ialg(7).eq.1)call relvar(ind,n,x,binf,bsup,x2,g,diag,
     &  imp,io,ibloc,izag,iter,nfac,irit)
c
c
c     section 4 expression de dir
      if (np.eq.0) then
         do 400 i=1,n
           dir(i)=-g(i)/diag(i)
400      continue
      else
         do 410 i=1,n
            dir(i)=-scal*g(i)
410      continue
         call gcp(n,index,ibloc,np,nt,y,s,z,ys,zs,diag,g,dir,wk1,
     &            wk2,epsgcp)
      endif
c
c     section 5  redemarrage
c
      if(ialg(8).eq.4.or.ialg(8).eq.5) then
         ired=0
         if(ialg(9).eq.2.and.ind.eq.1) ired=1
         if(ialg(9).eq.10.and.diff.lt.dfred1*tetaq) ired=1
         if(ialg(9).eq.11.and.diff.lt.difred*tetaq)  ired=1
         if(ialg(9).eq.12.and.znglib.le.tetaq*zngred) ired=1
         if(ired.eq.1) then
            icycl=icycl+1
            np=0
            lb=1
            if(imp.gt.2)write(io,1000)icycl
1000        format ('   redemarrage. icycl=',i5)
         endif
      endif
c
c     section 6 annulation de d(i) , i dans ib
      if(ialg(6).eq.1)go to 640
      do 630 i=1,n
630      if(ibloc(i).gt.0) dir(i)=0.0d+0
640   continue
c
c     recherche lineaire
c     conservation de x et g dans wk1 et wk2
      call dcopy(n,x,1,wk1,1)
      call dcopy(n,g,1,wk2,1)
c     calcul de la derivee dans la direction dir
      ifp=0
      fn=f
      znog0=zng
702   dfp=0.0d+0
      do 710 i=1,n
      epsxi=epsx(i)
      xi=x(i)
      diri=dir(i)
      if(xi-binf(i).le.epsxi.and.diri.lt.0.0d+0)dir(i)=0.0d+0
710   if(bsup(i)-xi.le.epsxi.and.diri.gt.0.0d+0)dir(i)=0.0d+0
      dfp=ddot(n,g,1,dir,1)
      if(-dfp.gt.0)go to 715
      if(ifp.eq.1) then
         indgc=6
         go to 900
      endif
c     restauration dir
      if(imp.ge.3)write(io,1712)dfp,zero
1712  format(' gcbd : restauration dir ; fp,zero',2d11.4)
      do 712 i=1,n
712   dir(i)=-scal*g(i)
      ifp=1
      go to 702
715   continue
c     pas initial suivant idee fletcher
      t=-2.0d+0*diff/dfp
      if(iter.eq.1)t=-2.0d+0*df0/dfp
      tmax=1.0d+10
      t=min(t,tmax)
      t=max(t,1.0d+10*zero)
      napm=15
      napm1=nap + napm
      if(napm1.gt.napmax) napm1=napmax
      napav=nap
      amd=0.70d+0
      amf=0.10d+0
c
      call rlbd(indrl,n,simul,x,binf,bsup,f,dfp,t,tmax,dir,g,tproj,
     &  amd,amf,imp,io,zero,nap,napm1,x2,izs,rzs,dzs)
      if(imp.gt.2)write(io,750)indrl,t,f
750   format(' retour mlibd indrl=',i6,' pas= ',d11.4,' f= ',d11.4)
      if(nap-napav.ge.5) irl=irl+1
      if(indrl.ge.10)then
         indsim=4
         nap=nap + 1
         call simul(indsim,n,x,f,g,izs,rzs,dzs)
         if(indsim.le.0)then
            indgc=-3
            if(indsim.eq.0)indgc=0
            if(imp.gt.0)write(io,123)indgc
123         format(' gcbd : retour avec indgc=',i8)
            go to 900
         endif
      endif
      if(indrl.le.0)then
         indgc=10
         if(indrl.eq.0)indgc=0
         if(indrl.eq.-3)indgc=13
         if(indrl.eq.-4)indgc=12
         if(indrl.le.-1000)indgc=11
         if(imp.gt.0)write(io,123)indgc
         go to 900
      endif
      if(imp.ge.5) then
         do 760 i=1,n
760      if(imp.gt.2)write(io,777)i,x(i),g(i),dir(i)
777      format(' i=',i2,' xgd ',3f11.4)
      endif
c
      if(nap.lt.napmax)go to 758
      if(imp.gt.0)  write(io,755)
755   format(' gcbd max appels simul')
      indgc=4
      go to 900
758   continue
c
c     section 8 test de convergence
      do 805 i=1,n
      if(abs(x(i)-wk1(i)).gt.epsx(i))go to 806
805   continue
      if(imp.gt.0) write(io,1805)
1805  format(' gcbd. retour apres convergence sur x')
      indgc=3
      go to 900
c     calcul grad residuel,norme l2
806   continue
      difg=rednor(n,binf,bsup,x,epsx,g)
      diff=fn-f
      if(imp.ge.2)write(io,860)epsg,difg,epsf,diff,nap
860   format(' gcbd. epsg,difg=',2d11.4,'  epsf,diff=',2d11.4
     &,'  nap=',i3)
c
      if(diff.le.epsf) then
         indgc=2
         go to 900
      endif
      if(difg.le.epsg) then
         indgc=1
         go to 900
      endif
c
c     -----mise a jour de difrit,dfrit1,difred,dfred1
      if(irit.eq.1) then
         difrit=diff
         dfrit1=diff
       else
         difrit=difrit+diff
      endif
      if(ired.eq.1) then
         difred=diff
         dfred1=diff
       else
         difred=difred + diff
      endif
c
      znglib=0.0d+0
      do 884 i=1,n
         if(ibloc(i).gt.0)go to 884
         aa=g(i)
         if(x(i)-binf(i).le.epsx(i)) aa=min(0.0d+0,aa)
         if(bsup(i)-x(i).le.epsx(i)) aa=max(0.0d+0,aa)
         znglib=znglib+aa**2
884   continue
      znglib=sqrt(znglib)
      if(ired.eq.1)zngred=znglib
      if(irit.eq.1)zngrit=znglib
      go to 200
c
c      fin des calculs
900   if(indrl.eq.0)indgc=0
      if(indgc.eq.1.and.indrl.le.0)  indgc=indrl
      if(imp.gt.0)write(io,123)indgc
      if(imp.ge.1.and.indrl.le.zero)write(io,1910) indrl
1910  format(' arret impose par la recherche lineaire. cf notice rlbd',
     & /,' indicateur de rlbd=',i6)
      if(imp.ge.1)write(io,950)f,difg,nap,iter,indgc
950   format(' f,norme grad,nap,iter,indgc=',2e11.4,3i6)
c
c     autres impressions finales
      if(indgc1.lt.100) return
      zrl=0.
      if(iter.gt.0) zrl=dble(nap)/dble(iter)
2000  format('     nom    n       f        norm2g   nf   iter  rl/it ',
     &          ' irl   cpu  cycl   red')
      write(io,2001)nomf,f,difg,nap,iter,zrl,irl
2001  format(1x,a6,2e11.4,2i5,f6.2,i5)
      end

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    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.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``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 PSBLAS GROUP OR ITS 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.
!
!
!
!  Sorting routines
!  References:
!  D. Knuth
!  The Art of Computer Programming, vol. 3
!  Addison-Wesley
!
!  Aho, Hopcroft, Ullman
!  Data Structures and Algorithms
!  Addison-Wesley
!
module psb_z_qsort_mod
  use psb_const_mod

  interface psb_qsort
    subroutine psb_zqsort(x,ix,dir,flag)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), optional, intent(in)    :: dir, flag
      integer(psb_ipk_), optional, intent(inout) :: ix(:)
    end subroutine psb_zqsort
  end interface psb_qsort

  interface
    subroutine psi_zlqsrx_up(n,x,ix)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zlqsrx_up
    subroutine psi_zlqsrx_dw(n,x,ix)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zlqsrx_dw
    subroutine psi_zlqsr_up(n,x)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zlqsr_up
    subroutine psi_zlqsr_dw(n,x)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zlqsr_dw
    subroutine psi_zalqsrx_up(n,x,ix)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zalqsrx_up
    subroutine psi_zalqsrx_dw(n,x,ix)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zalqsrx_dw
    subroutine psi_zalqsr_up(n,x)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zalqsr_up
    subroutine psi_zalqsr_dw(n,x)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zalqsr_dw
    subroutine psi_zaqsrx_up(n,x,ix)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zaqsrx_up
    subroutine psi_zaqsrx_dw(n,x,ix)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(inout) :: ix(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zaqsrx_dw
    subroutine psi_zaqsr_up(n,x)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zaqsr_up
    subroutine psi_zaqsr_dw(n,x)
      import
      complex(psb_dpk_), intent(inout)  :: x(:)
      integer(psb_ipk_), intent(in)   :: n
    end subroutine psi_zaqsr_dw
  end interface

end module psb_z_qsort_mod


text:
  CALL MergeMatrixLocalBlocks(this, local)

  !! Merge Columns
  CALL TransposeMatrix(local, localT)
  CALL ReduceAndComposeMatrix(localT, merged_columns, &
       & this%process_grid%column_comm)

  !! Merge Rows
  CALL TransposeMatrix(merged_columns, merged_columnsT)
  CALL ReduceAndComposeMatrix(merged_columnsT, gathered, &
       & this%process_grid%row_comm)

  !! Remove the excess rows and columns that come from the logical size.
  CALL ConstructEmptyMatrix(local_mat, this%actual_matrix_dimension, &
       & this%actual_matrix_dimension)
  local_mat%outer_index = gathered%outer_index(:this%actual_matrix_dimension+1)
  ALLOCATE(local_mat%inner_index(SIZE(gathered%inner_index)))
  local_mat%inner_index = gathered%inner_index
  ALLOCATE(local_mat%values(SIZE(gathered%values)))
  local_mat%values = gathered%values

max_stars_repo_name: Kokookster/NTPoly

max_issues_repo_name: Kokookster/NTPoly

max_forks_repo_name: Kokookster/NTPoly


text:
C    @(#)sen_jac2.f	20.2 6/27/97
      subroutine sen_jac2 (kt,ksw)

C      This subroutine computes Jacobian elements for bus
C      constraints. KSW determines the type.

C        0 - Normal Jacobian row for node
C        1 - Jacobian for Q-contraint only, regardless of bus type

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/alpha2.inc'
      include 'ipfinc/amtrx.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/beta2.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/dc.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/factor.inc'
      include 'ipfinc/gamma.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/optim1.inc'
      include 'ipfinc/pctvr2.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/tran.inc'
      include 'ipfinc/xtran.inc'
C
      real ikr,iki,imr,imi
      integer senchk

      lp=0
      kl=0
      kp = 0
      isw=0
      i1 = iflag(kt)
      i2 = iflag(kt+1) - 1
      if (ksw .eq. 0 .and. kspare(20) .eq. 1) then

C        Check for slack bus (conditional upon AI_CONTROL option
C        on > SENSITIVITY < not enabled).

         do 90 i = i1,i2
         if (jflag(1,i) .eq. 3) then
           isw = jflag(2,i)
           if (senchk(kt) .eq. 0) then
              go to 100
           else

C             Area slack constraint only

              call nrpqv (kt,pk,dp,qk,dq,vk)
              lp=1
              kl = lp
              korder(0)=1
              kolum(0)=0
              korder(1)=0
              kolum(kl) = kt
              rowh(kl) = 0.0
              rown(kl) = 0.0
              rowj(kl) = 0.0
              rowl(kl) = 1.0
              max = kt
              mend=2
              go to 294
           endif
         endif
   90 continue
      endif
  100 ek=e(kt)
      fk=f(kt)
      vksq=ek*ek+fk*fk
      ikr=0.0
      iki=0.0

C     Process branches

      kp = 0
      do l=km(kt), km(kt)-1+kmlen(kt)
         mt=ikmu(l)
         lp=lp+1
         if (kp .eq. 0)   then
           if (mt .ge. kt)   then
              kp = l
              kl = lp
              kolum(lp) = kt
              korder(lp) = lp +1
              lp = lp +1
           endif
         endif
         em=e(mt)
         fm=f(mt)
         g12 = gkmu(l)
         b12 = bkmu(l)
         imr=em*g12-fm*b12
         imi=em*b12+fm*g12
         ikr=ikr+imr
         iki=iki+imi
         kolum(lp)=mt
         korder(lp)=lp+1
         rh=imr*fk-imi*ek
         rn=imr*ek+imi*fk
         rj=-rn
         rl=rh
         rowh(lp)=rh
         rowj(lp)=rj
         rown(lp)=rn
         rowl(lp)=rl
      enddo

C     Process diagonal

      if (kp.eq.0) then
         kl=lp+1
         lp=kl
         kolum(lp)=kt
         mt=kt
      endif

C     Add identity row for passive BM nodes

      if (lp.eq.1) then
         rowh(kl) = 1.0
         rown(kl) = 0.0
         rowj(kl) = 0.0
         rowl(kl) = 1.0
      else
         g12 = gkku(kt)
         b12 = bkku(kt)
         imr = ek*g12 - fk*b12
         imi = ek*b12 + fk*g12
         ikr = ikr + imr
         iki = iki + imi

         dvk=sqrt(vksq)
         imr=inetr(kt)*dvk
         imi=ineti(kt)*dvk
         dpk=ikr*ek+iki*fk+imr
         dqk=ikr*fk-iki*ek-imi
         pk = dpk
         qk = dqk
         g12 = gkku(kt)*vksq
         b12 = bkku(kt)*vksq
         rh=-dqk-b12-imi
         rn=dpk+g12
         rj=dpk-g12-imr
         rl=dqk-b12
         rowh(kl)=rh
         rown(kl)=rn
         rowj(kl)=rj
         rowl(kl)=rl
      endif
      korder(lp)=0
      max=mt
      mend=lp+1
      ko=lp
      korder(0)=1
      if (ksw.ne.0) return
  294 mt = senchk(kt)

C     Set up contraints for PV node

      if (mt.ne.0) then
          do 300 i = 1,mend-1
            rowj(i) = 0.0
            rowl(i) = 0.0
  300     continue
          if (mt.eq.kt) then

C            PV node

             rowl(kl) = 1.0
          else

C            PV node controlling MT

             ip = 0
             ko = korder(ip)
  302        if (kolum(ko)-mt) 304,308,306
  304        ip=ko
             ko=korder(ko)
             if (ko.ne.0) go to 302
             max = mt
             lp=mend
  306        korder(mend)=korder(ip)
             korder(ip)=mend
             ko=mend
             mend=mend+1
             kolum(ko)=mt
             rowh(ko)=0.0
             rowj(ko)=0.0
             rown(ko)=0.0
             rowl(ko)=0.0
  308        rowl(ko)=1.0
             rowl(kl) = 0.0001
          endif
          if (isw.eq.0) go to 480
          lp = 1
          ko = 0
          go to 500
      endif

C     Implement phase shifter constraints

      do i = i1,i2
         if (jflag(1,i) .eq. 11 .or. jflag(1,i) .eq. 12) then
            if (jflag(1,i) .eq. 11) then
               ax = 1.0
            else
               ax = -1.0
            endif
            ji = jflag(2,i)
            mt = ntot + ntopt + dabs(txtie(8,ji))
            max = mt
            korder(lp) = mend
            lp = mend
            mend = mend+1
            kolum(lp) = mt
            rowh(lp) = ax
            rowj(lp) = 0.0
            rown(lp) = 0.0
            rowl(lp) = ax
            korder(lp) = 0
            go to 320
         endif
      enddo
  320 continue

  480 if (isw .eq. 0) go to 780
C
C     Nullify P constraint in Jacobian. It will be replaced with
C     an Area Interchange constraint.
C
      do i = 1,mend-1
         rowh(i) = 0.0
         rown(i) = 0.0
      enddo
      lp=korder(0)
      ko=0

C     Define area interchange constraints

  500 xatot=0.0
      jt=isw
      if(jt.eq.0) call erexit
      j1=karea(3,jt)
      js=karea(4,jt)+j1-1
      do 700 j=j1,js
      ix=kaloc(j)
      iax=iabs(ix)
      k1=tie(1,iax)
      k2=tie(7,iax)
      mt=k1
      if (ix.lt.0) mt=k2
      ka1=tie(2,iax)
      kdc=tie(9,iax)
      if (kdc.eq.0) go to 580
      if (ix.lt.0) go to 700
      kd = kdc
  530 k1x = dmin1 (dcline(1,kd),dcline(2,kd))
      k2x = dmax1 (dcline(1,kd),dcline(2,kd))
      if (k1x.ne.min0(k1,k2)) then
         if (kd.ne.kdc) call erexit
         if (mtdcln.eq.0) call erexit
         kd = kdc + mtdcln
         go to 530
      else if (k2x.ne.max0(k1,k2)) then
         call erexit
      endif
      if (k1.eq.dcline(1,kd)) then
         l1 = dcline(8,kd)
         l2 = dcline(9,kd)
      else
         l1 = dcline(9,kd)
         l2 = dcline(8,kd)
      endif
      v1=dcbus(20,l1)
      v2=dcbus(20,l2)
      pin = v1*(v1 - v2)/(dcline(4,kd)*bmva)
      if (jt.ne.ka1) pin = -pin
      xatot = xatot + pin
      if (idswb.ne.0) write (dbug,570)kt,pin
  570 format(' DC export from area slack node ',i4,' is ',f7.3)
      go to 700

  580 continue
  590 if (kolum(lp)-mt) 600,620,610
  620 ko=lp
      go to 630

  600 ko=lp
      lp=korder(lp)
      if (lp.ne.0) go to 590
      max=mt
  610 korder(mend)=korder(ko)
      korder(ko)=mend
      kolum(mend)=mt
      ko=mend
      lp=mend
      mend=mend+1
      rowh(lp)=0.0
      rown(lp)=0.0
      rowj(lp)=0.0
      rowl(lp)=0.0
  630 ek=e(k1)
      fk=f(k1)
      em=e(k2)
      fm=f(k2)
      vksq=ek*ek+fk*fk
      g12=tie(5,iax)
      b12=tie(6,iax)
      ikr=g12*em-b12*fm
      iki=b12*em+g12*fm
      rh=-ek*iki+fk*ikr
      rn=ek*ikr+fk*iki
      pin=rn+vksq*tie(3,iax)
      if (ix.gt.0) then
         rh=-rh
         rn=rn+2.0*vksq*tie(3,iax)
      endif
      if (ka1.ne.jt) then
         pin=-pin
         rh=-rh
         rn=-rn
      endif
      rowh(lp)=rowh(lp)+rh
      rown(lp)=rown(lp)+rn
  690 if (ix.gt.0) xatot=xatot+pin
  700 continue
  710 if (kolum(lp)-max) 720,740,730
  720 lp=korder(lp)
      go to 710

  730 call erexit

  740 continue
C
C     The following logic projects a target PNET for
c     the area slack bus such that the ensuing convergence
c     test in SWITCH is applicable.
C
      if (idswc.gt.0) then
         write (dbug,760) jt,kt,area(2,jt),xatot
  760    format ('0 Area interchange ',2i6,4e20.6)
         l = 0
         write (dbug,770) l,kolum(l),korder(l)
  770    format ('  Debug of area Jacobian '/(3i8,4e20.6))
         ix=mend-1
         write (dbug,772) (l,kolum(l),korder(l),rowh(l),rown(l),rowj(l),
     1    rowl(l),l=1,ix)
  772    format (3i8,4e20.6)
      endif
  780 continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2015 NCAR/RAL
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! 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
! (at your option) 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 diagn_evar_module
! data types
USE nrtype

! physical constants
USE multiconst,only:&
                    iden_air,    & ! intrinsic density of air      (kg m-3)
                    iden_ice,    & ! intrinsic density of ice      (kg m-3)
                    iden_water,  & ! intrinsic density of water    (kg m-3)
                    ! specific heat
                    Cp_air,      & ! specific heat of air          (J kg-1 K-1)
                    Cp_ice,      & ! specific heat of ice          (J kg-1 K-1)
                    Cp_soil,     & ! specific heat of soil         (J kg-1 K-1)
                    Cp_water,    & ! specific heat of liquid water (J kg-1 K-1)
                    ! thermal conductivity
                    lambda_air,  & ! thermal conductivity of air   (J s-1 m-1)
                    lambda_ice,  & ! thermal conductivity of ice   (J s-1 m-1)
                    lambda_water   ! thermal conductivity of water (J s-1 m-1)

! access missing values
USE globalData,only:integerMissing  ! missing integer
USE globalData,only:realMissing     ! missing real number

! named variables that define the layer type
USE globalData,only:iname_snow     ! snow
USE globalData,only:iname_soil     ! soil

! model decisions
USE mDecisions_module,only:Smirnova2000  ! option for temporally constant thermal conductivity

implicit none
private
public::diagn_evar
! algorithmic parameters
real(dp),parameter     :: valueMissing=-9999._dp  ! missing value, used when diagnostic or state variables are undefined
real(dp),parameter     :: verySmall=1.e-6_dp   ! used as an additive constant to check if substantial difference among real numbers
real(dp),parameter     :: mpe=1.e-6_dp         ! prevents overflow error if division by zero
real(dp),parameter     :: dx=1.e-6_dp          ! finite difference increment
contains

 ! **********************************************************************************************************
 ! public subroutine diagn_evar: compute diagnostic energy variables (thermal conductivity and heat capacity)
 ! **********************************************************************************************************
 subroutine diagn_evar(&
                       ! input: control variables
                       computeVegFlux,          & ! intent(in): flag to denote if computing the vegetation flux
                       canopyDepth,             & ! intent(in): canopy depth (m)
                       ! input/output: data structures
                       mpar_data,               & ! intent(in):    model parameters
                       indx_data,               & ! intent(in):    model layer indices
                       prog_data,               & ! intent(in):    model prognostic variables for a local HRU
                       diag_data,               & ! intent(inout): model diagnostic variables for a local HRU
                       ! output: error control
                       err,message)               ! intent(out): error control
 ! --------------------------------------------------------------------------------------------------------------------------------------
 ! provide access to the derived types to define the data structures
 USE data_types,only:&
                     var_d,            & ! data vector (dp)
                     var_ilength,      & ! data vector with variable length dimension (i4b)
                     var_dlength         ! data vector with variable length dimension (dp)
 ! provide access to named variables defining elements in the data structures
 USE var_lookup,only:iLookPARAM,iLookPROG,iLookDIAG,iLookINDEX  ! named variables for structure elements
 USE var_lookup,only:iLookDECISIONS               ! named variables for elements of the decision structure
 ! provide access to named variables for thermal conductivity of soil
 USE globalData,only:model_decisions        ! model decision structure
 USE mDecisions_module,only: funcSoilWet, & ! function of soil wetness
                             mixConstit,  & ! mixture of constituents
                             hanssonVZJ     ! test case for the mizoguchi lab experiment, Hansson et al. VZJ 2004
 ! provide access to external subroutines
 USE snow_utils_module,only:tcond_snow            ! compute thermal conductivity of snow
 ! --------------------------------------------------------------------------------------------------------------------------------------
 ! input: model control
 logical(lgt),intent(in)         :: computeVegFlux         ! logical flag to denote if computing the vegetation flux
 real(dp),intent(in)             :: canopyDepth            ! depth of the vegetation canopy (m)
 ! input/output: data structures
 type(var_dlength),intent(in)    :: mpar_data              ! model parameters
 type(var_ilength),intent(in)    :: indx_data              ! model layer indices
 type(var_dlength),intent(in)    :: prog_data              ! model prognostic variables for a local HRU
 type(var_dlength),intent(inout) :: diag_data              ! model diagnostic variables for a local HRU
 ! output: error control
 integer(i4b),intent(out)        :: err                    ! error code
 character(*),intent(out)        :: message                ! error message
 ! --------------------------------------------------------------------------------------------------------------------------------
 ! local variables
 character(LEN=256)                :: cmessage               ! error message of downwind routine
 integer(i4b)                      :: iLayer                 ! index of model layer
 integer(i4b)                      :: iSoil                  ! index of soil layer
 real(dp)                          :: TCn                    ! thermal conductivity below the layer interface (W m-1 K-1)
 real(dp)                          :: TCp                    ! thermal conductivity above the layer interface (W m-1 K-1)
 real(dp)                          :: zdn                    ! height difference between interface and lower value (m)
 real(dp)                          :: zdp                    ! height difference between interface and upper value (m)
 real(dp)                          :: bulkden_soil           ! bulk density of soil (kg m-3)
 real(dp)                          :: lambda_drysoil         ! thermal conductivity of dry soil (W m-1)
 real(dp)                          :: lambda_wetsoil         ! thermal conductivity of wet soil (W m-1)
 real(dp)                          :: lambda_wet             ! thermal conductivity of the wet material
 real(dp)                          :: relativeSat            ! relative saturation (-)
 real(dp)                          :: kerstenNum             ! the Kersten number (-), defining weight applied to conductivity of the wet medium
 real(dp)                          :: den                    ! denominator in the thermal conductivity calculations
 ! local variables to reproduce the thermal conductivity of Hansson et al. VZJ 2005
 real(dp),parameter                :: c1=0.55_dp             ! optimized parameter from Hansson et al. VZJ 2005 (W m-1 K-1)
 real(dp),parameter                :: c2=0.8_dp              ! optimized parameter from Hansson et al. VZJ 2005 (W m-1 K-1)
 real(dp),parameter                :: c3=3.07_dp             ! optimized parameter from Hansson et al. VZJ 2005 (-)
 real(dp),parameter                :: c4=0.13_dp             ! optimized parameter from Hansson et al. VZJ 2005 (W m-1 K-1)
 real(dp),parameter                :: c5=4._dp               ! optimized parameter from Hansson et al. VZJ 2005 (-)
 real(dp),parameter                :: f1=13.05_dp            ! optimized parameter from Hansson et al. VZJ 2005 (-)
 real(dp),parameter                :: f2=1.06_dp             ! optimized parameter from Hansson et al. VZJ 2005 (-)
 real(dp)                          :: fArg,xArg              ! temporary variables (see Hansson et al. VZJ 2005 for details)
 ! --------------------------------------------------------------------------------------------------------------------------------
 ! associate variables in data structure
 associate(&
 ! input: model decisions
 ixThCondSnow            => model_decisions(iLookDECISIONS%thCondSnow)%iDecision,      & ! intent(in): choice of method for thermal conductivity of snow
 ixThCondSoil            => model_decisions(iLookDECISIONS%thCondSoil)%iDecision,      & ! intent(in): choice of method for thermal conductivity of soil
 ! input: state variables
 scalarCanopyIce         => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1),           & ! intent(in): canopy ice content (kg m-2)
 scalarCanopyLiquid      => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1),           & ! intent(in): canopy liquid water content (kg m-2)
 mLayerVolFracIce        => prog_data%var(iLookPROG%mLayerVolFracIce)%dat,             & ! intent(in): volumetric fraction of ice at the start of the sub-step (-)
 mLayerVolFracLiq        => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat,             & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
 ! input: coordinate variables
 nSnow                   => indx_data%var(iLookINDEX%nSnow)%dat(1),                    & ! intent(in): number of snow layers
 nSoil                   => indx_data%var(iLookINDEX%nSoil)%dat(1),                    & ! intent(in): number of soil layers
 nLayers                 => indx_data%var(iLookINDEX%nLayers)%dat(1),                  & ! intent(in): total number of layers
 layerType               => indx_data%var(iLookINDEX%layerType)%dat,                   & ! intent(in): layer type (iname_soil or iname_snow)
 mLayerHeight            => prog_data%var(iLookPROG%mLayerHeight)%dat,                 & ! intent(in): height at the mid-point of each layer (m)
 iLayerHeight            => prog_data%var(iLookPROG%iLayerHeight)%dat,                 & ! intent(in): height at the interface of each layer (m)
 ! input: heat capacity and thermal conductivity
 specificHeatVeg         => mpar_data%var(iLookPARAM%specificHeatVeg)%dat(1),          & ! intent(in): specific heat of vegetation (J kg-1 K-1)
 maxMassVegetation       => mpar_data%var(iLookPARAM%maxMassVegetation)%dat(1),        & ! intent(in): maximum mass of vegetation (kg m-2)
 fixedThermalCond_snow   => mpar_data%var(iLookPARAM%fixedThermalCond_snow)%dat(1),    & ! intent(in): temporally constant thermal conductivity of snow (W m-1 K-1)
 ! input: depth varying soil parameters
 iden_soil               => mpar_data%var(iLookPARAM%soil_dens_intr)%dat,              & ! intent(in): intrinsic density of soil (kg m-3)
 thCond_soil             => mpar_data%var(iLookPARAM%thCond_soil)%dat,                 & ! intent(in): thermal conductivity of soil (W m-1 K-1)
 theta_sat               => mpar_data%var(iLookPARAM%theta_sat)%dat,                   & ! intent(in): soil porosity (-)
 frac_sand               => mpar_data%var(iLookPARAM%frac_sand)%dat,                   & ! intent(in): fraction of sand (-)
 frac_silt               => mpar_data%var(iLookPARAM%frac_silt)%dat,                   & ! intent(in): fraction of silt (-)
 frac_clay               => mpar_data%var(iLookPARAM%frac_clay)%dat,                   & ! intent(in): fraction of clay (-)
 ! output: diagnostic variables
 scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1),   & ! intent(out): volumetric heat capacity of the vegetation (J m-3 K-1)
 mLayerVolHtCapBulk      => diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat,           & ! intent(out): volumetric heat capacity in each layer (J m-3 K-1)
 mLayerThermalC          => diag_data%var(iLookDIAG%mLayerThermalC)%dat,               & ! intent(out): thermal conductivity at the mid-point of each layer (W m-1 K-1)
 iLayerThermalC          => diag_data%var(iLookDIAG%iLayerThermalC)%dat,               & ! intent(out): thermal conductivity at the interface of each layer (W m-1 K-1)
 mLayerVolFracAir        => diag_data%var(iLookDIAG%mLayerVolFracAir)%dat              & ! intent(out): volumetric fraction of air in each layer (-)
 )  ! end associate statement
 ! --------------------------------------------------------------------------------------------------------------------------------
 ! initialize error control
 err=0; message="diagn_evar/"

 ! initialize the soil layer
 iSoil=integerMissing

 ! compute the bulk volumetric heat capacity of vegetation (J m-3 K-1)
 if(computeVegFlux)then
  scalarBulkVolHeatCapVeg = specificHeatVeg*maxMassVegetation/canopyDepth + & ! vegetation component
                            Cp_water*scalarCanopyLiquid/canopyDepth       + & ! liquid water component
                            Cp_ice*scalarCanopyIce/canopyDepth                ! ice component
 else
  scalarBulkVolHeatCapVeg = valueMissing
 end if
 !print*, 'diagn_evar: scalarBulkVolHeatCapVeg = ', scalarBulkVolHeatCapVeg

 ! loop through layers
 do iLayer=1,nLayers

  ! get the soil layer
  if(iLayer>nSnow) iSoil = iLayer-nSnow

  ! compute the thermal conductivity of dry and wet soils (W m-1)
  ! NOTE: this is actually constant over the simulation, and included here for clarity
  if(ixThCondSoil == funcSoilWet .and. layerType(iLayer)==iname_soil)then
   bulkden_soil   = iden_soil(iSoil)*( 1._dp - theta_sat(iSoil) )
   lambda_drysoil = (0.135_dp*bulkden_soil + 64.7_dp) / (iden_soil(iSoil) - 0.947_dp*bulkden_soil)
   lambda_wetsoil = (8.80_dp*frac_sand(iSoil) + 2.92_dp*frac_clay(iSoil)) / (frac_sand(iSoil) + frac_clay(iSoil))
  end if

  ! *****
  ! * compute the volumetric fraction of air in each layer...
  ! *********************************************************
  select case(layerType(iLayer))
   case(iname_soil); mLayerVolFracAir(iLayer) = theta_sat(iSoil) - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer))
   case(iname_snow); mLayerVolFracAir(iLayer) = 1._dp - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer))
   case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute volumetric fraction of air'; return
  end select

  ! *****
  ! * compute the volumetric heat capacity of each layer (J m-3 K-1)...
  ! *******************************************************************
  select case(layerType(iLayer))
   ! * soil
   case(iname_soil)
    mLayerVolHtCapBulk(iLayer) = iden_soil(iSoil)  * Cp_soil  * ( 1._dp - theta_sat(iSoil) ) + & ! soil component
                                 iden_ice          * Cp_Ice   * mLayerVolFracIce(iLayer)     + & ! ice component
                                 iden_water        * Cp_water * mLayerVolFracLiq(iLayer)     + & ! liquid water component
                                 iden_air          * Cp_air   * mLayerVolFracAir(iLayer)         ! air component
   ! * snow
   case(iname_snow)
    mLayerVolHtCapBulk(iLayer) = iden_ice          * Cp_ice   * mLayerVolFracIce(iLayer)     + & ! ice component
                                 iden_water        * Cp_water * mLayerVolFracLiq(iLayer)     + & ! liquid water component
                                 iden_air          * Cp_air   * mLayerVolFracAir(iLayer)         ! air component
   case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute olumetric heat capacity'; return
  end select

  ! *****
  ! * compute the thermal conductivity of snow and soil at the mid-point of each layer...
  ! *************************************************************************************
  select case(layerType(iLayer))

   ! ***** soil
   case(iname_soil)

    ! select option for thermal conductivity of soil
    select case(ixThCondSoil)

     ! ** function of soil wetness
     case(funcSoilWet)

      ! compute the thermal conductivity of the wet material (W m-1)
      lambda_wet  = lambda_wetsoil**( 1._dp - theta_sat(iSoil) ) * lambda_water**theta_sat(iSoil) * lambda_ice**(theta_sat(iSoil) - mLayerVolFracLiq(iLayer))
      relativeSat = (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer))/theta_sat(iSoil)  ! relative saturation
      ! compute the Kersten number (-)
      if(relativeSat > 0.1_dp)then ! log10(0.1) = -1
       kerstenNum = log10(relativeSat) + 1._dp
      else
       kerstenNum = 0._dp  ! dry thermal conductivity
      endif
      ! ...and, compute the thermal conductivity
      mLayerThermalC(iLayer) = kerstenNum*lambda_wet + (1._dp - kerstenNum)*lambda_drysoil

     ! ** mixture of constituents
     case(mixConstit)
      mLayerThermalC(iLayer) = thCond_soil(iSoil) * ( 1._dp - theta_sat(iSoil) ) + & ! soil component
                               lambda_ice         * mLayerVolFracIce(iLayer)     + & ! ice component
                               lambda_water       * mLayerVolFracLiq(iLayer)     + & ! liquid water component
                               lambda_air         * mLayerVolFracAir(iLayer)         ! air component

     ! ** test case for the mizoguchi lab experiment, Hansson et al. VZJ 2004
     case(hanssonVZJ)
      fArg  = 1._dp + f1*mLayerVolFracIce(iLayer)**f2
      xArg  = mLayerVolFracLiq(iLayer) + fArg*mLayerVolFracIce(iLayer)
      mLayerThermalC(iLayer) = c1 + c2*xArg + (c1 - c4)*exp(-(c3*xArg)**c5)

     ! ** check
     case default; err=20; message=trim(message)//'unable to identify option for thermal conductivity of soil'; return

    end select  ! option for the thermal conductivity of soil

   ! ***** snow
   case(iname_snow)
    ! temporally constant thermal conductivity
    if(ixThCondSnow==Smirnova2000)then
     mLayerThermalC(iLayer) = fixedThermalCond_snow
    ! thermal conductivity as a function of snow density
    else
     call tcond_snow(mLayerVolFracIce(iLayer)*iden_ice,  & ! input: snow density (kg m-3)
                     mLayerThermalC(iLayer),             & ! output: thermal conductivity (W m-1 K-1)
                     err,cmessage)                         ! output: error control
     if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
    endif

   ! * error check
   case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute thermal conductivity'; return

  end select
  !print*, 'iLayer, mLayerThermalC(iLayer) = ', iLayer, mLayerThermalC(iLayer)

 end do  ! looping through layers
 !pause

 ! *****
 ! * compute the thermal conductivity of snow at the interface of each layer...
 ! ****************************************************************************
 do iLayer=1,nLayers-1  ! (loop through layers)
  ! get temporary variables
  TCn = mLayerThermalC(iLayer)    ! thermal conductivity below the layer interface (W m-1 K-1)
  TCp = mLayerThermalC(iLayer+1)  ! thermal conductivity above the layer interface (W m-1 K-1)
  zdn = iLayerHeight(iLayer)   - mLayerHeight(iLayer) ! height difference between interface and lower value (m)
  zdp = mLayerHeight(iLayer+1) - iLayerHeight(iLayer) ! height difference between interface and upper value (m)
  den = TCn*zdp + TCp*zdn  ! denominator
  ! compute thermal conductivity
  if(TCn+TCp > epsilon(TCn))then
   iLayerThermalC(iLayer) = (TCn*TCp*(zdn + zdp)) / den
  else
   iLayerThermalC(iLayer) = (TCn*zdn +  TCp*zdp) / (zdn + zdp)
  endif
  !write(*,'(a,1x,i4,1x,10(f9.3,1x))') 'iLayer, TCn, TCp, zdn, zdp, iLayerThermalC(iLayer) = ', iLayer, TCn, TCp, zdn, zdp, iLayerThermalC(iLayer)
 end do  ! looping through layers

 ! special case of hansson
 if(ixThCondSoil==hanssonVZJ)then
  iLayerThermalC(0) = 28._dp*(0.5_dp*(iLayerHeight(1) - iLayerHeight(0)))
 else
  iLayerThermalC(0) = mLayerThermalC(1)
 end if

 ! assume the thermal conductivity at the domain boundaries is equal to the thermal conductivity of the layer
 iLayerThermalC(nLayers) = mLayerThermalC(nLayers)

 ! end association to variables in the data structure
 end associate

 end subroutine diagn_evar

end module diagn_evar_module

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
*DECK HWSPLR
      SUBROUTINE HWSPLR (A, B, M, MBDCND, BDA, BDB, C, D, N, NBDCND,
     +   BDC, BDD, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
C***BEGIN PROLOGUE  HWSPLR
C***PURPOSE  Solve a finite difference approximation to the Helmholtz
C            equation in polar coordinates.
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B1A1A
C***TYPE      SINGLE PRECISION (HWSPLR-S)
C***KEYWORDS  ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, POLAR
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C     Subroutine HWSPLR solves a finite difference approximation to the
C     Helmholtz equation in polar coordinates:
C
C          (1/R)(d/dR)(R(dU/dR)) + (1/R**2)(d/dTHETA)(dU/dTHETA)
C
C                                + LAMBDA*U = F(R,THETA).
C
C
C
C
C     * * * * * * * *    Parameter Description     * * * * * * * * * *
C
C             * * * * * *   On Input    * * * * * *
C
C     A,B
C       The range of R, i.e., A .LE. R .LE. B.  A must be less than B
C       and A must be non-negative.
C
C     M
C       The number of panels into which the interval (A,B) is
C       subdivided.  Hence, there will be M+1 grid points in the
C       R-direction given by R(I) = A+(I-1)DR, for I = 1,2,...,M+1,
C       where DR = (B-A)/M is the panel width. M must be greater than 3.
C
C     MBDCND
C       Indicates the type of boundary condition at R = A and R = B.
C
C       = 1  If the solution is specified at R = A and R = B.
C       = 2  If the solution is specified at R = A and the derivative of
C            the solution with respect to R is specified at R = B.
C       = 3  If the derivative of the solution with respect to R is
C            specified at R = A (see note below) and R = B.
C       = 4  If the derivative of the solution with respect to R is
C            specified at R = A (see note below) and the solution is
C            specified at R = B.
C       = 5  If the solution is unspecified at R = A = 0 and the
C            solution is specified at R = B.
C       = 6  If the solution is unspecified at R = A = 0 and the
C            derivative of the solution with respect to R is specified
C            at R = B.
C
C       NOTE:  If A = 0, do not use MBDCND = 3 or 4, but instead use
C              MBDCND = 1,2,5, or 6  .
C
C     BDA
C       A one-dimensional array of length N+1 that specifies the values
C       of the derivative of the solution with respect to R at R = A.
C       When MBDCND = 3 or 4,
C
C            BDA(J) = (d/dR)U(A,THETA(J)), J = 1,2,...,N+1  .
C
C       When MBDCND has any other value, BDA is a dummy variable.
C
C     BDB
C       A one-dimensional array of length N+1 that specifies the values
C       of the derivative of the solution with respect to R at R = B.
C       When MBDCND = 2,3, or 6,
C
C            BDB(J) = (d/dR)U(B,THETA(J)), J = 1,2,...,N+1  .
C
C       When MBDCND has any other value, BDB is a dummy variable.
C
C     C,D
C       The range of THETA, i.e., C .LE. THETA .LE. D.  C must be less
C       than D.
C
C     N
C       The number of panels into which the interval (C,D) is
C       subdivided.  Hence, there will be N+1 grid points in the
C       THETA-direction given by THETA(J) = C+(J-1)DTHETA for
C       J = 1,2,...,N+1, where DTHETA = (D-C)/N is the panel width.  N
C       must be greater than 3.
C
C     NBDCND
C       Indicates the type of boundary conditions at THETA = C and
C       at THETA = D.
C
C       = 0  If the solution is periodic in THETA, i.e.,
C            U(I,J) = U(I,N+J).
C       = 1  If the solution is specified at THETA = C and THETA = D
C            (see note below).
C       = 2  If the solution is specified at THETA = C and the
C            derivative of the solution with respect to THETA is
C            specified at THETA = D (see note below).
C       = 4  If the derivative of the solution with respect to THETA is
C            specified at THETA = C and the solution is specified at
C            THETA = D (see note below).
C
C       NOTE:  When NBDCND = 1,2, or 4, do not use MBDCND = 5 or 6
C              (the former indicates that the solution is specified at
C              R = 0, the latter indicates the solution is unspecified
C              at R = 0).  Use instead MBDCND = 1 or 2  .
C
C     BDC
C       A one-dimensional array of length M+1 that specifies the values
C       of the derivative of the solution with respect to THETA at
C       THETA = C.  When NBDCND = 3 or 4,
C
C            BDC(I) = (d/dTHETA)U(R(I),C), I = 1,2,...,M+1  .
C
C       When NBDCND has any other value, BDC is a dummy variable.
C
C     BDD
C       A one-dimensional array of length M+1 that specifies the values
C       of the derivative of the solution with respect to THETA at
C       THETA = D.  When NBDCND = 2 or 3,
C
C            BDD(I) = (d/dTHETA)U(R(I),D), I = 1,2,...,M+1  .
C
C       When NBDCND has any other value, BDD is a dummy variable.
C
C     ELMBDA
C       The constant LAMBDA in the Helmholtz equation.  If
C       LAMBDA .LT. 0, a solution may not exist.  However, HWSPLR will
C       attempt to find a solution.
C
C     F
C       A two-dimensional array that specifies the values of the right
C       side of the Helmholtz equation and boundary values (if any).
C       For I = 2,3,...,M and J = 2,3,...,N
C
C            F(I,J) = F(R(I),THETA(J)).
C
C       On the boundaries F is defined by
C
C            MBDCND   F(1,J)            F(M+1,J)
C            ------   -------------     -------------
C
C              1      U(A,THETA(J))     U(B,THETA(J))
C              2      U(A,THETA(J))     F(B,THETA(J))
C              3      F(A,THETA(J))     F(B,THETA(J))
C              4      F(A,THETA(J))     U(B,THETA(J))   J = 1,2,...,N+1
C              5      F(0,0)            U(B,THETA(J))
C              6      F(0,0)            F(B,THETA(J))
C
C            NBDCND   F(I,1)            F(I,N+1)
C            ------   ---------         ---------
C
C              0      F(R(I),C)         F(R(I),C)
C              1      U(R(I),C)         U(R(I),D)
C              2      U(R(I),C)         F(R(I),D)   I = 1,2,...,M+1
C              3      F(R(I),C)         F(R(I),D)
C              4      F(R(I),C)         U(R(I),D)
C
C       F must be dimensioned at least (M+1)*(N+1).
C
C       NOTE
C
C       If the table calls for both the solution U and the right side F
C       at a corner then the solution must be specified.
C
C
C     IDIMF
C       The row (or first) dimension of the array F as it appears in the
C       program calling HWSPLR.  This parameter is used to specify the
C       variable dimension of F.  IDIMF must be at least M+1  .
C
C     W
C       A one-dimensional array that must be provided by the user for
C       work space.  W may require up to 4*(N+1) +
C       (13 + INT(log2(N+1)))*(M+1) locations.  The actual number of
C       locations used is computed by HWSPLR and is returned in location
C       W(1).
C
C
C             * * * * * *   On Output     * * * * * *
C
C     F
C       Contains the solution U(I,J) of the finite difference
C       approximation for the grid point (R(I),THETA(J)),
C       I = 1,2,...,M+1, J = 1,2,...,N+1  .
C
C     PERTRB
C       If a combination of periodic, derivative, or unspecified
C       boundary conditions is specified for a Poisson equation
C       (LAMBDA = 0), a solution may not exist.  PERTRB is a constant,
C       calculated and subtracted from F, which ensures that a solution
C       exists.  HWSPLR then computes this solution, which is a least
C       squares solution to the original approximation.  This solution
C       plus any constant is also a solution.  Hence, the solution is
C       not unique.  PERTRB should be small compared to the right side.
C       Otherwise, a solution is obtained to an essentially different
C       problem.  This comparison should always be made to insure that a
C       meaningful solution has been obtained.
C
C     IERROR
C       An error flag that indicates invalid input parameters.  Except
C       for numbers 0 and 11, a solution is not attempted.
C
C       =  0  No error.
C       =  1  A .LT. 0  .
C       =  2  A .GE. B.
C       =  3  MBDCND .LT. 1 or MBDCND .GT. 6  .
C       =  4  C .GE. D.
C       =  5  N .LE. 3
C       =  6  NBDCND .LT. 0 or .GT. 4  .
C       =  7  A = 0, MBDCND = 3 or 4  .
C       =  8  A .GT. 0, MBDCND .GE. 5  .
C       =  9  MBDCND .GE. 5, NBDCND .NE. 0 and NBDCND .NE. 3  .
C       = 10  IDIMF .LT. M+1  .
C       = 11  LAMBDA .GT. 0  .
C       = 12  M .LE. 3
C
C       Since this is the only means of indicating a possibly incorrect
C       call to HWSPLR, the user should test IERROR after the call.
C
C     W
C       W(1) contains the required length of W.
C
C *Long Description:
C
C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   BDA(N+1),BDB(N+1),BDC(M+1),BDD(M+1),F(IDIMF,N+1),
C     Arguments      W(see argument list)
C
C     Latest         June 1, 1976
C     Revision
C
C     Subprograms    HWSPLR,GENBUN,POISD2,POISN2,POISP2,COSGEN,MERGE,
C     Required       TRIX,TRI3,PIMACH
C
C     Special        None
C     Conditions
C
C     Common         NONE
C     Blocks
C
C     I/O
C
C     Precision      Single
C
C     Specialist     Roland Sweet
C
C     Language       FORTRAN
C
C     History        Standardized April 1, 1973
C                    Revised January 1, 1976
C
C     Algorithm      The routine defines the finite difference
C                    equations, incorporates boundary data, and adjusts
C                    the right side of singular systems and then calls
C                    GENBUN to solve the system.
C
C     Space          13430(octal) = 5912(decimal)  locations on the NCAR
C     Required       Control Data 7600
C
C     Timing and        The execution time T on the NCAR Control Data
C     Accuracy       7600 for subroutine HWSPLR is roughly proportional
C                    to M*N*log2(N), but also depends on the input
C                    parameters NBDCND and MBDCND.  Some typical values
C                    are listed in the table below.
C                       The solution process employed results in a loss
C                    of no more than three significant digits for N and
C                    M as large as 64.  More detailed information about
C                    accuracy can be found in the documentation for
C                    subroutine GENBUN which is the routine that
C                    solves the finite difference equations.
C
C
C                       M(=N)    MBDCND    NBDCND    T(MSECS)
C                       -----    ------    ------    --------
C
C                        32        1         0          31
C                        32        1         1          23
C                        32        3         3          36
C                        64        1         0         128
C                        64        1         1          96
C                        64        3         3         142
C
C     Portability    American National Standards Institute FORTRAN.
C                    The machine dependent constant PI is defined in
C                    function PIMACH.
C
C     Required       COS
C     Resident
C     Routines
C
C     Reference      Swarztrauber, P. and R. Sweet, 'Efficient FORTRAN
C                    Subprograms For The Solution Of Elliptic Equations'
C                    NCAR TN/IA-109, July, 1975, 138 pp.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES  P. N. Swarztrauber and R. Sweet, Efficient Fortran
C                 subprograms for the solution of elliptic equations,
C                 NCAR TN/IA-109, July 1975, 138 pp.
C***ROUTINES CALLED  GENBUN
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  HWSPLR
C
C
      DIMENSION       F(IDIMF,*)
      DIMENSION       BDA(*)     ,BDB(*)     ,BDC(*)     ,BDD(*)     ,
     1                W(*)
C***FIRST EXECUTABLE STATEMENT  HWSPLR
      IERROR = 0
      IF (A .LT. 0.) IERROR = 1
      IF (A .GE. B) IERROR = 2
      IF (MBDCND.LE.0 .OR. MBDCND.GE.7) IERROR = 3
      IF (C .GE. D) IERROR = 4
      IF (N .LE. 3) IERROR = 5
      IF (NBDCND.LE.-1 .OR. NBDCND.GE.5) IERROR = 6
      IF (A.EQ.0. .AND. (MBDCND.EQ.3 .OR. MBDCND.EQ.4)) IERROR = 7
      IF (A.GT.0. .AND. MBDCND.GE.5) IERROR = 8
      IF (MBDCND.GE.5 .AND. NBDCND.NE.0 .AND. NBDCND.NE.3) IERROR = 9
      IF (IDIMF .LT. M+1) IERROR = 10
      IF (M .LE. 3) IERROR = 12
      IF (IERROR .NE. 0) RETURN
      MP1 = M+1
      DELTAR = (B-A)/M
      DLRBY2 = DELTAR/2.
      DLRSQ = DELTAR**2
      NP1 = N+1
      DELTHT = (D-C)/N
      DLTHSQ = DELTHT**2
      NP = NBDCND+1
C
C     DEFINE RANGE OF INDICES I AND J FOR UNKNOWNS U(I,J).
C
      MSTART = 2
      MSTOP = MP1
      GO TO (101,105,102,103,104,105),MBDCND
  101 MSTOP = M
      GO TO 105
  102 MSTART = 1
      GO TO 105
  103 MSTART = 1
  104 MSTOP = M
  105 MUNK = MSTOP-MSTART+1
      NSTART = 1
      NSTOP = N
      GO TO (109,106,107,108,109),NP
  106 NSTART = 2
      GO TO 109
  107 NSTART = 2
  108 NSTOP = NP1
  109 NUNK = NSTOP-NSTART+1
C
C     DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
C
      ID2 = MUNK
      ID3 = ID2+MUNK
      ID4 = ID3+MUNK
      ID5 = ID4+MUNK
      ID6 = ID5+MUNK
      A1 = 2./DLRSQ
      IJ = 0
      IF (MBDCND.EQ.3 .OR. MBDCND.EQ.4) IJ = 1
      DO 110 I=1,MUNK
         R = A+(I-IJ)*DELTAR
         J = ID5+I
         W(J) = R
         J = ID6+I
         W(J) = 1./R**2
         W(I) = (R-DLRBY2)/(R*DLRSQ)
         J = ID3+I
         W(J) = (R+DLRBY2)/(R*DLRSQ)
         J = ID2+I
         W(J) = -A1+ELMBDA
  110 CONTINUE
      GO TO (114,111,112,113,114,111),MBDCND
  111 W(ID2) = A1
      GO TO 114
  112 W(ID2) = A1
  113 W(ID3+1) = A1
  114 CONTINUE
C
C     ENTER BOUNDARY DATA FOR R-BOUNDARIES.
C
      GO TO (115,115,117,117,119,119),MBDCND
  115 A1 = W(1)
      DO 116 J=NSTART,NSTOP
         F(2,J) = F(2,J)-A1*F(1,J)
  116 CONTINUE
      GO TO 119
  117 A1 = 2.*DELTAR*W(1)
      DO 118 J=NSTART,NSTOP
         F(1,J) = F(1,J)+A1*BDA(J)
  118 CONTINUE
  119 GO TO (120,122,122,120,120,122),MBDCND
  120 A1 = W(ID4)
      DO 121 J=NSTART,NSTOP
         F(M,J) = F(M,J)-A1*F(MP1,J)
  121 CONTINUE
      GO TO 124
  122 A1 = 2.*DELTAR*W(ID4)
      DO 123 J=NSTART,NSTOP
         F(MP1,J) = F(MP1,J)-A1*BDB(J)
  123 CONTINUE
C
C     ENTER BOUNDARY DATA FOR THETA-BOUNDARIES.
C
  124 A1 = 1./DLTHSQ
      L = ID5-MSTART+1
      LP = ID6-MSTART+1
      GO TO (134,125,125,127,127),NP
  125 DO 126 I=MSTART,MSTOP
         J = I+LP
         F(I,2) = F(I,2)-A1*W(J)*F(I,1)
  126 CONTINUE
      GO TO 129
  127 A1 = 2./DELTHT
      DO 128 I=MSTART,MSTOP
         J = I+LP
         F(I,1) = F(I,1)+A1*W(J)*BDC(I)
  128 CONTINUE
  129 A1 = 1./DLTHSQ
      GO TO (134,130,132,132,130),NP
  130 DO 131 I=MSTART,MSTOP
         J = I+LP
         F(I,N) = F(I,N)-A1*W(J)*F(I,NP1)
  131 CONTINUE
      GO TO 134
  132 A1 = 2./DELTHT
      DO 133 I=MSTART,MSTOP
         J = I+LP
         F(I,NP1) = F(I,NP1)-A1*W(J)*BDD(I)
  133 CONTINUE
  134 CONTINUE
C
C     ADJUST RIGHT SIDE OF EQUATION FOR UNKNOWN AT POLE WHEN HAVE
C     DERIVATIVE SPECIFIED BOUNDARY CONDITIONS.
C
      IF (MBDCND.GE.5 .AND. NBDCND.EQ.3)
     1    F(1,1) = F(1,1)-(BDD(2)-BDC(2))*4./(N*DELTHT*DLRSQ)
C
C     ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A
C     SOLUTION.
C
      PERTRB = 0.
      IF (ELMBDA) 144,136,135
  135 IERROR = 11
      GO TO 144
  136 IF (NBDCND.NE.0 .AND. NBDCND.NE.3) GO TO 144
      S2 = 0.
      GO TO (144,144,137,144,144,138),MBDCND
  137 W(ID5+1) = .5*(W(ID5+2)-DLRBY2)
      S2 = .25*DELTAR
  138 A2 = 2.
      IF (NBDCND .EQ. 0) A2 = 1.
      J = ID5+MUNK
      W(J) = .5*(W(J-1)+DLRBY2)
      S = 0.
      DO 140 I=MSTART,MSTOP
         S1 = 0.
         IJ = NSTART+1
         K = NSTOP-1
         DO 139 J=IJ,K
            S1 = S1+F(I,J)
  139    CONTINUE
         J = I+L
         S = S+(A2*S1+F(I,NSTART)+F(I,NSTOP))*W(J)
  140 CONTINUE
      S2 = M*A+DELTAR*((M-1)*(M+1)*.5+.25)+S2
      S1 = (2.+A2*(NUNK-2))*S2
      IF (MBDCND .EQ. 3) GO TO 141
      S2 = N*A2*DELTAR/8.
      S = S+F(1,1)*S2
      S1 = S1+S2
  141 CONTINUE
      PERTRB = S/S1
      DO 143 I=MSTART,MSTOP
         DO 142 J=NSTART,NSTOP
            F(I,J) = F(I,J)-PERTRB
  142    CONTINUE
  143 CONTINUE
  144 CONTINUE
C
C     MULTIPLY I-TH EQUATION THROUGH BY (R(I)*DELTHT)**2.
C
      DO 146 I=MSTART,MSTOP
         K = I-MSTART+1
         J = I+LP
         A1 = DLTHSQ/W(J)
         W(K) = A1*W(K)
         J = ID2+K
         W(J) = A1*W(J)
         J = ID3+K
         W(J) = A1*W(J)
         DO 145 J=NSTART,NSTOP
            F(I,J) = A1*F(I,J)
  145    CONTINUE
  146 CONTINUE
      W(1) = 0.
      W(ID4) = 0.
C
C     CALL GENBUN TO SOLVE THE SYSTEM OF EQUATIONS.
C
      CALL GENBUN (NBDCND,NUNK,1,MUNK,W(1),W(ID2+1),W(ID3+1),IDIMF,
     1             F(MSTART,NSTART),IERR1,W(ID4+1))
      IWSTOR = W(ID4+1)+3*MUNK
      GO TO (157,157,157,157,148,147),MBDCND
C
C     ADJUST THE SOLUTION AS NECESSARY FOR THE PROBLEMS WHERE A = 0.
C
  147 IF (ELMBDA .NE. 0.) GO TO 148
      YPOLE = 0.
      GO TO 155
  148 CONTINUE
      J = ID5+MUNK
      W(J) = W(ID2)/W(ID3)
      DO 149 IP=3,MUNK
         I = MUNK-IP+2
         J = ID5+I
         LP = ID2+I
         K = ID3+I
         W(J) = W(I)/(W(LP)-W(K)*W(J+1))
  149 CONTINUE
      W(ID5+1) = -.5*DLTHSQ/(W(ID2+1)-W(ID3+1)*W(ID5+2))
      DO 150 I=2,MUNK
         J = ID5+I
         W(J) = -W(J)*W(J-1)
  150 CONTINUE
      S = 0.
      DO 151 J=NSTART,NSTOP
         S = S+F(2,J)
  151 CONTINUE
      A2 = NUNK
      IF (NBDCND .EQ. 0) GO TO 152
      S = S-.5*(F(2,NSTART)+F(2,NSTOP))
      A2 = A2-1.
  152 YPOLE = (.25*DLRSQ*F(1,1)-S/A2)/(W(ID5+1)-1.+ELMBDA*DLRSQ*.25)
      DO 154 I=MSTART,MSTOP
         K = L+I
         DO 153 J=NSTART,NSTOP
            F(I,J) = F(I,J)+YPOLE*W(K)
  153    CONTINUE
  154 CONTINUE
  155 DO 156 J=1,NP1
         F(1,J) = YPOLE
  156 CONTINUE
  157 CONTINUE
      IF (NBDCND .NE. 0) GO TO 159
      DO 158 I=MSTART,MSTOP
         F(I,NP1) = F(I,1)
  158 CONTINUE
  159 CONTINUE
      W(1) = IWSTOR
      RETURN
      END


text:
Address(Eucalyptus Street) is a residential street in South Davis.

 Intersecting Streets

  San Marino Drive
  Boxelder Place
  Dogwood Place

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!==============================================================================!
  program Test_Metis_01
!------------------------------------------------------------------------------!
!   First demonstration program which calls Metis library.                     !
!------------------------------------------------------------------------------!
  implicit none
!----------------------------------[Locals]------------------------------------!
  integer              :: nels   = 2,  &  ! number of elements
                          nnds   = 6,  &  ! number of nodes
                          npel   = 4,  &  ! nodes per element
                          nparts = 2,  &  ! number of desired partitions
                          objval          ! returns something after METIS
  integer, allocatable :: eptr(:),   &
                          nodes(:),  &
                          epart(:),  &
                          npart(:)
  integer, pointer     :: vwgt   => null(),  &
                          vsize  => null()
  real, pointer        :: tpwgts => null()
  integer              :: mopts(41)       ! options passed to METIS

  integer, parameter :: METIS_OPTION_PTYPE     =  1
  integer, parameter :: METIS_OPTION_OBJTYPE   =  2
  integer, parameter :: METIS_OPTION_CTYPE     =  3
  integer, parameter :: METIS_OPTION_IPTYPE    =  4
  integer, parameter :: METIS_OPTION_RTYPE     =  5
  integer, parameter :: METIS_OPTION_DBGLVL    =  6
  integer, parameter :: METIS_OPTION_NITER     =  7
  integer, parameter :: METIS_OPTION_NCUTS     =  8
  integer, parameter :: METIS_OPTION_SEED      =  9
  integer, parameter :: METIS_OPTION_NO2HOP    = 10
  integer, parameter :: METIS_OPTION_MINCONN   = 11
  integer, parameter :: METIS_OPTION_CONTIG    = 12
  integer, parameter :: METIS_OPTION_COMPRESS  = 13
  integer, parameter :: METIS_OPTION_CCORDER   = 14
  integer, parameter :: METIS_OPTION_PFACTOR   = 15
  integer, parameter :: METIS_OPTION_NSEPS     = 16
  integer, parameter :: METIS_OPTION_UFACTOR   = 17
  integer, parameter :: METIS_OPTION_NUMBERING = 18
  integer, parameter :: METIS_OPTION_HELP      = 19
  integer, parameter :: METIS_OPTION_TPWGTS    = 20
  integer, parameter :: METIS_OPTION_NCOMMON   = 21
  integer, parameter :: METIS_OPTION_NOOUTPUT  = 22
  integer, parameter :: METIS_OPTION_BALANCE   = 23
  integer, parameter :: METIS_OPTION_GTYPE     = 24
  integer, parameter :: METIS_OPTION_UBVEC     = 25
!------------------------------------------------------------------------------!

  allocate(eptr(nels+1))
  allocate(nodes(nels*npel))
  allocate(epart(nels))
  allocate(npart(nnds))

  ! Initialize data to be sent to Metis
  eptr  = (/0, 4, 8/)
  nodes = (/1, 2, 3, 4, 2, 5, 6, 3/)  ! element 1 has nodes 1 2 3 4
                                      ! element 2 has nodes 2 5 6 3

  ! Set options to default values
  mopts = -1
  mopts(METIS_OPTION_NUMBERING) = 1
  mopts(METIS_OPTION_DBGLVL)    = 1

  ! Call Metis
  call METIS_PartMeshNodal(nels,    &
                           nnds,    &
                           eptr,    &
                           nodes,   &
                           vwgt,    &
                           vsize,   &
                           nparts,  &
                           tpwgts,  &
                           mopts,   &
                           objval,  &
                           epart,   &
                           npart)

  print *, '# Node partitioning: ', npart
  print *, '# Elem partitioning: ', epart

  end program

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE FACTT
!                                                                      *
!   Calculates the logs  of factorials required by the Racah coeffi-   *
!   cient routine DRACAH.                                              *
!                                                                      *
!   Written by N S Scott                    Last update: 15 Oct 1992   *
!                                                                      *
!***********************************************************************
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  10:47:26   2/14/04
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE FACTS_C
      IMPLICIT NONE
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I
      REAL(DOUBLE) :: X
!-----------------------------------------------
!
      GAM(1) = 1.0D00
      GAM(2) = 1.0D00
      X = 2.0D00
!
      DO I = 3, 30
         GAM(I) = GAM(I-1)*X
         X = X + 1.0D00
      END DO
!
      DO I = 1, 30
         GAM(I) = LOG(GAM(I))
      END DO
!
      X = 3.0D01
!
      DO I = 31, MFACT
         GAM(I) = GAM(I-1) + LOG(X)
         X = X + 1.0D00
      END DO
!
      RETURN
      END SUBROUTINE FACTT

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
	SUBROUTINE GDAINP  ( gdfile, gdatim, gpack, iret )
C************************************************************************
C* GDAINP								*
C*									*
C* This subroutine gets the input parameters for GDOMEG.		*
C*									*
C* GDAINP  ( GDFILE, GDATIM, GPACK, IRET )				*
C**									*
C* Log:									*
C* K. Brill/NMC		10/90						*
C8 M. desJardins/NMC	 1/92	GDOINP-->GDAINP				*
C************************************************************************
	CHARACTER*(*)	gdfile, gdatim, gpack
C------------------------------------------------------------------------
	CALL IP_STR  ( 'GDFILE',  gdfile, ier1 )
	CALL IP_STR  ( 'GDATTIM', gdatim, ier2 )
	CALL IP_STR  ( 'GPACK',   gpack,  ier3 )
	iret = ier1 + ier2 + ier3
	IF  ( iret .ne. 0 )  iret = -2
C*
	RETURN
	END


text:
! Anantua 2019 semianalytic jet model
! Check conversion to cgs units!

module fluid_model_rrjet

      use class_four_vector
      use interpolate, only: interp
      use kerr, only: kerr_metric, calc_u0
      use phys_constants, only: pi, c2, m
      implicit none

      namelist /fluiddata/ betaeconst, ximax

      real :: betaeconst, ximax

      interface init_rrjet
        module procedure init_rrjet
      end interface

      interface rrjet_vals
        module procedure rrjet_vals
      end interface

      contains

        subroutine rrjet_vals(x0,a,rho,p,b,u,bmag)
        type (four_Vector), intent(in), dimension(:) :: x0
        real, intent(in) :: a
        real(kind=8), dimension(size(x0)) :: done
        real, dimension(size(x0)) :: x2,x1,zm,zr,theta,s,z,logz,xi,fone
        real, dimension(size(x0)) :: omega, phi, phiprime, current, Bs, Bz, Bphi,Vs,Vz,Vphi,V0z,Vr,Vth,Br,Bth
        real, dimension(size(x0)) ::  u0,b0,vpl0,vrl0,vtl0,p0,rho0,bph0,bth0,rd,td,dzero,vr0,vth0,vph0,dummy
        real(8), dimension(size(x0),10) :: metric
        integer :: npts,nx1,nx2
        real, dimension(size(x0)), intent(out) :: rho,p,bmag
        type (four_Vector), intent(out), dimension(size(x0)) :: u,b

!        write(6,*) 'rrjet: ',size(x0)
        dzero=0d0; done=1d0; fone=1.;

        ! Get spherical coordinates
        npts=size(x0)
        theta=x0%data(3)
        zm=cos(theta)

        ! Jet solution is symmetric about xy plane
        ! ANDREW ?? -- is this ok?
        x2=acos(abs(zm)); theta=x2
        zr=x0%data(2)

        ! Get cylindrical coordinates
        s = zr*sin(theta)
        z = zr*cos(theta)
        xi = s*s / z
        where(z.eq.0)
          xi=1.e10
        endwhere

        ! Get magnetic flux, fluid line angular velocity, and current
        omega = 0.15*exp(-0.3*xi*xi)
        phi = tanh(0.3*xi)
        phiprime = 0.3*(1.- phi*phi)
        current = -2.*xi*omega*phiprime

        ! Get the B-field vector in cylindrical coordinates
        ! !AC - in dimensionless units, the flux through the horizon is 1
        ! !AC - we convert to physical units in  convert_fluidvars_rrjet in fluid.f90
        Bs = 1.02e4 * s * phiprime / (2.*pi*z*z)
        Bz = 1.02e4 * phiprime / (pi*z)
        !Bz = 1.02e4 * phiprime / (2.*pi*z)
        Bphi = 1.02e4 * current / (2*pi*s)
        Bphi = Bphi / s !AC -- put in coordinate basis

        ! Get the three-velocity in cylindrical coordinates
        logz = log10(z)
        V0z = -1.2 + 2.2*tanh(0.84*logz)
        Vz = V0z * exp(-0.001 * xi**4)
        Vs = s*Vz / (2.* z)
        Vphi = s*omega*(1.-Vz)
        Vphi = Vphi / s !AC -- put in coordinate basis

        ! Transform the velocity and B-field to spherical coordinates
        Br = s*Bs/zr + z*Bz/zr
        Bth = z*Bs/(zr*zr) - s*Bz/(zr*zr)
        Vr = s*Vs/zr + z*Vz/zr
        Vth = z*Vs/(zr*zr) - s*Vz/(zr*zr)

        ! Find the timelike u0 and convert V -> u
        metric = kerr_metric(zr,real(x0%data(3)),a)
        u0 = calc_u0(metric,dble(Vr),dble(Vth),dble(Vphi))

        u%data(1)=u0
        u%data(2)=Vr*u0
        u%data(3)=Vth*u0
        u%data(4)=Vphi*u0
        call assign_metric(u,transpose(kerr_metric(zr,real(x0%data(3)),a)))

        ! Find b0 and convert B -> b
        b%data(1)=dzero
        b%data(2)=Br
        b%data(3)=Bth
        b%data(4)=Bphi
        call assign_metric(b,transpose(kerr_metric(zr,real(x0%data(3)),a)))

        b0 = b*u

        b%data(1)=b0
        b%data(2)=(Br + b0*u%data(2))/u0
        b%data(3)=(Bth + b0*u%data(3))/u0
        b%data(4)=(Bphi + b0*u%data(4))/u0
        !b%data(2)=u0*(Br + b0*u%data(2))
        !b%data(3)=u0*(Bth + b0*u%data(3))
        !b%data(4)=u0*(Bphi + b0*u%data(4))

        ! Protect azimuthal velocities at poles
        ! Compute magnitude of interpolated b-field
        ! and force b^2 > 0 (need to look into numerical issues here):
        call assign_metric(b,transpose(kerr_metric(zr,real(x0%data(3)),a)))
        bmag=b*b;
        bmag=merge(bmag,dzero,bmag.ge.0d0)
        bmag=sqrt(bmag)
!        write(6,*) 'maxr: ',maxval(uniqr),maxval(zr)
!        write(6,*) 'bmag: ',maxval(bmag),minval(bmag)
!        write(6,*) 'bmag: ',bmag

        ! Find the nonthermal electron pressure (betae * magnetic pressure)
        p = 0.5 * betaeconst * (bmag*bmag) !/ (8*pi)
        rho = p ! Convert pressure to number density in  convert_fluidvars_rrjet using emis params

        ! Zero everything where xi>ximax and z<2
        where((xi.gt.ximax).or.(abs(z).le.2d0))
           !u%data(1) = dzero;
           !u%data(2) = dzero;
           !u%data(3) = dzero;
           !u%data(4) = dzero;
           !b%data(1) = dzero;
           !b%data(2) = dzero;
           !b%data(3) = dzero;
           !b%data(4) = dzero;
           rho = dzero;
           bmag=dzero;
           p=dzero;
        endwhere

        ! !AC Cut off emission from below the eq. plane:
        rho=merge(rho,rho*0.,zm.ge.0)
        p = merge(p,p*0.,zm.ge.0)
        end subroutine rrjet_vals

        subroutine init_rrjet(betaeconst0, ximax0)
        !character(len=20), intent(in) :: ifile
        !open(unit=8,file=ifile,form='formatted',status='old')
          !read(8,nml=fluiddata)
          real(8), intent(in) :: betaeconst0,  ximax0
          betaeconst = betaeconst0
          ximax=ximax0
        write(6,*) 'read: ',betaeconst,ximax
        close(unit=8)
        end subroutine init_rrjet

      end module fluid_model_rrjet

max_stars_repo_name: RAnantua/grtrans

max_issues_repo_name: RAnantua/grtrans

max_forks_repo_name: RAnantua/grtrans


text:
!*******************************************************************************
!  File ipsl_cdf.f
!  Contains the module ipsl_cdf
!    Module for defining variables and writing netCDF files, and reading
!    netCDF files with the derived types ipsl_desc and ipsl_data
!    (from the ipsl_T module).
!
!    Information about the  EZcdf module is at:
!       http://w3.pppl.gov/NTCC/EZcdf/
!
!-------------------------------------------------------------------------------
!   DEPENDENCIES
!-------------------------------------------------------------------------------
!
!    This module uses the following modules:
!       stel_kinds
!       stel_constants
!       ipsl_T
!       bsc
!       bsc_cdf
!       ezcdf
!       v3_utilities
!
!-------------------------------------------------------------------------------
!   CHANGE HISTORY
!-------------------------------------------------------------------------------
!
!  See Section X at the end of the module.
!
!-------------------------------------------------------------------------------
!   USAGE
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!   COMMENTS
!
! 1) All cdf_define calls must be completed before the first cdf_write call.
!-------------------------------------------------------------------------------
!
!*******************************************************************************

!*******************************************************************************
!  MODULE ipsl_cdf
!
! SECTION I. VARIABLE DECLARATIONS
! SECTION II. INTERFACE BLOCKS
! SECTION III. DEFINE SUBROUTINES
! SECTION IV. WRITE SUBROUTINES
! SECTION V. READ SUBROUTINES
! SECTION VI. AUXILIARY FUNCTIONS AND SUBROUTINES

! SECTION X.  COMMENTS FOR DIFFERENT REVISIONS
!*******************************************************************************

      MODULE ipsl_cdf

!*******************************************************************************
! SECTION I. VARIABLE DECLARATIONS
!*******************************************************************************

!-------------------------------------------------------------------------------
!  Type declarations - lengths of reals, integers, and complexes.
!  Frequently used mathematical constants, lots of extra precision.
!-------------------------------------------------------------------------------
      USE stel_kinds, only : rprec, iprec, cprec
      USE stel_constants, only : pi, twopi, one, zero

!-------------------------------------------------------------------------------
!  Modules to USE
!-------------------------------------------------------------------------------
      USE ipsl_T
      USE ip_beamline
      USE intpol_cdf
      USE ezcdf
      USE v3_utilities

!-------------------------------------------------------------------------------
!  Implicit None comes after USE statements, before other declarations
!-------------------------------------------------------------------------------
      IMPLICIT NONE

!-------------------------------------------------------------------------------
!  Make type declarations and constants Private, so there are no conflicts.
!-------------------------------------------------------------------------------
      PRIVATE rprec, iprec, cprec, pi, twopi, one, zero

!-------------------------------------------------------------------------------
!  Variable Names for netCDF. Make them Private.
!-------------------------------------------------------------------------------

      CHARACTER (LEN=*), PRIVATE, PARAMETER ::                                 &
     &  vn_s_name = 'ipsl_desc_s_name',                                        &
     &  vn_l_name = 'ipsl_desc_l_name',                                        &
     &  vn_ip_sname = 'ipsl_desc_ip_sname',                                        &
     &  vn_ip_lname = 'ipsl_desc_ip_lname',                                        &
     &  vn_units = 'ipsl_desc_units',                                          &
     &  vn_sigma_default = 'ipsl_desc_sigma_default',                          &
     &  vn_l_ipbeam_def = 'ipsl_desc_l_ipbeam_def',                            &
     &  vn_ipsl_type = 'ipsl_desc_ipsl_type'

      CHARACTER (LEN=64), PRIVATE ::                                           &
     &  vn_s_name_use,                                                         &
     &  vn_l_name_use,                                                         &
     &  vn_ip_sname_use,                                                         &
     &  vn_ip_lname_use,                                                         &
     &  vn_units_use,                                                          &
     &  vn_sigma_default_use,                                                  &
     &  vn_l_ipbeam_def_use,                                                   &
     &  vn_ipsl_type_use

!-------------------------------------------------------------------------------
!  Lengths of Character Variables
!-------------------------------------------------------------------------------
      INTEGER(iprec), PARAMETER, PRIVATE :: type_len=10
      INTEGER(iprec), PARAMETER, PRIVATE :: sn_len=30
      INTEGER(iprec), PARAMETER, PRIVATE :: ln_len=80
      INTEGER(iprec), PARAMETER, PRIVATE :: units_len=30

!-------------------------------------------------------------------------------
!  Generic Define
!-------------------------------------------------------------------------------
      INTERFACE ipsl_cdf_define
         MODULE PROCEDURE ipsl_cdf_define_desc
      END INTERFACE

!-------------------------------------------------------------------------------
!  Generic Write
!-------------------------------------------------------------------------------
      INTERFACE ipsl_cdf_write
         MODULE PROCEDURE ipsl_cdf_write_desc
       END INTERFACE

!-------------------------------------------------------------------------------
!  Generic Read
!-------------------------------------------------------------------------------
      INTERFACE ipsl_cdf_read
         MODULE PROCEDURE ipsl_cdf_read_desc
       END INTERFACE
!-------------------------------------------------------------------------------

      CONTAINS
!*******************************************************************************
! SECTION III. DEFINE SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
      SUBROUTINE ipsl_cdf_define_desc(this,iou,prefix)
!  Subroutine to do the appropriate netCDF definition calls for a ipsl_desc
!
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      TYPE (ipsl_desc), INTENT (in)            :: this
      INTEGER, INTENT(in)                      :: iou
      CHARACTER (len=*), INTENT(in), OPTIONAL  :: prefix

!  this        ipsl_desc - this is what gets defined and written.
!  iou         i/o unit number of the netCDF file.
!  prefix      character -  prefixed to variable names, so that netCDF
!                  doesn't have problems with repeated, identical names
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'ipsl_cdf_define_desc: '
      CHARACTER(len=32) :: prefix_use

!-----------------------------------------------
!  Start of Executable Code
!-----------------------------------------------

! Define the prefix to actually use
      IF (PRESENT(prefix)) THEN
         prefix_use = TRIM(ADJUSTL(prefix))
      ELSE
         prefix_use = ' '
      ENDIF

! Define all vn_--_use variable names
      CALL ipsl_cdf_defvn_desc(prefix_use)

! Scalar Components
      CALL cdf_define(iou, TRIM(vn_s_name_use), this % s_name)
      CALL cdf_define(iou, TRIM(vn_l_name_use), this % l_name)
      CALL cdf_define(iou, TRIM(vn_ip_sname_use), this % ip_sname)
      CALL cdf_define(iou, TRIM(vn_ip_lname_use), this % ip_lname)
      CALL cdf_define(iou, TRIM(vn_units_use), this % units)
      CALL cdf_define(iou, TRIM(vn_sigma_default_use),                         &
     &   this % sigma_default)
      CALL cdf_define(iou, TRIM(vn_l_ipbeam_def_use),                          &
     &      this % l_ipbeam_def)
      CALL cdf_define(iou, TRIM(vn_ipsl_type_use),                             &
     &      this % ipsl_type)

!  define ip_beam
      CALL intpol_cdf_define(this % ipbeam,iou,prefix)

      RETURN

      END SUBROUTINE ipsl_cdf_define_desc

!*******************************************************************************
! SECTION III. WRITE SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
      SUBROUTINE ipsl_cdf_write_desc(this,iou,prefix)
!  Subroutine to do the appropriate netCDF write calls for a ipsl_desc
!
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      TYPE (ipsl_desc), INTENT (in)        :: this
      INTEGER, INTENT (in)                       :: iou
      CHARACTER (len=*), INTENT (in), OPTIONAL   :: prefix

!  this        ipsl_desc - this is what gets written.
!  iou         i/o unit number of the netCDF file.
!  prefix      character -  prefixed to variable names, so that netCDF
!                  doesn't have problems with repeated, identical names
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'ipsl_cdf_write_desc: '
      CHARACTER(len=32) :: prefix_use

!-----------------------------------------------
!  Start of Executable Code
!-----------------------------------------------
! Define the prefix to actually use
      IF (PRESENT(prefix)) THEN
         prefix_use = TRIM(ADJUSTL(prefix))
      ELSE
         prefix_use = ' '
      ENDIF

! Define all vn_--_use variable names
      CALL ipsl_cdf_defvn_desc(prefix_use)

! Scalar Components
      CALL cdf_write(iou, TRIM(vn_s_name_use), this % s_name)
      CALL cdf_write(iou, TRIM(vn_l_name_use), this % l_name)
      CALL cdf_write(iou, TRIM(vn_ip_sname_use), this % ip_sname)
      CALL cdf_write(iou, TRIM(vn_ip_lname_use), this % ip_lname)
      CALL cdf_write(iou, TRIM(vn_units_use), this % units)
      CALL cdf_write(iou, TRIM(vn_sigma_default_use),                         &
     &   this % sigma_default)
      CALL cdf_write(iou, TRIM(vn_l_ipbeam_def_use),                          &
     &      this % l_ipbeam_def)
      CALL cdf_write(iou, TRIM(vn_ipsl_type_use),                             &
     &      this % ipsl_type)

!  ip_beam part
      CALL intpol_cdf_write(this % ipbeam,iou,prefix)

      RETURN

      END SUBROUTINE ipsl_cdf_write_desc

!*******************************************************************************
! SECTION IV. READ SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
      SUBROUTINE ipsl_cdf_read_desc(this,iou,prefix)
!  Subroutine to do the appropriate netCDF read calls for a ipsl_desc
!
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      TYPE (ipsl_desc), INTENT (inout)        :: this
      INTEGER, INTENT (in)                       :: iou
      CHARACTER (len=*), INTENT (in), OPTIONAL   :: prefix

!  this        ipsl_desc - this is what gets defined and written.
!  iou         i/o unit number of the netCDF file.
!  prefix      character -  prefixed to variable names, so that netCDF
!                  doesn't have problems with repeated, identical names
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'ipsl_cdf_read_desc: '
      CHARACTER(len=32) :: prefix_use
      INTEGER(iprec), DIMENSION(3) :: dimlens

      CHARACTER (len=sn_len)    :: s_name
      CHARACTER (len=ln_len)    :: l_name
      CHARACTER (len=sn_len)    :: ip_sname
      CHARACTER (len=ln_len)    :: ip_lname
      CHARACTER (len=units_len) :: units
      CHARACTER (len=30)        :: ipsl_type
      LOGICAL                   :: l_ipbeam_def
      REAL(rprec)               :: sigma_default
      TYPE (ip_beam)            :: ipbeam

!-----------------------------------------------
!  Start of Executable Code
!-----------------------------------------------

! Define the prefix to actually use
      IF (PRESENT(prefix)) THEN
         prefix_use = TRIM(ADJUSTL(prefix))
      ELSE
         prefix_use = ' '
      ENDIF

! Define all vn_--_use variable names
      CALL ipsl_cdf_defvn_desc(prefix_use)

! Read Scalar pieces
! Note: Read in to variables local to this subroutine.
      CALL cdf_read(iou, TRIM(vn_s_name_use), s_name)
      CALL cdf_read(iou, TRIM(vn_l_name_use), l_name)
      CALL cdf_read(iou, TRIM(vn_ip_sname_use), ip_sname)
      CALL cdf_read(iou, TRIM(vn_ip_lname_use), ip_lname)
      CALL cdf_read(iou, TRIM(vn_units_use), units)
      CALL cdf_read(iou, TRIM(vn_sigma_default_use), sigma_default)
      CALL cdf_read(iou, TRIM(vn_l_ipbeam_def_use),l_ipbeam_def)
      CALL cdf_read(iou, TRIM(vn_ipsl_type_use),ipsl_type)

!  Read Derived Type
      CALL intpol_cdf_read(ipbeam,iou,prefix)

! Create the ipsl_desc, this
      CALL ipsl_desc_construct(this,s_name,l_name,ip_sname,ip_lname,           &
     &   units, sigma_default,ipsl_type,ipbeam)

!  Destroy the local int/pol beamline ipbeam  to avoid memory leakage
      CALL ip_beam_destroy(ipbeam)

      RETURN

      END SUBROUTINE ipsl_cdf_read_desc
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------

!*******************************************************************************
! SECTION IV. AUXILIARY FUNCTIONS AND SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
      SUBROUTINE ipsl_cdf_defvn_desc(prefix_use)
!  Subroutine to do define the character variable names for a ipsl_desc,
!  using the prefix. All the vn_ variables are module variables, and so do not
!   need to be declared here

!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      CHARACTER (len=*), INTENT (in)   :: prefix_use

!  prefix_use      character -  prefixed to variable names, so that netCDF
!                  doesn't have problems with repeated, identical names

!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'ipsl_cdf_defvn_desc: '

!-----------------------------------------------
!  Start of Executable Code
!-----------------------------------------------

! Define all variable names
      vn_s_name_use = ipsl_cdf_mknam(prefix_use,vn_s_name)
      vn_l_name_use = ipsl_cdf_mknam(prefix_use,vn_l_name)
      vn_ip_sname_use = ipsl_cdf_mknam(prefix_use,vn_ip_sname)
      vn_ip_lname_use = ipsl_cdf_mknam(prefix_use,vn_ip_lname)
      vn_units_use = ipsl_cdf_mknam(prefix_use,vn_units)
      vn_sigma_default_use = ipsl_cdf_mknam(prefix_use,                        &
     &   vn_sigma_default)
      vn_l_ipbeam_def_use = ipsl_cdf_mknam(prefix_use,                         &
     &   vn_l_ipbeam_def)
      vn_ipsl_type_use = ipsl_cdf_mknam(prefix_use,                            &
     &   vn_ipsl_type)

      RETURN

      END SUBROUTINE ipsl_cdf_defvn_desc

!========================================================
      FUNCTION ipsl_cdf_mknam(c1,c2)
! A simple function to help in the generation of names

!-----------------------------------------------
!   F u n c t i o n   N a m e
!-----------------------------------------------
      CHARACTER(LEN=64) ipsl_cdf_mknam

!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      CHARACTER(LEN=*), INTENT (in) :: c1,c2

!-----------------------------------------------
!  Start of Executable Code
!-----------------------------------------------
      IF (LEN_TRIM(c1) .eq. 0) THEN
         ipsl_cdf_mknam = TRIM(c2)
      ELSE
         ipsl_cdf_mknam = ADJUSTL(TRIM(c1) // '_' // TRIM(c2))
      ENDIF

      RETURN

      END FUNCTION ipsl_cdf_mknam

!-----------------------------------------------
!-----------------------------------------------
!*******************************************************************************
! SECTION X.  COMMENTS FOR DIFFERENT REVISIONS
!*******************************************************************************
!
! JMS 2007-06-14. First version of ipsl_cdf.
!    Based on diagnostic_cdf
!

      END MODULE ipsl_cdf