! @Copyright 2007 Kristjan Haule
!

SUBROUTINE cmp_overlap(projector,Olapm, SOlapm, Qcomplex, nsymop, csort, iorbital, cix_orb, nindo, cixdim, iSx, noccur, cfX, crotloc_x_rotij, maxucase, maxdim2, norbitals, pr_proc, SIMPLE, Qsymmetrize)
  USE com_mpi,  ONLY: myrank, master, vectors, nvector, vector_para, fvectors, AllReduce_MPI, Qprint
  USE com,      ONLY: MINWAV, MAXWAV, iso, emin, emax
  USE abc,      ONLY: KX, KY, KZ, BK3, E, A, ALM, ALML, aK
  USE param,    ONLY: nmat , LMAX2, iblock, nrf, NRAD, LOMAX, nloat, nemin0, nemax0, max_nl
  USE w_atpar,  ONLY: w_FJ, w_DFJ, w_jatom, w_alo, w_nlo, w_nlov, w_nlon, w_ilo, w_lapw, w_ri_mat, w_P, w_DP, nnlo
  USE structure,ONLY: RMT, VOL, pos, tau, nat, rotij, tauij, iord, iz, BR1, rot_spin_quantization
  USE case,     ONLY: iatom, isort, nl, ll, crotloc, shft, maxdim, ncix, natom, Sigind, csize, maxsize, cix, ifirst
  USE sym2,     ONLY: tmat, idet, iz_cartesian
  USE kpts,     ONLY: mweight, tweight, numkpt
  USE p_project,ONLY: P_rfi, rix_mat, n_al_ucase, max_lcase, al_ucase, l_al_ucase, j_al_ucase, kNri, rix, dri, al_interstitial
  USE matpar,   only: atpar, alo, nlo, nlov, nlon, ilo, lapw, RI_MAT, P, DP
  IMPLICIT NONE
  INTEGER, intent(in)     :: projector
  COMPLEX*16, intent(out) :: Olapm(maxdim,maxdim,ncix), SOlapm(maxdim,maxdim,ncix)
  LOGICAL, intent(in)     :: Qcomplex
  INTEGER, intent(in)     :: nsymop, csort(nat), pr_proc
  INTEGER, intent(in)     :: iorbital(natom,lmax2+1), cix_orb(norbitals), nindo(norbitals)
  INTEGER, intent(in)     :: cixdim(ncix), iSx(maxdim2, norbitals), noccur(maxsize,ncix)
  COMPLEX*16, intent(in)  :: cfX(maxdim2,maxdim2,norbitals,norbitals)
  INTEGER, intent(in)     :: maxucase, maxdim2, norbitals
  REAL*8,  intent(in)     :: crotloc_x_rotij(3,3,max_nl,natom)
  LOGICAL, intent(in)     :: SIMPLE, Qsymmetrize
  ! interfaces
  interface
     REAL*8 Function Romb(y,N,dx)
       IMPLICIT NONE
       REAL*8, intent(in) :: y(N)
       REAL*8, intent(in) :: dx
       INTEGER, intent(in):: N
     end Function Romb
     REAL*8 FUNCTION detx(a)
       IMPLICIT NONE
       REAL*8, intent(in) :: a(3,3)
     end FUNCTION detx
  end interface
  ! locals
  real*8,     allocatable :: a_real(:)
  complex*16, allocatable :: DMFTU(:,:,:), olp(:,:)
  complex*16, allocatable :: Olapmk(:,:,:)
  COMPLEX*16, allocatable :: URx(:,:), tmp(:,:)
  !COMPLEX*16, allocatable :: Uu(:,:), Ud(:,:)
  COMPLEX*16, allocatable :: a_interstitial(:,:,:), h_interstitial(:,:,:)
  REAL*8,     allocatable  :: phi_jl(:,:)
  REAL*8,     allocatable :: aKR(:), jlr(:), jlp(:)
  COMPLEX*16, allocatable :: olocf(:,:), work(:),tmp1(:,:), solap(:,:)
  REAL*8,     allocatable :: ws(:), rwork(:)
  INTEGER,    allocatable :: iwork(:)
  !
  INTEGER, allocatable :: cind(:), cini(:)
  INTEGER :: cixdms
  !
  INTEGER      :: ikp, iikp, iks, nkp, ivector, N, NE, NEMIN, NEMAX, nbands, NUM, isym, icase, lcase, lfirst
  INTEGER      :: latom, jatom, iucase, i, IND_YL, ibb, lda, ldb, ldc, l1
  INTEGER      :: nind, idt, icix, iorb, m1, lms1, ind1, ind2, is1, iorb1, iorb2, nind1, nind2
  INTEGER      :: cixdm, ip, iq, it, is, itape, L, M, ii, i3, N1, num1, irf1, irf2
  logical      :: Tcompute
  character*10 :: KNAME
  real*8       :: S, T, Z, exxx, FAC, PI, TWOPI, ARG123, ARG2, ARGT, EMIST, rx
  complex*16   :: h_yl(2*LMAX2+1,iblock)
  complex*16   :: h_alyl(2*LMAX2+1,iblock,2), h_blyl(2*LMAX2+1,iblock,2)
  complex*16   :: YL((LMAX2+1)*(LMAX2+1))
  complex*16   :: PHSHEL, CFAC, IMAG, cc, corc, sm1, sm2, ff
  real*8       :: Olapc(maxsize,ncix), ARGT2, tmp3(3,3)
!  COMPLEX*16   :: tuu, tud, tdu, tdd
  CHARACTER*200:: FNAME
  !
  real*8       :: BK(3), BKROT(3), BKROT2(3), BKRLOC(3)
  REAL*8       :: Kn(3), crotloc_x_BR1(3,3)
  INTEGER      :: lwork, lrwork, liwork, info, iind, ir, Nri, isize, isize_
  LOGICAL      :: nonzero_interst
  REAL*8       :: Det, phi1, the1, psi1
  COMPLEX*16   :: Rispin(2,2)
  INTEGER      :: Nstar
  REAL*8       :: k_star(3,iord), Trans3(3,3)
  INTEGER      :: gind(iord), nsymop2, igi, j
  !------------------------------------------------------------------

  DATA IMAG/(0.0D0,1.0D0)/

  PI=ACOS(-1.0D0)
  TWOPI=2.D0*PI

  if (Qprint) WRITE(6,'(A,I3,2x,A,I3)') 'pr_proc=', pr_proc, 'tot-k=', numkpt

  ALLOCATE( a_real(nmat) ) !---  for eigenvectors -------------!
  allocate( olp(maxdim2,maxdim2) )
  allocate( Olapmk(maxdim,maxdim,ncix) )
  Olapm=0

  iikp=0
  DO ivector=1,nvector

     DO is=1,iso    !------ over up/dn ---------------------!
        itape=8+is
        if (vector_para) then
           FNAME = fvectors(ivector,is)
           open(itape,FILE=FNAME,STATUS='old',FORM='unformatted')
        else
           !------------ vector files need to be read from the beginning -------------------!
           rewind(itape) !--- both vector files: 9 and 10 rewind -------------------------!
        endif
        DO I=1,NAT
           READ(itape) EMIST !---- At the beginninge we have linearization energies --!
           READ(itape) EMIST !---- We just skip them ---------------------------------!
        ENDDO
     END DO

     if (vector_para) then
        nkp = vectors(ivector,2)
     else
        nkp= numkpt
     endif

     DO iks=1,nkp   !------ Over all irreducible k-points ------!
        if (vector_para) then
           ikp = vectors(ivector,3)+iks  ! successive index in k-point table from case.klist
           iikp = iikp+1                 ! successive index in k-point table on this processor
        else
           ikp = iks                     ! successive index in k-point table from case.klist
           !--- We need to go over all k-points even though we will compute only some of them on this processor.
           !--- This is because we need to read vector file sequentially.
           iikp = ikp-myrank*pr_proc            ! The index of the point to compute. If negative, do not compute!
        endif

        Tcompute=.FALSE.
        if (iikp.gt.0) Tcompute=.TRUE.       ! If Tcompute is true, the point needs to be computed.
        if (iikp.GT.pr_proc) EXIT            ! Processor finished. The rest of the points will be taken care of by the other processors.

        if (Tcompute .and. abs(projector).eq.5) allocate( phi_jl(nmat, n_al_ucase) )

        !------- reading from vector for both spins -----------!
        DO is=1,iso    !------ over up/dn ---------------------!
           itape=8+is
           READ(itape,END=998) S,T,Z,KNAME,N,NE  !--- Here we can jupm out of loop 4 and stop at 998 -----------------------------------!
           IF(N.GT.MAXWAV) MAXWAV=N
           IF(N.LT.MINWAV) MINWAV=N
           READ(itape) (KX(I),KY(I),KZ(I),I=1,N) !--- Reads all reciprocal vectors -----------------------------------------------------!
           NEMIN=1
           NEMAX=0
           NUM=0
           DO WHILE (NUM.NE.NE)
              READ(itape) NUM,exxx                !----- eigenvalues read -----!
              E(NUM)=exxx
              if (.not.Qcomplex) then
                 READ(itape) (A_real(I),I=1,N)    !----- eigenvector read -----!
                 a(:,num,is) = a_real(:)
              else
                 READ(itape) (A(I,NUM,is),I=1,N)  !--- eigenvector complex of size n (where n is number of reciprocal vectors) --!
              end if
              if (abs(projector).LT.4) then
                 IF(E(NUM).LT.EMIN) NEMIN=NEMIN+1
                 IF(E(NUM).LT.EMAX) NEMAX=NEMAX+1
              !else
              !   nemin=nemin0
              !   nemax=nemax0
              endif
           ENDDO

           if (abs(projector).GE.4) then
              nemin=nemin0
              nemax=min(nemax0,NE)
           endif

           IF (Tcompute .and. is.eq.1) THEN
              isize=N-nnlo
              DO I=1,N                                !----- Over reciprocal vectors --------------------------------------------------------!
                 BK3(1,I)=(S+KX(I))                   !----  Creates K+k where K is reciprocal vector and k is irreducible k-point ----------!
                 BK3(2,I)=(T+KY(I))
                 BK3(3,I)=(Z+KZ(I))
                 Kn(:) = matmul(BR1,BK3(:,I))
                 aK(I) = sqrt(Kn(1)**2+Kn(2)**2+Kn(3)**2)
              ENDDO
              DO iucase=1,maxucase
                 jatom = w_jatom(iucase)
                 !--- Computes Bessel functions for all needed l's and atoms. --------------------------------!
                 !---  For optimization purposes, done only for irreducible k-points -------------------------!
                 !---  output: FJ, DFJ / Spherical Bessel functions and derivative, uses common block BR1/general ---!
                 CALL HARMON(isize,aK(:isize),LMAX2,w_FJ(:,:isize,iucase),w_DFJ(:,:isize,iucase),RMT(jatom))
              ENDDO
              ! New
              if (abs(projector).EQ.5) then
                 Nri=2**kNri+1    ! Number of radial points in the interstitials
                 phi_jl(:,:)=0
                 allocate( aKR(Nri), jlr(Nri), jlp(Nri) )
                 DO iind=1,n_al_ucase
                    if (abs(dri(iind)).gt.1e-10) then
                       l     = l_al_ucase(iind)
                       jatom = j_al_ucase(iind)
                       DO i=1,isize         ! over all reciprocal vectors K
                          rx = RMT(jatom)
                          do ir=1,Nri
                             aKR(ir)=rx*aK(i)    ! |k+K|*r
                             rix(ir)=rx          !  r
                             rx = rx + dri(iind)
                          enddo
                          CALL sphbes2(l,Nri,aKR,jlr)  ! spherical bessel : j_l(|k+K|*r)
                          jlr(:) = jlr(:)*rix(:)       ! r*j_l . We need to do that, because the integral is Int[ (phi(r)/r)*j_l(r)*r^2]=Int[phi(r)*j_l(r)*r]
                          jlp(:) = jlr(:)*P_rfi(:,iind)
                          phi_jl(i,iind) = romb(jlp, Nri, dri(iind))  ! Integration over r on rix mesh
                       ENDDO
                    endif
                 ENDDO
                 deallocate( aKR, jlr, jlp )
              endif
              ! New
           ENDIF
        ENDDO  !
        if (.not.Tcompute) CYCLE  ! This k-points was read, but will not be computed on this processor
        !---  Finished reading eigenvectors and eigenvalues for both spins -----!
        !---   E(iband) contains eigenvalues   ---------------------------------!
        !---   A(:,iband,is) contains eigenvectors -----------------------------!
        nbands = nemax-nemin+1
        allocate( DMFTU(nbands,maxdim2,norbitals) )
        allocate( URx(nbands,maxdim2), tmp(nbands,maxdim2) )

        if (abs(projector).eq.5) then
           allocate( h_interstitial(2*LMAX2+1,iblock,2) )
           allocate( a_interstitial(2*LMAX2+1,nbands,iso), al_interstitial(2*LMAX2+1,nbands,iso,max_lcase) )
        endif

        if (.False.) then
           ! Finds all star members and corresponding group operations
           ! This is switched off, as there are some cases of selection of G vectors, which are not completely symmetric
           ! Taking only the star members can then lead to less symmetry as expected.
           BK(1)=S
           BK(2)=T
           BK(3)=Z
           call STERN(BK, Nstar, k_star, gind, iz, tau, iord, Qcomplex)
           nsymop2 = min(nsymop,Nstar)
        else
           nsymop2=nsymop
           do igi=1,iord
              gind(igi)=igi
           enddo
        endif
        !print *, 'tweight=', tweight, 'nsymop=', nsymop2, 'mweight=', mweight(ikp)

        !----------- sum over all k-points in the star of the irreducible k-point ---------------------------------!
        DO igi=1,nsymop2  !-- Over all symmetry operarations -> all k-points in the star of irreducible k point --!
           isym = gind(igi)
           DMFTU=0
           DO icase=1,natom  !--------------- over all atoms requested in the input ------------------------!
              latom = iatom(icase)   ! The succesive number of atom (all atoms counted)
              jatom = isort(latom)   ! The sort of the atom  ( == w_jatom(iucase) )
              lfirst = ifirst(latom)
              iucase = csort(jatom)  ! The renumbert sorts, such that the required atoms from the input give continuous index
              !----  Setting all values of common-blocks which depend on atom and were computed by atpar ---!
              alo(:,:,:,:) = w_alo(:,:,:,:,iucase)
              nlo = w_nlo(iucase)
              nlov = w_nlov(iucase)
              nlon = w_nlon(iucase)
              ilo(:) = w_ilo(:,iucase)
              lapw(:) = w_lapw(:,iucase)
              ri_mat(:,:,:,:) = w_ri_mat(:,:,:,:,iucase)
              P(:,:,:) = w_p(:,:,:,iucase)
              DP(:,:,:) = w_dp(:,:,:,iucase)
              if ((nlo+nlon+nlov).NE.nnlo) then
                 WRITE(6,*) 'ERROR: nlo+nlon+nlov should be equal to nnlo but is not', nlo+nlon+nlov, nnlo
                 STOP
              endif
              isize = N-(nlo+nlon+nlov)

              if (abs(projector).eq.5) al_interstitial(:,:,:,:)=0.d0

              FAC=4.0D0*PI*RMT(jatom)**2/SQRT(VOL)
              do lcase=1,nl(icase) !----------- loop over L(jatom) requested in the ionput ---------------!
                 l=ll(icase,lcase) !------ current L --!

                 if (iso.eq.2) then
                    !!  local_axis_defined_by_locrot  <- local_axis_of_equivalent_atom <- group_operation_symmetry <- from_spin_quantization_to_global_cartesian
                    !!* Trans3 = crotloc(:,:,icase) * rotij_cartesian * iz_cartesian(:,:,isym) * rot_spin_quantization
                    tmp3 = matmul(iz_cartesian(:,:,isym), rot_spin_quantization)
                    Trans3 = matmul(crotloc_x_rotij(:,:,lcase,icase),tmp3)
                    Det = detx(Trans3)
                    Trans3 = transpose(Trans3*Det)
                    CALL Angles_zxz(phi1,the1,psi1, Trans3 )
                    CALL Spin_Rotation(Rispin,phi1,the1,psi1)
                 endif
                 crotloc_x_BR1(:,:) = matmul( crotloc(:,:,lcase,icase),BR1 )

                 ALM = 0.0         !------  ALM(m,band,nrf,is) will hold product of eigenvectors and a/b expansion coefficients --!
                 nonzero_interst=.false.
                 if (abs(projector).eq.5) then
                    a_interstitial=0
                    iind=al_ucase(icase,lcase)
                    if (abs(dri(iind)).gt.1e-10) nonzero_interst=.true.
                 endif
                 !--------- blocks are for efficiency. Matrix is multiplied in block form. This must be important in the past, while modern BLAS should do that better. I think it is obsolete.
                 DO ii=1,isize,iblock !------ iblock is 128 for 32-bit system -------!
                    !-------- nlo-number of local orbitals -----!
                    i3=0
                    do i=ii,min(ii+iblock-1,isize)  ! 121
                       !---------  rotates ylm(k+K) to ylm(k'+K) where k' is in star of irreducible k. ------------!
                       i3=i3+1
                       ! BKROT = R_a.(k+K) transforms to the reducible k-point
                       BKROT = matmul(TMAT(:,:,isym), BK3(:,I))
                       ! BKROT2 = R_n.R_a.(k+K), transformation from the first atom to an equivalent atom
                       BKROT2 = matmul(rotij(:,:,latom), BKROT)
                       !---- BR1 transforms integer reciprocal lattice vectors, as given in the VECTORLIST of LAPW1, into cartesian system ----!
                       ! BKROT3 = R_n.R_a.(k+K), but in cartesian coordinate system
                       !BKROT3 = matmul(BR1,BKROT2)
                       !!---- BKRLOC = crotloc.R_n.R_a.(k+K),  rotates according to the user specified local coordinate system.
                       !BKRLOC = matmul(crotloc(:,:,icase), BKROT3)
                       BKRLOC = matmul(crotloc_x_BR1, BKROT2)
                       !---- YLM = Y_{L}(Rotloc.R_g.(k+K))
                       CALL YLM (BKRLOC,LMAX2,YL)  !
                       ! ARG123 = (R_g.(k+K)) *  R(iatom) * 2pi
                       ARG123 = dot_product( BKROT2, POS(:,lfirst) )*TWOPI
                       ! ARGT = (k+K)*tau(isym) * 2pi
                       ARGT = dot_product(BK3(:3,i), TAU(:3,isym))*TWOPI
                       ! ARGT2 = (R_a.(k+K)).tau_n * 2pi
                       ARGT2 = dot_product(BKROT, tauij(:3,latom))*TWOPI
                       ! ARG2 = (R_a.(k+K)) *  shft * 2pi
                       ARG2  = dot_product(BKROT, shft(:3,latom) )*TWOPI
                       ! PHSEHL = e^{I*2pi*( (R_a.(k+K))*tau_n + (K+k)*tau(isym) + (R_n.R_a.(k+K)*R(first)))}
                       PHSHEL=EXP(IMAG*(ARG123+ARG2+ARGT+ARGT2))
                       DO  M=1,2*L+1
                          IND_YL=M+L*L
                          h_yl(M,i3)=conjg(yl(ind_yl))*phshel !----- h_yl is rotated yl when k is rotated to k' -----!
                       END DO
                    enddo

                    DO is=1,iso  !--- over both spins
                       i3=0
                       do i=ii,min(ii+iblock-1,isize)
                          i3=i3+1
                          DO M=1,2*L+1
                             if (lapw(l)) then
                                h_ALYL(m,i3,is)=(w_DFJ(L,I,iucase)*P(l,is,2)-w_FJ(L,I,iucase)*DP(l,is,2))* h_yl(M,i3)!*ph_spin(is) ! derivatives of bessel functions and spheric harmonics
                                h_BLYL(m,i3,is)=(w_FJ(L,I,iucase)*DP(l,is,1)-w_DFJ(L,I,iucase)*P(l,is,1))* h_yl(M,i3)!*ph_spin(is)
                             else
                                h_ALYL(m,i3,is)=w_FJ(L,I,iucase)/P(l,is,1)/RMT(jatom)**2*h_yl(M,i3)!*ph_spin(is)
                                h_BLYL(m,i3,is)=(0.d0,0.d0)
                             end if
                          END DO
                       enddo
                       !!!  New
                       if (abs(projector).eq.5 .and. nonzero_interst) then
                          do i=ii,min(ii+iblock-1,isize)
                             DO M=1,2*L+1
                                !    h_interstitial = <P_phi|j_l> Y_L^*(R(k+K))*exp(i*(k+K)*r_latom)*exp(i*spin_phase)/Rmt^2
                                h_interstitial(M,i-ii+1,is) = phi_jl(i,iind)*h_yl(M,i-ii+1)/RMT(jatom)**2!*ph_spin(is)
                             ENDDO
                          enddo
                       endif
                       !!!  New
                       ibb=min(iblock,isize-ii+1)
                       ldb=nmat
                       lda=2*LMAX2+1
                       ldc=lda
                       isize_=2*l+1
                       !---- h_alyl(2*lmax+1,iblock,is)  contains rotated Apw's, such that chi(r) = (Apw*u(r) + Bpw*udot(r))*Ylm(r)
                       !---- h_blyl(2*lmax+1,iblock,is)  contains rotated Bpw's
                       !---- A(:N,iband,is)              contains eigenvectors, also named C(k+G,inad,is)
                       !---- alm[2*l+1,iband][irf=1,is] += sum_{iK\in block} h_alyl[2*l+1,iK][is]*A[iK,iband][is]
                       !---- alm[2*l+1,iband][irf=2,is] += sum_{iK\in block} h_blyl[2*l+1,iK][is]*A[iK,iband][is]
                       !----
                       !---- The results is:
                       !---- alm[lm,iband][1,is] = sum_G Apw(lm,is,K+G) * C(k+G,iband,is)
                       !---- alm[lm,iband][2,is] = sum_G Bpw(lm,is,K+G) * C(k+G,iband,is)
                       !---- Where C(k+G,iband,is) are eigenvectors, and Apw and Bpw are expansion coefficients defined in Shick et.al., PRB 60, 10763 (1999).
                       call zgemm('N','N',isize_,nemax-nemin+1,ibb,(1.d0,0.d0), h_alyl(1,1,is),lda,a(ii,nemin,is),ldb,(1.d0,0.d0), alm(1,nemin,1,is),ldc)
                       call zgemm('N','N',isize_,nemax-nemin+1,ibb,(1.d0,0.d0), h_blyl(1,1,is),lda,a(ii,nemin,is),ldb,(1.d0,0.d0), alm(1,nemin,2,is),ldc)
                       if (abs(projector).eq.5 .and. nonzero_interst) then !! The new
                          call zgemm('N','N',isize_,nbands,ibb,(1.d0,0.d0),h_interstitial(1,1,is),lda,a(ii,nemin,is),ldb,(1.d0,0.d0), a_interstitial(1,1,is),ldc)
                       endif !! The new
                    ENDDO !----------- over both spins ---------!
                 ENDDO    !----------- over iblock -------------!

                 !-------------- Adds localized orbitals to alm. -------------------!
                 if (nlo.ne.0) then
                    call lomain(nemin,nemax,lfirst,latom,n,jatom,isym,L,iso,crotloc_x_BR1)
                 end if

                 !do m=1,2*l+1
                 !   do num=nemin,nemax
                 !      do irf=1,nrf
                 !         WRITE(6,'(3I4,2x,4F10.6)') m,num,irf, alm(m,num,irf,1), alm(m,num,irf,2)
                 !      enddo
                 !   enddo
                 !enddo

                 !--------- FACCFAC = (i)^l*4.*PI*RMT(JATOM)**2/SQRT(VOL) --------!
                 CFAC = FAC*IMAG**L      !------  (i)^l*fac -----!

                 if (iso.eq.2) then
                    ALML(l,:(2*L+1),nemin:nemax,:nrf,1) = ALM(:(2*L+1),nemin:nemax,:nrf,1)*(CFAC*Rispin(1,1)) + ALM(:(2*L+1),nemin:nemax,:nrf,2)*(CFAC*Rispin(2,1))
                    ALML(l,:(2*L+1),nemin:nemax,:nrf,2) = ALM(:(2*L+1),nemin:nemax,:nrf,1)*(CFAC*Rispin(1,2)) + ALM(:(2*L+1),nemin:nemax,:nrf,2)*(CFAC*Rispin(2,2))
                    !ALML(l,:(2*L+1),nemin:nemax,:nrf,:iso) = ALM(:(2*L+1),nemin:nemax,:nrf,:iso)*CFAC
                    if (abs(projector).eq.5 .and. nonzero_interst) then
                       al_interstitial(:(2*L+1),:nbands,1,lcase) = a_interstitial(:(2*L+1),:nbands,1)*(CFAC*Rispin(1,1)) + a_interstitial(:(2*L+1),:nbands,2)*(CFAC*Rispin(2,1))
                       al_interstitial(:(2*L+1),:nbands,2,lcase) = a_interstitial(:(2*L+1),:nbands,1)*(CFAC*Rispin(1,2)) + a_interstitial(:(2*L+1),:nbands,2)*(CFAC*Rispin(2,2))
                       !al_interstitial(:(2*L+1),:nbands,:iso,lcase) = a_interstitial(:(2*L+1),:nbands,:iso)*CFAC
                    endif
                 else
                    ALML(l,:(2*L+1),nemin:nemax,:nrf,1) = ALM(:(2*L+1),nemin:nemax,:nrf,1)*CFAC
                    if (abs(projector).eq.5 .and. nonzero_interst) al_interstitial(:(2*L+1),:nbands,1,lcase) = a_interstitial(:(2*L+1),:nbands,1)*CFAC
                 endif
              enddo  !--------- end of atom L loop (alm)  ----------------!
              !----------------  ALML(l,m,iband,ifr,ispin) contains (A*C,B*C) for all l,m,iband,ifr=(1,2),ispin --------!
              !---------- Computing the DMFT transformation --------------!
              do lcase=1,nl(icase)

                 l1=ll(icase,lcase)
                 nind=(2*l1+1)*iso
                 idt = idet(isym)
                 icix = cix(icase,lcase)
                 iorb = iorbital(icase,lcase)

                 if (icix.gt.0) then
                    URx=0
                    N1=2*l1+1

                    if (abs(projector).eq.1) then
                       ! For correlated orbitals, we will also create transformation matrix U, such that P(LL'ij)=U^\dagger*U
                       !$OMP PARALLEL DO SHARED(URx) PRIVATE(num1,is1,m1,lms1,ind1,cc,irf1) SCHEDULE(STATIC)
                       do num1=1,nbands      !nemin,nemax     ! over bands
                          do is1=1,iso            ! over spin-1
                             do m1=-l1,l1           ! over m-1
                                lms1=l1+1+m1
                                ind1=l1+1+idt*m1+N1*(is1-1)
                                cc=0
                                do irf1=1,nrf
                                   cc = cc + alml(l1,lms1,num1+nemin-1,irf1,is1)*ri_mat(irf1,1,l1,is1)*(idt)**m1
                                enddo
                                if (idt.gt.0) then
                                   URx(num1,ind1) = cc
                                else
                                   URx(num1,ind1) = conjg(cc)
                                endif
                             enddo
                          enddo
                       enddo
                       !$OMP END PARALLEL DO
                    elseif (abs(projector).eq.2 .or. abs(projector).eq.4) then
                       ! L -> l1
                       ! nrfmax -> nrf
                       ! For correlated orbitals, we will also create transformation matrix U, such that P(LL'ij)=U^\dagger*U
                       !$OMP PARALLEL DO SHARED(URx) PRIVATE(num1,is1,m1,lms1,ind1,cc,sm1,sm2,irf1,irf2,ff) SCHEDULE(STATIC)
                       do num1=1,nbands      !nemin,nemax     ! over bands
                          do is1=1,iso            ! over spin-1
                             do m1=-l1,l1           ! over m-1
                                lms1=l1+1+m1
                                ind1=l1+1+idt*m1+N1*(is1-1)
                                cc=0
                                sm1=0
                                sm2=0
                                do irf1=1,nrf
                                   do irf2=1,nrf
                                      ff = alml(l1,lms1,num1+nemin-1,irf1,is1)*conjg(alml(l1,lms1,num1+nemin-1,irf2,is1))
                                      sm1 = sm1 + ff * ri_mat(irf1,irf2,l1,is1)
                                      sm2 = sm2 + ff * ri_mat(irf1,1,l1,is1) * ri_mat(1,irf2,l1,is1)
                                   enddo
                                   cc = cc + alml(l1,lms1,num1+nemin-1,irf1,is1)*ri_mat(irf1,1,l1,is1)*(idt)**m1
                                enddo
                                cc = cc * sqrt(abs(sm1/sm2))
                                if (idt.gt.0) then
                                   URx(num1,ind1) = cc
                                else
                                   URx(num1,ind1) = conjg(cc)
                                endif
                             enddo
                          enddo
                       enddo
                       !$OMP END PARALLEL DO
                    elseif (abs(projector).eq.3) then
                       ! For correlated orbitals, we will also create transformation matrix U, such that P(LL'ij)=U^\dagger*U
                       !$OMP PARALLEL DO SHARED(URx) PRIVATE(num1,is1,m1,lms1,ind1,cc,corc,irf1,irf2) SCHEDULE(STATIC)
                       do num1=1,nbands      !nemin,nemax     ! over bands
                          do is1=1,iso            ! over spin-1
                             do m1=-l1,l1           ! over m-1
                                lms1=l1+1+m1
                                ind1=l1+1+idt*m1+N1*(is1-1)
                                irf1 = 1
                                cc = alml(l1,lms1,num1+nemin-1,irf1,is1)*sqrt(ri_mat(irf1,irf1,l1,is1))*(idt)**m1
                                ! Adding correction due to \cdot{u} and local orbitals
                                ! The results is  A(i,L,irf=1)*sqrt(o(L)) * sqrt( 1 + \sum_{irf>1} |A(i,L,irf)|^2*o(L,irf))
                                if ( abs(cc)>1e-6 ) then
                                   corc = 0
                                   do irf1=1,nrf
                                      do irf2=1,nrf
                                         if (irf1.eq.1 .and. irf2.eq.1) cycle
                                         corc = corc + alml(l1,lms1,num1+nemin-1,irf1,is1)*conjg(alml(l1,lms1,num1+nemin-1,irf2,is1))*ri_mat(irf1,irf2,l1,is1)
                                      enddo
                                   enddo
                                   cc = cc*sqrt(1 + abs(corc)/abs(cc)**2)
                                endif
                                if (idt.gt.0) then
                                   URx(num1,ind1) = cc
                                else
                                   URx(num1,ind1) = conjg(cc)
                                endif
                             enddo
                          enddo
                       enddo
                       !$OMP END PARALLEL DO
                    elseif (abs(projector).eq.5) then  ! fixed projector
                       ii = al_ucase(icase,lcase)
                       !$OMP PARALLEL DO SHARED(URx,rix_mat) PRIVATE(num1,is1,m1,lms1,ind1,cc,irf1) SCHEDULE(STATIC)
                       do num1=1,nbands
                          do is1=1,iso
                             do m1=-l1,l1
                                lms1=l1+1+m1
                                ind1=l1+1+idt*m1+N1*(is1-1)
                                cc=0
                                do irf1=1,nrf
                                   cc = cc + alml(l1,lms1,num1+nemin-1,irf1,is1)*rix_mat(irf1,ii)*(idt)**m1
                                enddo
                                !!! New
                                cc = cc + al_interstitial(lms1,num1,is1,lcase)*(idt)**m1
                                !!! New
                                if (idt.gt.0) then
                                   URx(num1,ind1) = cc
                                else
                                   URx(num1,ind1) = conjg(cc)
                                endif
                             enddo
                          enddo
                       enddo
                       !$OMP END PARALLEL DO
                    else
                       print *, 'Only projector=[1,2,3,4,5,-1,-2,-3,-4,-5] is allowed!'
                       stop
                    endif

                    nind1 = nind
                    iorb1=iorb

                    !if (iso.eq.2) then
                    !   no1 = nind1/iso
                    !   ALLOCATE( Uu(nbands,no1), Ud(nbands,no1))
                    !   tuu = Rspin(1,1,iorb1)
                    !   tud = Rspin(1,2,iorb1)
                    !   tdu = Rspin(2,1,iorb1)
                    !   tdd = Rspin(2,2,iorb1)
                    !   Uu(:,:) = URx(:,1:no1)
                    !   Ud(:,:) = URx(:,no1+1:2*no1)
                    !   URx(:,1:no1)       = Uu*tuu + Ud*tdu
                    !   URx(:,no1+1:2*no1) = Uu*tud + Ud*tdd
                    !   DEALLOCATE( Uu, Ud )
                    !endif

                    DO iorb2=1,norbitals
                       if ( cix_orb(iorb2).NE. icix ) CYCLE
                       nind2 = nindo(iorb2)
                       call zgemm('N','C', nbands, nind2, nind1, (1.d0,0.d0), URx,nbands, cfX(:,:,iorb2,iorb1),maxdim2, (0.d0,0.d0), tmp,nbands)
                       DMFTU(:,:nind2,iorb2) = DMFTU(:,:nind2,iorb2) + conjg(tmp(:,:nind2))
                    ENDDO

                 endif
              enddo  ! over L-case
           ENDDO     ! over atom-cases
           !-------- computing correlated green's function --------!

           Olapmk=0
           DO iorb1=1,norbitals
              icix = cix_orb(iorb1)
              if ( icix.EQ.0 ) CYCLE
              nind1 = nindo(iorb1)
              DO iorb2=1,norbitals
                 if ( icix.NE.cix_orb(iorb2) ) CYCLE
                 nind2 = nindo(iorb2)

                 call zgemm('C','N', nind1, nind2, nbands, (1.d0,0.d0), DMFTU(:,:,iorb1),nbands, DMFTU(:,:,iorb2),nbands, (0.d0,0.d0), olp,maxdim2)

                 do ind1=1,nind1
                    do ind2=1,nind2
                       Olapmk( iSx(ind1,iorb1), iSx(ind2,iorb2), icix ) = olp(ind1,ind2)
                    enddo
                 enddo
              enddo
           ENDDO
           Olapm(:,:,:) = Olapm(:,:,:) + Olapmk(:,:,:) * (mweight(ikp)/tweight/nsymop2) ! proper weight for the reducible k-point
        ENDDO !------  isym: over star of the irreducible k-point ----!

        deallocate( DMFTU)
        deallocate( URx, tmp)

        if (abs(projector).eq.5) then
           deallocate( phi_jl )
           deallocate( h_interstitial )
           deallocate( a_interstitial, al_interstitial )
        endif

        WRITE(*,'(I3,A,1x,I5,1x,A,I4,A,I2)') myrank, ') Finished k-point number', ikp, 'with #bands=', nbands, ' Nstar=', nsymop2
        call flush(6)
     ENDDO !---- over reducible k-point: ikp

998  CONTINUE !---- irreducible k-points end (jump) from reading
     !--- end k-points

     DO is=1,iso    !------ over up/dn ---------------------!
        itape=8+is
        if (vector_para) then
           close(itape)
        else
           rewind(itape)
        endif
     ENDDO

  ENDDO ! over different vector files

  CALL AllReduce_MPI(Olapm, maxdim, ncix)
  if (myrank.eq.master) then
     WRITE(*,*) 'Renormalizing Gloc to account for the interstitials'
  endif

  if (Qsymmetrize) Call SymmetrizeOverlap(Olapm, cfX, cix_orb, cixdim, iSx, iorbital, nindo, norbitals, maxdim2)

  if (SIMPLE) then
     Olapc=0
     DO icix=1,ncix
        cixdm = cixdim(icix)
        !---- Olap to vector form, and s_oo to matrix form
        DO ip=1,cixdm
           do iq=1,cixdm
              it = abs(Sigind(ip,iq,icix))
              if (it.gt.0) then
                 Olapc(it, icix) =  Olapc(it, icix) + real(Olapm(ip,iq,icix))
              endif
           enddo
        ENDDO
        DO it=1,csize(icix)
           Olapc(it, icix) = Olapc(it, icix)/noccur(it,icix)
        ENDDO
     ENDDO
     if (myrank.eq.master) then
        WRITE(*,'(A)',advance='no') '  Z due to finite E-window='
        do icix=1,ncix
           do it=1,csize(icix)
              WRITE(*,'(F15.12)',advance='no') Olapc(it,icix)
           enddo
        enddo
        WRITE(*,*)
     endif
  else if (.True.) then
     !!! We copute SOlapm = 1/sqrt(Olapm) for later use
     lwork  = 4*maxdim + maxdim*maxdim
     lrwork = 10*maxdim + 2*maxdim*maxdim + 1
     liwork = 3+10*maxdim
     allocate( work(lwork), rwork(lrwork), iwork(liwork) )

     SOlapm(:,:,:) = 0.d0
     DO icix=1,ncix
        cixdm = cixdim(icix)

        allocate( cind(cixdm) )
        cixdms=0 ! real dimension, excluding states which must be projected out, because they are treated as non-correlated
        cind=0   ! In most cases, cind(i)=i, however, when some states are projected out, cind points to smaller block
        DO ip=1,cixdm
           if (Sigind(ip,ip,icix) .ne. 0) then
              cixdms = cixdms + 1
              cind(ip) = cixdms
           endif
        ENDDO
        allocate( cini(cixdms))
        do ip=1,cixdm
           if (cind(ip).gt.0) cini(cind(ip))=ip
        enddo
        deallocate( cind )

        allocate( olocf(cixdms,cixdms), ws(cixdms), tmp1(cixdms,cixdms), solap(cixdms,cixdms) )
        do ip=1,cixdms
           do iq=1,cixdms
              olocf(ip,iq) = Olapm(cini(ip),cini(iq),icix)
           enddo
        enddo
        CALL ZHEEVD('V','U', cixdms, olocf, cixdms, ws, work, lwork, rwork, lrwork, iwork, liwork, info )
        if (info .ne. 0) then
           print *, 'Diagonalization of renormalization factor failed. Info-zheevd=', info
        endif
        do ip=1,cixdms
           tmp1(:,ip) = olocf(:,ip)*(1./sqrt(abs(ws(ip))))
        enddo
        call zgemm('N','C', cixdms, cixdms, cixdms, (1.d0,0.d0), tmp1, cixdms, olocf, cixdms, (0.d0,0.d0), solap, cixdms)

        do ip=1,cixdm
           SOlapm(ip,ip,icix)=1.0d0
        enddo
        do ip=1,cixdms
           do iq=1,cixdms
              SOlapm(cini(ip),cini(iq),icix) = solap(ip,iq)
           enddo
        enddo

        if (myrank.eq.master) then
           !WRITE(*,'(A)',advance='no') '  Z due to finite E-window='
           !do it=1,cixdm
           !   WRITE(*,'(F15.12)',advance='no') ws(it)
           !enddo
           !WRITE(*,*)
           WRITE(*,'(A)',advance='no') 'Z to renormalize='
           do ip=1,cixdms
              WRITE(*,'(F16.10)',advance='no') 1./dble(SOlapm(cini(ip),cini(ip),icix))
           enddo
           WRITE(*,*)
        endif

        deallocate( cini )
        deallocate( olocf, ws, tmp1, solap )
     ENDDO

     deallocate( work, rwork, iwork )

     ! Save for later for dmft2 step
     open(996, FILE='SOlapm.dat', status='unknown')
     DO icix=1,ncix
        cixdm = cixdim(icix)
        WRITE(996,*) icix
        do ind1=1,cixdm
           do ind2=1,cixdm
              WRITE(996,'(2F20.16)',advance='no')  SOlapm(ind1,ind2,icix)
           enddo
           WRITE(996,*)
        enddo
     ENDDO
     close(996)

  else
     !!! We copute SOlapm = 1/sqrt(Olapm) for later use
     lwork  = 4*maxdim + maxdim*maxdim
     lrwork = 10*maxdim + 2*maxdim*maxdim + 1
     liwork = 3+10*maxdim
     allocate( work(lwork), rwork(lrwork), iwork(liwork) )
     allocate( olocf(maxdim,maxdim), ws(maxdim), tmp1(maxdim,maxdim) )

     DO iorb1=1,norbitals
        icix = cix_orb(iorb1)
        if ( icix.EQ.0 ) CYCLE
        nind1 = nindo(iorb1)
        olocf=0
        do ind1=1,nind1
           do ind2=1,nind1
              olocf(ind1,ind2) = Olapm( iSx(ind1,iorb1), iSx(ind2,iorb1), icix )
           enddo
        enddo

        CALL ZHEEVD('V','U', nind1, olocf, maxdim, ws, work, lwork, rwork, lrwork, iwork, liwork, info )
        if (info .ne. 0) then
           print *, 'Diagonalization of renormalization factor failed. Info-zheevd=', info
        endif

        do ip=1,nind1
           tmp1(:,ip) = olocf(:,ip)*(1./sqrt(abs(ws(ip))))
        enddo
        call zgemm('N','C', nind1, nind1, nind1, (1.d0,0.d0), tmp1, maxdim, olocf,maxdim, (0.d0,0.d0), SOlapm(:,:,icix),maxdim)

        if (myrank.eq.master) then
           WRITE(*,*) 'Z due to finite E-window=', ws(:nind1)
        endif

     ENDDO
     deallocate( work, rwork, iwork )
     deallocate( olocf, ws, tmp1 )
  endif

  DEALLOCATE( a_real ) !---  for eigenvectors -------------!
  deallocate( olp )
  deallocate( Olapmk )

  !CALL cputim(t2c)
  !CALL walltim(t2w)
  !time1c = time1c + t2c-t1c
  !time1w = time1w + t2w-t1w

  !print *, 'total time overlap:', time1c, time1w
  !print *, 'overlap interst. proj:', time2c, time2w

END SUBROUTINE cmp_overlap


text:
#compdef $CAL_HOME/shell/system/calfind.zsh

local -a commands

while {read alias} {
    commands+=(${alias%%:*})
} < $CAL_HPATH

compadd $commands

ext: f

max_stars_repo_name: calject/cal-shell

max_issues_repo_name: calject/cal-shell

max_forks_repo_name: calject/cal-shell


text:
      SUBROUTINE reset_params
      USE precon2d, ONLY: ictrl_prec2d
      USE vmec_main, ONLY: iequi, ivac, ftolv, fsqr, fsqz, fsq, dp,
     1                     res0, delt0r, iter1, iter2, ijacob, irst
      USE vmec_input, ONLY: delt
      USE timer_sub, ONLY: timer
      IMPLICIT NONE
C-----------------------------------------------

!     2d preconditioner
      ictrl_prec2d = 0

      iequi = 0
      ivac  = -1

      fsqr = 1
      fsqz = 1
      ftolv = fsqr

      fsq   = 1
      iter2 = 1
      iter1 = iter2
      ijacob = 0
      irst = 1
      res0 = -1       !move to vmec_main, remove from runvmec SAVE
      delt0r = delt   !move delt0 to vmec_main from runvmec SAVE (renamed delt0r)

      timer =  0

      END SUBROUTINE reset_params

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT


text:
PROGRAM stencilmark
      IMPLICIT NONE
      real(kind=8) :: A0(258,258,258)
      real(kind=8) :: A1(258,258,258)
      real(kind=8) :: A2(258,258,258)
      real(kind=8) :: A3(258,258,258)
      real(kind=8) :: T0(258,258,258)
      real(kind=8) :: T1(258,258,258)
      real :: tstart, tfinish, tmin
      integer i,j,k,trial

      ! initialize arrays
      do i=1,258
        do j=1,258
          do k=1,258
            A0(i,j,k) = RAND()
            A1(i,j,k) = RAND()
            A2(i,j,k) = RAND()
            A3(i,j,k) = RAND()
          enddo
        enddo
      enddo

      tmin = 99999999.0_8
      do trial=1,5

        ! clear the cache
        call clear_cache(T0,T1)

        ! call the stencil
        call cpu_time(tstart)
        call stencil(A0, A1, A2, A3)
        call cpu_time(tfinish)
        print '("Time = ",f6.3," seconds.")',tfinish-tstart

        if (tfinish-tstart < tmin) then
          tmin = tfinish-tstart
        endif
      enddo
      print '("Min Time = ",f6.3," seconds.")',tmin

CONTAINS
      SUBROUTINE clear_cache(T0,T1)
        real(kind=8) :: T0(258,258,258)
        real(kind=8) :: T1(258,258,258)
        integer :: i,j,k

        do i=1,258
          do j=1,258
            do k=1,258
              T0(i,j,k) = RAND()
              T1(i,j,k) = RAND()
            enddo
          enddo
        enddo
      END SUBROUTINE

END PROGRAM stencilmark

max_stars_repo_name: uwplse/stng

max_issues_repo_name: uwplse/stng

max_forks_repo_name: uwplse/stng


!---------------------------------------------------!
! Copyright (c) 2016 Shunsuke A. Sato               !
! Released under the MIT license                    !
! https://opensource.org/licenses/mit-license.php   !
!---------------------------------------------------!
!-------10--------20--------30--------40--------50--------60--------70--------80--------90
subroutine energy_calc(Eelec,Eelemag,Etot)
  use global_variables_ms
  implicit none
  real(8),intent(out) :: Eelec,Eelemag,Etot
  real(8) :: Eelec_l,Et,Bt
  integer :: ix,im,ikz,ikr

  Eelec_l = 0d0
  do im = Mx_s,Mx_e
    do ikz = -NKz,NKz
      do ikr = 1,NKr
        Eelec_l = Eelec_l + deps(ikr,ikz,im)*abs(zCt(2,ikr,ikz,im))**2*kr(ikr)
      end do
    end do
  end do
  Eelec_l = Eelec_l * 2d0/((2d0*pi)**3)*(2d0*pi*dkr*dkz)
  call MPI_ALLREDUCE(Eelec_l,Eelec,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,ierr)

  Eelemag = 0d0
  do ix = Nx_L+1,Nx_R-1
    Et = -0.5d0*(Az_new(ix)-Az_old(ix))/dt
    Bt = -0.5d0*(Az(ix+1)-Az(ix-1))*c_light/Hx
    Eelemag = Eelemag + (Et**2+Bt**2)
  end do

  Etot = Eelec + Eelemag
!!== Unit should be consistent between Eelec and Eelemag ==!!

end subroutine energy_calc


MODULE module_cbmz_addemiss
   integer, parameter :: cbmz_addemiss_masscheck = -1
CONTAINS
   subroutine cbmz_addemiss_anthro( id, dtstep, dz8w, config_flags, &
               rho_phy, chem,emis_ant,alt, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte )
  USE module_configure
  USE module_state_description
  USE module_data_radm2
  IMPLICIT NONE
   TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
   INTEGER, INTENT(IN ) :: id, &
                                  ids,ide, jds,jde, kds,kde, &
                                  ims,ime, jms,jme, kms,kme, &
                                  its,ite, jts,jte, kts,kte
   REAL, INTENT(IN ) :: dtstep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
         INTENT(INOUT ) :: chem
   REAL, DIMENSION( ims:ime, kms:config_flags%kemit, jms:jme,num_emis_ant),&
         INTENT(IN ) :: &
                         emis_ant
   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
          INTENT(IN ) :: dz8w, rho_phy,alt
   integer :: i,j,k
   integer :: iflg
   real, parameter :: efact1 = 1.0/60.0
   real :: conv,conv2,conv3
   double precision :: chem_sum(num_chem)
      if (cbmz_addemiss_masscheck > 0) call addemiss_masscheck( &
               id, config_flags, 1, 'cbmz_addemiss', &
               dtstep, efact1, dz8w, chem, chem_sum, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte, &
               21, &
               emis_ant(ims,kms,jms,p_e_so2),emis_ant(ims,kms,jms,p_e_no), &
               emis_ant(ims,kms,jms,p_e_co),emis_ant(ims,kms,jms,p_e_eth), &
               emis_ant(ims,kms,jms,p_e_hc3),emis_ant(ims,kms,jms,p_e_hc5), &
               emis_ant(ims,kms,jms,p_e_hc8),emis_ant(ims,kms,jms,p_e_xyl), &
               emis_ant(ims,kms,jms,p_e_ol2),emis_ant(ims,kms,jms,p_e_olt), &
               emis_ant(ims,kms,jms,p_e_oli),emis_ant(ims,kms,jms,p_e_tol), &
               emis_ant(ims,kms,jms,p_e_csl),emis_ant(ims,kms,jms,p_e_hcho), &
               emis_ant(ims,kms,jms,p_e_ald),emis_ant(ims,kms,jms,p_e_ket), &
               emis_ant(ims,kms,jms,p_e_ora2),emis_ant(ims,kms,jms,p_e_nh3), &
               emis_ant(ims,kms,jms,p_e_no2),emis_ant(ims,kms,jms,p_e_ch3oh), &
               emis_ant(ims,kms,jms,p_e_c2h5oh))
      do 100 j=jts,jte
      do 100 i=its,ite
      DO k=kts,min(config_flags%kemit,kte)
        conv = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
        conv2 = (dtstep/dz8w(i,k,j))/rho_phy(i,k,j)
        conv3 = (dtstep/dz8w(i,k,j))*alt(i,k,j)*28/250*1e-3
        chem(i,k,j,p_csl) = chem(i,k,j,p_csl) &
                         +emis_ant(i,k,j,p_e_csl)*conv
        chem(i,k,j,p_so2) = chem(i,k,j,p_so2) &
                         +emis_ant(i,k,j,p_e_so2)*conv
        chem(i,k,j,p_no) = chem(i,k,j,p_no) &
                         +emis_ant(i,k,j,p_e_no)*conv
        chem(i,k,j,p_ald) = chem(i,k,j,p_ald) &
                         +emis_ant(i,k,j,p_e_ald)*conv
        chem(i,k,j,p_hcho) = chem(i,k,j,p_hcho) &
                         +emis_ant(i,k,j,p_e_hcho)*conv
        chem(i,k,j,p_ora2) = chem(i,k,j,p_ora2) &
                         +emis_ant(i,k,j,p_e_ora2)*conv
        chem(i,k,j,p_nh3) = chem(i,k,j,p_nh3) &
                         +emis_ant(i,k,j,p_e_nh3)*conv
        chem(i,k,j,p_eth) = chem(i,k,j,p_eth) &
                         +emis_ant(i,k,j,p_e_eth)*conv
        chem(i,k,j,p_co) = chem(i,k,j,p_co) &
                         +emis_ant(i,k,j,p_e_co)*conv
        chem(i,k,j,p_ol2) = chem(i,k,j,p_ol2) &
                         +emis_ant(i,k,j,p_e_ol2)*conv
        chem(i,k,j,p_olt) = chem(i,k,j,p_olt) &
                         +emis_ant(i,k,j,p_e_olt)*conv
        chem(i,k,j,p_oli) = chem(i,k,j,p_oli) &
                         +emis_ant(i,k,j,p_e_oli)*conv
        chem(i,k,j,p_tol) = chem(i,k,j,p_tol) &
                         +emis_ant(i,k,j,p_e_tol)*conv
        chem(i,k,j,p_xyl) = chem(i,k,j,p_xyl) &
                         +emis_ant(i,k,j,p_e_xyl)*conv
        chem(i,k,j,p_ket) = chem(i,k,j,p_ket) &
                         +emis_ant(i,k,j,p_e_ket)*conv
        chem(i,k,j,p_no2) = chem(i,k,j,p_no2) &
                         +emis_ant(i,k,j,p_e_no2)*conv
        chem(i,k,j,p_ch3oh) = chem(i,k,j,p_ch3oh) &
                         +emis_ant(i,k,j,p_e_ch3oh)*conv
        chem(i,k,j,p_c2h5oh) = chem(i,k,j,p_c2h5oh) &
                         +emis_ant(i,k,j,p_e_c2h5oh)*conv
        chem(i,k,j,p_iso) = chem(i,k,j,p_iso) &
                         +emis_ant(i,k,j,p_e_iso)*conv
         chem_select_2 : SELECT CASE( config_flags%chem_opt )
         CASE(CBMZ_MOSAIC_8bIN_VBS9_KPP,CBMZ_MOSAIC_4BIN_VBS9_KPP)
        chem(i,k,j,p_pcg1_b_c) = chem(i,k,j,p_pcg1_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*1.06
        chem(i,k,j,p_pcg2_b_c) = chem(i,k,j,p_pcg2_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.12
        chem(i,k,j,p_pcg3_b_c) = chem(i,k,j,p_pcg3_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.18
        chem(i,k,j,p_pcg4_b_c) = chem(i,k,j,p_pcg4_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.28
        chem(i,k,j,p_pcg5_b_c) = chem(i,k,j,p_pcg5_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.36
        chem(i,k,j,p_pcg6_b_c) = chem(i,k,j,p_pcg6_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.6
        chem(i,k,j,p_pcg7_b_c) = chem(i,k,j,p_pcg7_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*1.0
        chem(i,k,j,p_pcg8_b_c) = chem(i,k,j,p_pcg8_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*1.5
        chem(i,k,j,p_pcg9_b_c) = chem(i,k,j,p_pcg9_b_c) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*2.4
        chem(i,k,j,p_pcg1_f_c) = chem(i,k,j,p_pcg1_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*1.06
        chem(i,k,j,p_pcg2_f_c) = chem(i,k,j,p_pcg2_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.12
        chem(i,k,j,p_pcg3_f_c) = chem(i,k,j,p_pcg3_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.18
        chem(i,k,j,p_pcg4_f_c) = chem(i,k,j,p_pcg4_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.28
        chem(i,k,j,p_pcg5_f_c) = chem(i,k,j,p_pcg5_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.36
        chem(i,k,j,p_pcg6_f_c) = chem(i,k,j,p_pcg6_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.6
        chem(i,k,j,p_pcg7_f_c) = chem(i,k,j,p_pcg7_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*1.0
        chem(i,k,j,p_pcg8_f_c) = chem(i,k,j,p_pcg8_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*1.5
        chem(i,k,j,p_pcg9_f_c) = chem(i,k,j,p_pcg9_f_c) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*2.4
        chem(i,k,j,p_pcg1_b_o) = chem(i,k,j,p_pcg1_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.318
        chem(i,k,j,p_pcg2_b_o) = chem(i,k,j,p_pcg2_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.036
        chem(i,k,j,p_pcg3_b_o) = chem(i,k,j,p_pcg3_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.054
        chem(i,k,j,p_pcg4_b_o) = chem(i,k,j,p_pcg4_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.084
        chem(i,k,j,p_pcg5_b_o) = chem(i,k,j,p_pcg5_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.108
        chem(i,k,j,p_pcg6_b_o) = chem(i,k,j,p_pcg6_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.18
        chem(i,k,j,p_pcg7_b_o) = chem(i,k,j,p_pcg7_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.30
        chem(i,k,j,p_pcg8_b_o) = chem(i,k,j,p_pcg8_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.45
        chem(i,k,j,p_pcg9_b_o) = chem(i,k,j,p_pcg9_b_o) &
        +(emis_ant(i,k,j,p_e_orgi_bb)/1.4+emis_ant(i,k,j,p_e_orgj_bb)/1.4)*conv3*0.72
        chem(i,k,j,p_pcg1_f_o) = chem(i,k,j,p_pcg1_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.053
        chem(i,k,j,p_pcg2_f_o) = chem(i,k,j,p_pcg2_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.006
        chem(i,k,j,p_pcg3_f_o) = chem(i,k,j,p_pcg3_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.009
        chem(i,k,j,p_pcg4_f_o) = chem(i,k,j,p_pcg4_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.014
        chem(i,k,j,p_pcg5_f_o) = chem(i,k,j,p_pcg5_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.018
        chem(i,k,j,p_pcg6_f_o) = chem(i,k,j,p_pcg6_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.03
        chem(i,k,j,p_pcg7_f_o) = chem(i,k,j,p_pcg7_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.05
        chem(i,k,j,p_pcg8_f_o) = chem(i,k,j,p_pcg8_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.075
        chem(i,k,j,p_pcg9_f_o) = chem(i,k,j,p_pcg9_f_o) &
        +(emis_ant(i,k,j,p_e_orgi_a)/1.4+emis_ant(i,k,j,p_e_orgj_a)/1.4)*conv3*0.12
         END SELECT chem_select_2
        iflg=0
        if(id.eq.1) then
          if(i.ge.82.and.i.le.85.and.j.ge.47.and.j.le.50.and.k.le.1) then
            iflg=1
          elseif(k.eq.1) then
            iflg=2
          endif
        endif
        if(id.eq.2) then
          if(i.ge.43.and.i.le.59.and.j.ge.30.and.j.le.56.and.k.le.1) then
            iflg=1
          elseif(k.eq.1) then
            iflg=2
          endif
        endif
        if(iflg.eq.0) chem(i,k,j,p_cobb) = chem(i,k,j,p_cobb) &
                         +emis_ant(i,k,j,p_e_co)*conv
        if(iflg.eq.1) chem(i,k,j,p_comc) = chem(i,k,j,p_comc) &
                         +emis_ant(i,k,j,p_e_co)*conv
        if(iflg.eq.2) chem(i,k,j,p_comx) = chem(i,k,j,p_comx) &
                         +emis_ant(i,k,j,p_e_co)*conv
        if(iflg.eq.0) chem(i,k,j,p_poabb) = chem(i,k,j,p_poabb) &
                         +emis_ant(i,k,j,p_e_orgi)*conv2 &
                         +emis_ant(i,k,j,p_e_orgj)*conv2
        if(iflg.eq.1) chem(i,k,j,p_poamc) = chem(i,k,j,p_poamc) &
                         +emis_ant(i,k,j,p_e_orgi)*conv2 &
                         +emis_ant(i,k,j,p_e_orgj)*conv2
        if(iflg.eq.2) chem(i,k,j,p_poamx) = chem(i,k,j,p_poamx) &
                         +emis_ant(i,k,j,p_e_orgi)*conv2 &
                         +emis_ant(i,k,j,p_e_orgj)*conv2
        if ( (config_flags%emiss_inpt_opt == EMISS_INPT_DEFAULT) .or. &
             (config_flags%emiss_inpt_opt == EMISS_INPT_PNNL_RS) ) then
            chem(i,k,j,p_par) = chem(i,k,j,p_par) &
                + conv* &
                  ( 0.4*emis_ant(i,k,j,p_e_ald) + 2.9*emis_ant(i,k,j,p_e_hc3) &
                  + 4.8*emis_ant(i,k,j,p_e_hc5) + 7.9*emis_ant(i,k,j,p_e_hc8) &
                  + 0.9*emis_ant(i,k,j,p_e_ket) + 2.8*emis_ant(i,k,j,p_e_oli) &
                  + 1.8*emis_ant(i,k,j,p_e_olt) + 1.0*emis_ant(i,k,j,p_e_ora2) )
        else
            chem(i,k,j,p_par) = chem(i,k,j,p_par) &
                + conv*emis_ant(i,k,j,p_e_hc5)
            chem(i,k,j,p_no2) = chem(i,k,j,p_no2) &
                + conv*emis_ant(i,k,j,p_e_no2)
            chem(i,k,j,p_ch3oh) = chem(i,k,j,p_ch3oh) &
                + conv*emis_ant(i,k,j,p_e_ch3oh)
            chem(i,k,j,p_c2h5oh) = chem(i,k,j,p_c2h5oh) &
                + conv*emis_ant(i,k,j,p_e_c2h5oh)
        end if
        if ( (config_flags%emiss_inpt_opt == EMISS_INPT_PNNL_MAM)) then
           chem(i,k,j,p_dms) = chem(i,k,j,p_dms) &
                + conv*emis_ant(i,k,j,p_e_dms)
        end if
      END DO
 100 continue
      if (cbmz_addemiss_masscheck > 0) call addemiss_masscheck( &
               id, config_flags, 2, 'cbmz_addemiss', &
               dtstep, efact1, dz8w, chem, chem_sum, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte, &
               21, &
               emis_ant(ims,kms,jms,p_e_so2),emis_ant(ims,kms,jms,p_e_no), &
               emis_ant(ims,kms,jms,p_e_co),emis_ant(ims,kms,jms,p_e_eth), &
               emis_ant(ims,kms,jms,p_e_hc3),emis_ant(ims,kms,jms,p_e_hc5), &
               emis_ant(ims,kms,jms,p_e_hc8),emis_ant(ims,kms,jms,p_e_xyl), &
               emis_ant(ims,kms,jms,p_e_ol2),emis_ant(ims,kms,jms,p_e_olt), &
               emis_ant(ims,kms,jms,p_e_oli),emis_ant(ims,kms,jms,p_e_tol), &
               emis_ant(ims,kms,jms,p_e_csl),emis_ant(ims,kms,jms,p_e_hcho), &
               emis_ant(ims,kms,jms,p_e_ald),emis_ant(ims,kms,jms,p_e_ket), &
               emis_ant(ims,kms,jms,p_e_ora2),emis_ant(ims,kms,jms,p_e_nh3), &
               emis_ant(ims,kms,jms,p_e_no2),emis_ant(ims,kms,jms,p_e_ch3oh), &
               emis_ant(ims,kms,jms,p_e_c2h5oh))
   END subroutine cbmz_addemiss_anthro
  subroutine cbmz_addemiss_bio( id, dtstep, dz8w, config_flags, &
        rho_phy, chem, e_bio, ne_area, e_iso, &
        ids,ide, jds,jde, kds,kde, &
        ims,ime, jms,jme, kms,kme, &
        its,ite, jts,jte, kts,kte )
  USE module_configure
  USE module_state_description
  USE module_data_radm2
  USE module_aerosols_sorgam
  IMPLICIT NONE
   TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
   INTEGER, INTENT(IN ) :: id, ne_area, &
                                  ids,ide, jds,jde, kds,kde, &
                                  ims,ime, jms,jme, kms,kme, &
                                  its,ite, jts,jte, kts,kte
   REAL, INTENT(IN ) :: dtstep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
         INTENT(INOUT ) :: chem
   REAL, DIMENSION( ims:ime, jms:jme,ne_area ), &
         INTENT(IN ) :: e_bio
   REAL, DIMENSION( ims:ime, kms:config_flags%kemit, jms:jme ), &
         INTENT(IN ) :: e_iso
   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
          INTENT(IN ) :: dz8w, rho_phy
   integer i,j,k,n
   real, parameter :: efact1 = 1.0/60.0
   double precision :: chem_sum(num_chem)
   if (config_flags%bio_emiss_opt == GUNTHER1) then
      if (cbmz_addemiss_masscheck > 0) call addemiss_masscheck( &
               id, config_flags, 1, 'cbmz_addemiss_bioaa', &
               dtstep, efact1, dz8w, chem, chem_sum, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte, &
               13, &
               e_bio(ims,jms,lald), e_bio(ims,jms,lhc3), &
               e_bio(ims,jms,lhc5), e_bio(ims,jms,lhc8), &
               e_bio(ims,jms,lhcho), e_bio(ims,jms,liso), &
               e_bio(ims,jms,lket), e_bio(ims,jms,lno), &
               e_bio(ims,jms,loli), e_bio(ims,jms,lolt), &
               e_bio(ims,jms,lora1), e_bio(ims,jms,lora2), &
               e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl) )
      do j=jts,jte
      do i=its,ite
        chem(i,kts,j,p_so2) = chem(i,kts,j,p_so2) &
                          + e_bio(i,j,lso2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_sulf) = chem(i,kts,j,p_sulf) &
                          + e_bio(i,j,lsulf)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_no2) = chem(i,kts,j,p_no2) &
                          + e_bio(i,j,lno2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_no) = chem(i,kts,j,p_no) &
                          + e_bio(i,j,lno)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_o3) = chem(i,kts,j,p_o3) &
                          + e_bio(i,j,lo3)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_hno3) = chem(i,kts,j,p_hno3) &
                          + e_bio(i,j,lhno3)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_h2o2) = chem(i,kts,j,p_h2o2) &
                          + e_bio(i,j,lh2o2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_ald) = chem(i,kts,j,p_ald) &
                          + e_bio(i,j,lald)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_hcho) = chem(i,kts,j,p_hcho) &
                          + e_bio(i,j,lhcho)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_op1) = chem(i,kts,j,p_op1) &
                          + e_bio(i,j,lop1)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_op2) = chem(i,kts,j,p_op2) &
                          + e_bio(i,j,lop2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_ora1) = chem(i,kts,j,p_ora1) &
                          + e_bio(i,j,lora1)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_ora2) = chem(i,kts,j,p_ora2) &
                          + e_bio(i,j,lora2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_nh3) = chem(i,kts,j,p_nh3) &
                          + e_bio(i,j,lnh3)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_n2o5) = chem(i,kts,j,p_n2o5) &
                          + e_bio(i,j,ln2o5)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_no2) = chem(i,kts,j,p_no2) &
                          + e_bio(i,j,lno2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_pan) = chem(i,kts,j,p_pan) &
                          + e_bio(i,j,lpan)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_eth) = chem(i,kts,j,p_eth) &
                          + e_bio(i,j,leth)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_co) = chem(i,kts,j,p_co) &
                          + e_bio(i,j,lco)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_ol2) = chem(i,kts,j,p_ol2) &
                          + e_bio(i,j,lol2)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_olt) = chem(i,kts,j,p_olt) &
                          + e_bio(i,j,lolt)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_oli) = chem(i,kts,j,p_oli) &
                          + e_bio(i,j,loli)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_tol) = chem(i,kts,j,p_tol) &
                          + e_bio(i,j,ltol)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_xyl) = chem(i,kts,j,p_xyl) &
                          + e_bio(i,j,lxyl)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_hono) = chem(i,kts,j,p_hono) &
                          + e_bio(i,j,lhono)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_hno4) = chem(i,kts,j,p_hno4) &
                          + e_bio(i,j,lhno4)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_ket) = chem(i,kts,j,p_ket) &
                          + e_bio(i,j,lket)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_mgly) = chem(i,kts,j,p_mgly) &
                          + e_bio(i,j,lmgly)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_onit) = chem(i,kts,j,p_onit) &
                          + e_bio(i,j,lonit)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_csl) = chem(i,kts,j,p_csl) &
                          + e_bio(i,j,lcsl)/(dz8w(i,kts,j)*60.)*dtstep
        chem(i,kts,j,p_iso) = chem(i,kts,j,p_iso) &
                          + e_bio(i,j,liso)/(dz8w(i,kts,j)*60.)*dtstep
      end do
      end do
      do j = jts, jte
      do i = its, ite
         chem(i,kts,j,p_par) = chem(i,kts,j,p_par) &
             + (dtstep/(dz8w(i,kts,j)*60.))* &
               ( 0.4*e_bio(i,j,lald) + 2.9*e_bio(i,j,lhc3) &
               + 4.8*e_bio(i,j,lhc5) + 7.9*e_bio(i,j,lhc8) &
               + 0.9*e_bio(i,j,lket) + 2.8*e_bio(i,j,loli) &
               + 1.8*e_bio(i,j,lolt) + 1.0*e_bio(i,j,lora2) )
      end do
      end do
      if (cbmz_addemiss_masscheck > 0) call addemiss_masscheck( &
               id, config_flags, 2, 'cbmz_addemiss_bioaa', &
               dtstep, efact1, dz8w, chem, chem_sum, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte, &
               13, &
               e_bio(ims,jms,lald), e_bio(ims,jms,lhc3), &
               e_bio(ims,jms,lhc5), e_bio(ims,jms,lhc8), &
               e_bio(ims,jms,lhcho), e_bio(ims,jms,liso), &
               e_bio(ims,jms,lket), e_bio(ims,jms,lno), &
               e_bio(ims,jms,loli), e_bio(ims,jms,lolt), &
               e_bio(ims,jms,lora1), e_bio(ims,jms,lora2), &
               e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl), &
               e_bio(ims,jms,lxyl), e_bio(ims,jms,lxyl) )
   end if
   if (config_flags%bio_emiss_opt /= GUNTHER1) then
      if (cbmz_addemiss_masscheck > 0) call addemiss_masscheck( &
               id, config_flags, 1, 'cbmz_addemiss_biobb', &
               dtstep, efact1, dz8w, chem, chem_sum, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte, &
               1, &
               e_iso,e_iso,e_iso,e_iso,e_iso,e_iso,e_iso, &
               e_iso,e_iso,e_iso,e_iso,e_iso,e_iso,e_iso, &
               e_iso,e_iso,e_iso,e_iso,e_iso,e_iso,e_iso )
      do j = jts, jte
      do k = kts, min(config_flags%kemit,kte)
      do i = its, ite
         chem(i,k,j,p_iso) = chem(i,k,j,p_iso) + e_iso(i,k,j) &
              *4.828e-4/rho_phy(i,k,j)*(dtstep/(dz8w(i,k,j)*60.))
      end do
      end do
      end do
      if (cbmz_addemiss_masscheck > 0) call addemiss_masscheck( &
               id, config_flags, 2, 'cbmz_addemiss_biobb', &
               dtstep, efact1, dz8w, chem, chem_sum, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte, &
               1, &
               e_iso,e_iso,e_iso,e_iso,e_iso,e_iso,e_iso, &
               e_iso,e_iso,e_iso,e_iso,e_iso,e_iso,e_iso, &
               e_iso,e_iso,e_iso,e_iso,e_iso,e_iso,e_iso )
   end if
   END subroutine cbmz_addemiss_bio
END MODULE module_cbmz_addemiss


text:
LOGICAL :: lopen
INTEGER :: lfnord
INQUIRE (UNIT=lfnord, OPENED=lopen)
END PROGRAM

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
!  This MODULE holds the routines that calculate air-sea exchange coefficients

MODULE module_sf_exchcoef
CONTAINS

  SUBROUTINE znot_m_v1(uref,znotm)
  IMPLICIT NONE

! uref(m/s)   :   Reference level wind
! znotm(meter):   Roughness scale for momentum
! Author      :  Biju Thomas on 02/07/2014
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znotm
    REAL             :: bs0, bs1, bs2, bs3, bs4, bs5, bs6
    REAL             :: cf0, cf1, cf2, cf3, cf4, cf5, cf6

    bs0 = -8.367276172397277e-12
    bs1 = 1.7398510865876079e-09
    bs2 = -1.331896578363359e-07
    bs3 = 4.507055294438727e-06
    bs4 = -6.508676881906914e-05
    bs5 = 0.00044745137674732834
    bs6 = -0.0010745704660847233

    cf0 = 2.1151080765239772e-13
    cf1 = -3.2260663894433345e-11
    cf2 = -3.329705958751961e-10
    cf3 = 1.7648562021709124e-07
    cf4 = 7.107636825694182e-06
    cf5 = -0.0013914681964973246
    cf6 = 0.0406766967657759

    IF ( uref .LE. 5.0 ) THEN
      znotm = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
    ELSEIF (uref .GT. 5.0 .AND. uref .LT. 10.0) THEN
      znotm =.00000235*(uref**2 - 25 ) + 3.805129199617346e-05
    ELSEIF ( uref .GE. 10.0  .AND. uref .LT. 60.0) THEN
      znotm = bs6 + bs5*uref + bs4*uref**2 + bs3*uref**3 + bs2*uref**4 +  &
              bs1*uref**5 + bs0*uref**6
    ELSE
      znotm = cf6 + cf5*uref + cf4*uref**2 + cf3*uref**3 + cf2*uref**4 +  &
              cf1*uref**5 + cf0*uref**6

    END IF

  END SUBROUTINE znot_m_v1

  SUBROUTINE znot_m_v0(uref,znotm)
  IMPLICIT NONE

! uref(m/s)   :   Reference level wind
! znotm(meter):   Roughness scale for momentum
! Author      :  Biju Thomas on 02/07/2014

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znotm
    REAL             :: yz, y1, y2, y3, y4

    yz =  0.0001344
    y1 =  3.015e-05
    y2 =  1.517e-06
    y3 = -3.567e-08
    y4 =  2.046e-10

    IF ( uref .LT. 12.5 ) THEN
       znotm  = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
    ELSE IF ( uref .GE. 12.5 .AND. uref .LT. 30.0 ) THEN
       znotm = (0.0739793 * uref -0.58)/1000.0
    ELSE
       znotm = yz + uref*y1 + uref**2*y2 + uref**3*y3 + uref**4*y4
    END IF

  END SUBROUTINE znot_m_v0

  SUBROUTINE znot_t_v1(uref,znott)
  IMPLICIT NONE

! uref(m/s)   :   Reference level wind
! znott(meter):   Roughness scale for temperature/moisture
! Author      :  Biju Thomas on 02/07/2014

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znott
    REAL             :: to0, to1, to2, to3
    REAL             :: tr0, tr1, tr2, tr3
    REAL             :: tn0, tn1, tn2, tn3, tn4, tn5
    REAL             :: ta0, ta1, ta2, ta3, ta4, ta5, ta6
    REAL             :: tt0, tt1, tt2, tt3, tt4, tt5, tt6, tt7

    tr0 = 6.451939325286488e-08
    tr1 = -7.306388137342143e-07
    tr2 = -1.3709065148333262e-05
    tr3 = 0.00019109962089098182

    to0 = 1.4379320027061375e-08
    to1 = -2.0674525898850674e-07
    to2 = -6.8950970846611e-06
    to3 = 0.00012199648268521026

    tn0 = 1.4023940955902878e-10
    tn1 = -1.4752557214976321e-08
    tn2 = 5.90998487691812e-07
    tn3 = -1.0920804077770066e-05
    tn4 = 8.898205876940546e-05
    tn5 = -0.00021123340439418298

    tt0 = 1.92409564131838e-12
    tt1 = -5.765467086754962e-10
    tt2 = 7.276979099726975e-08
    tt3 = -5.002261599293387e-06
    tt4 = 0.00020220445539973736
    tt5 = -0.0048088230565883
    tt6 = 0.0623468551971189
    tt7 = -0.34019193746967424

    ta0 = -1.7787470700719361e-10
    ta1 = 4.4691736529848764e-08
    ta2 = -3.0261975348463414e-06
    ta3 = -0.00011680322286017206
    ta4 = 0.024449377821884846
    ta5 = -1.1228628619105638
    ta6 = 17.358026773905973

    IF ( uref .LE. 7.0 ) THEN
      znott = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
    ELSEIF ( uref  .GE. 7.0 .AND. uref .LT. 12.5 ) THEN
      znott =  tr3 + tr2*uref + tr1*uref**2 + tr0*uref**3
    ELSEIF ( uref  .GE. 12.5 .AND. uref .LT. 15.0 ) THEN
      znott =  to3 + to2*uref + to1*uref**2 + to0*uref**3
    ELSEIF ( uref .GE. 15.0  .AND. uref .LT. 30.0) THEN
      znott =  tn5 + tn4*uref + tn3*uref**2 + tn2*uref**3 + tn1*uref**4 +   &
                                                       tn0*uref**5
    ELSEIF ( uref .GE. 30.0  .AND. uref .LT. 60.0) THEN
      znott = tt7 + tt6*uref + tt5*uref**2  + tt4*uref**3 + tt3*uref**4 +   &
             tt2*uref**5 + tt1*uref**6 + tt0*uref**7
    ELSE
      znott =  ta6 + ta5*uref + ta4*uref**2  + ta3*uref**3 + ta2*uref**4 +  &
              ta1*uref**5 + ta0*uref**6
    END IF

  END SUBROUTINE znot_t_v1

  SUBROUTINE znot_t_v0(uref,znott)
  IMPLICIT NONE

! uref(m/s)   :   Reference level wind
! znott(meter):   Roughness scale for temperature/moisture
! Author      :  Biju Thomas on 02/07/2014

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znott

    IF ( uref .LT. 7.0 ) THEN
       znott = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
    ELSE
       znott = (0.2375*exp(-0.5250*uref) + 0.0025*exp(-0.0211*uref))*0.01
    END IF

  END SUBROUTINE znot_t_v0

  SUBROUTINE znot_t_v2(uu,znott)
  IMPLICIT NONE

! uu in MKS
! znott in m
! Biju Thomas on 02/12/2015
!

    REAL, INTENT(IN) :: uu
    REAL, INTENT(OUT):: znott
    REAL             :: ta0, ta1, ta2, ta3, ta4, ta5, ta6
    REAL             :: tb0, tb1, tb2, tb3, tb4, tb5, tb6
    REAL             :: tt0, tt1, tt2, tt3, tt4, tt5, tt6

    ta0 = 2.51715926619e-09
    ta1 = -1.66917514012e-07
    ta2 = 4.57345863551e-06
    ta3 = -6.64883696932e-05
    ta4 = 0.00054390175125
    ta5 = -0.00239645231325
    ta6 = 0.00453024927761

    tb0 = -1.72935914649e-14
    tb1 = 2.50587455802e-12
    tb2 = -7.90109676541e-11
    tb3 = -4.40976353607e-09
    tb4 = 3.68968179733e-07
    tb5 = -9.43728336756e-06
    tb6 = 8.90731312383e-05

    tt0 = 4.68042680888e-14
    tt1 = -1.98125754931e-11
    tt2 = 3.41357133496e-09
    tt3 = -3.05130605309e-07
    tt4 = 1.48243563819e-05
    tt5 = -0.000367207751936
    tt6 = 0.00357204479347

    IF ( uu .LE. 7.0 ) THEN
       znott = (0.0185 / 9.8*(7.59e-4*uu**2+2.46e-2*uu)**2)
    ELSEIF ( uu  .GE. 7.0 .AND. uu .LT. 15. ) THEN
       znott = ta6 + ta5*uu + ta4*uu**2  + ta3*uu**3 + ta2*uu**4 +     &
               ta1*uu**5 + ta0*uu**6
    ELSEIF ( uu .GE. 15.0  .AND. uu .LT. 60.0) THEN
       znott = tb6 + tb5*uu + tb4*uu**2 + tb3*uu**3 + tb2*uu**4 +      &
               tb1*uu**5 + tb0*uu**6
    ELSE
       znott = tt6 + tt5*uu + tt4*uu**2  + tt3*uu**3 + tt2*uu**4 +    &
               tt1*uu**5 + tt0*uu**6
    END IF

  END SUBROUTINE znot_t_v2

  SUBROUTINE znot_m_v6(uref,znotm)
  IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
! For high winds, try to fit available observational data
!
! Bin Liu, NOAA/NCEP/EMC 2017
!
! uref(m/s)   :   wind speed at 10-m height
! znotm(meter):   areodynamical roughness scale over water
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znotm
    REAL             :: p13, p12, p11, p10
    REAL             :: p25, p24, p23, p22, p21, p20
    REAL             :: p35, p34, p33, p32, p31, p30
    REAL             :: p40

    p13 = -1.296521881682694e-02
    p12 =  2.855780863283819e-01
    p11 = -1.597898515251717e+00
    p10 = -8.396975715683501e+00

    p25 =  3.790846746036765e-10
    p24 =  3.281964357650687e-09
    p23 =  1.962282433562894e-07
    p22 = -1.240239171056262e-06
    p21 =  1.739759082358234e-07
    p20 =  2.147264020369413e-05

    p35 =  1.840430200185075e-07
    p34 = -2.793849676757154e-05
    p33 =  1.735308193700643e-03
    p32 = -6.139315534216305e-02
    p31 =  1.255457892775006e+00
    p30 = -1.663993561652530e+01

    p40 =  4.579369142033410e-04

    if (uref >= 0.0 .and.  uref <= 6.5 ) then
      znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
    elseif (uref > 6.5 .and. uref <= 15.7) then
      znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
    elseif (uref > 15.7 .and. uref <= 53.0) then
      znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30 )
    elseif ( uref > 53.0) then
      znotm = p40
    else
      print*, 'Wrong input uref value:',uref
    endif

  END SUBROUTINE znot_m_v6

  SUBROUTINE znot_t_v6(uref,znott)
  IMPLICIT NONE
! Calculate scalar roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
!
! Bin Liu, NOAA/NCEP/EMC 2017
!
! uref(m/s)   :   wind speed at 10-m height
! znott(meter):   scalar roughness scale over water
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znott

    REAL             :: p00
    REAL             :: p15, p14, p13, p12, p11, p10
    REAL             :: p25, p24, p23, p22, p21, p20
    REAL             :: p35, p34, p33, p32, p31, p30
    REAL             :: p45, p44, p43, p42, p41, p40
    REAL             :: p56, p55, p54, p53, p52, p51, p50
    REAL             :: p60

    p00 =  1.100000000000000e-04

    p15 = -9.144581627678278e-10
    p14 =  7.020346616456421e-08
    p13 = -2.155602086883837e-06
    p12 =  3.333848806567684e-05
    p11 = -2.628501274963990e-04
    p10 =  8.634221567969181e-04

    p25 = -8.654513012535990e-12
    p24 =  1.232380050058077e-09
    p23 = -6.837922749505057e-08
    p22 =  1.871407733439947e-06
    p21 = -2.552246987137160e-05
    p20 =  1.428968311457630e-04

    p35 =  3.207515102100162e-12
    p34 = -2.945761895342535e-10
    p33 =  8.788972147364181e-09
    p32 = -3.814457439412957e-08
    p31 = -2.448983648874671e-06
    p30 =  3.436721779020359e-05

    p45 = -3.530687797132211e-11
    p44 =  3.939867958963747e-09
    p43 = -1.227668406985956e-08
    p42 = -1.367469811838390e-05
    p41 =  5.988240863928883e-04
    p40 = -7.746288511324971e-03

    p56 = -1.187982453329086e-13
    p55 =  4.801984186231693e-11
    p54 = -8.049200462388188e-09
    p53 =  7.169872601310186e-07
    p52 = -3.581694433758150e-05
    p51 =  9.503919224192534e-04
    p50 = -1.036679430885215e-02

    p60 =  4.751256171799112e-05

    if (uref >= 0.0 .and. uref < 5.9 ) then
      znott = p00
    elseif (uref >= 5.9 .and. uref <= 15.4) then
      znott = p15*uref**5 + p14*uref**4 + p13*uref**3 + p12*uref**2 + p11*uref + p10
    elseif (uref > 15.4 .and. uref <= 21.6) then
      znott = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
    elseif (uref > 21.6 .and. uref <= 42.2) then
      znott = p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30
    elseif ( uref > 42.2 .and. uref <= 53.3) then
      znott = p45*uref**5 + p44*uref**4 + p43*uref**3 + p42*uref**2 + p41*uref + p40
    elseif ( uref > 53.3 .and. uref <= 80.0) then
      znott = p56*uref**6 + p55*uref**5 + p54*uref**4 + p53*uref**3 + p52*uref**2 + p51*uref + p50
    elseif ( uref > 80.0) then
      znott = p60
    else
      print*, 'Wrong input uref value:',uref
    endif

  END SUBROUTINE znot_t_v6

  SUBROUTINE znot_m_v7(uref,znotm)
  IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
! For high winds, try to fit available observational data
! Comparing to znot_t_v6, slightly decrease Cd for higher wind speed
!
! Bin Liu, NOAA/NCEP/EMC 2018
!
! uref(m/s)   :   wind speed at 10-m height
! znotm(meter):   areodynamical roughness scale over water
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znotm
    REAL             :: p13, p12, p11, p10
    REAL             :: p25, p24, p23, p22, p21, p20
    REAL             :: p35, p34, p33, p32, p31, p30
    REAL             :: p40

    p13 = -1.296521881682694e-02
    p12 =  2.855780863283819e-01
    p11 = -1.597898515251717e+00
    p10 = -8.396975715683501e+00

    p25 =  3.790846746036765e-10
    p24 =  3.281964357650687e-09
    p23 =  1.962282433562894e-07
    p22 = -1.240239171056262e-06
    p21 =  1.739759082358234e-07
    p20 =  2.147264020369413e-05

    p35 =  1.897534489606422e-07
    p34 = -3.019495980684978e-05
    p33 =  1.931392924987349e-03
    p32 = -6.797293095862357e-02
    p31 =  1.346757797103756e+00
    p30 = -1.707846930193362e+01

    p40 =  3.371427455376717e-04

    if (uref >= 0.0 .and.  uref <= 6.5 ) then
      znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
    elseif (uref > 6.5 .and. uref <= 15.7) then
      znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
    elseif (uref > 15.7 .and. uref <= 53.0) then
      znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30 )
    elseif ( uref > 53.0) then
      znotm = p40
    else
      print*, 'Wrong input uref value:',uref
    endif

  END SUBROUTINE znot_m_v7

  SUBROUTINE znot_t_v7(uref,znott)
  IMPLICIT NONE
! Calculate scalar roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
! To be compatible with the slightly decreased Cd for higher wind speed
!
! Bin Liu, NOAA/NCEP/EMC 2018
!
! uref(m/s)   :   wind speed at 10-m height
! znott(meter):   scalar roughness scale over water
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znott

    REAL             :: p00
    REAL             :: p15, p14, p13, p12, p11, p10
    REAL             :: p25, p24, p23, p22, p21, p20
    REAL             :: p35, p34, p33, p32, p31, p30
    REAL             :: p45, p44, p43, p42, p41, p40
    REAL             :: p56, p55, p54, p53, p52, p51, p50
    REAL             :: p60

    p00 =  1.100000000000000e-04

    p15 = -9.193764479895316e-10
    p14 =  7.052217518653943e-08
    p13 = -2.163419217747114e-06
    p12 =  3.342963077911962e-05
    p11 = -2.633566691328004e-04
    p10 =  8.644979973037803e-04

    p25 = -9.402722450219142e-12
    p24 =  1.325396583616614e-09
    p23 = -7.299148051141852e-08
    p22 =  1.982901461144764e-06
    p21 = -2.680293455916390e-05
    p20 =  1.484341646128200e-04

    p35 =  7.921446674311864e-12
    p34 = -1.019028029546602e-09
    p33 =  5.251986927351103e-08
    p32 = -1.337841892062716e-06
    p31 =  1.659454106237737e-05
    p30 = -7.558911792344770e-05

    p45 = -2.694370426850801e-10
    p44 =  5.817362913967911e-08
    p43 = -5.000813324746342e-06
    p42 =  2.143803523428029e-04
    p41 = -4.588070983722060e-03
    p40 =  3.924356617245624e-02

    p56 = -1.663918773476178e-13
    p55 =  6.724854483077447e-11
    p54 = -1.127030176632823e-08
    p53 =  1.003683177025925e-06
    p52 = -5.012618091180904e-05
    p51 =  1.329762020689302e-03
    p50 = -1.450062148367566e-02

    p60 =  6.840803042788488e-05

    if (uref >= 0.0 .and. uref < 5.9 ) then
      znott = p00
    elseif (uref >= 5.9 .and. uref <= 15.4) then
      znott = p15*uref**5 + p14*uref**4 + p13*uref**3 + p12*uref**2 + p11*uref + p10
    elseif (uref > 15.4 .and. uref <= 21.6) then
      znott = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
    elseif (uref > 21.6 .and. uref <= 42.6) then
      znott = p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30
    elseif ( uref > 42.6 .and. uref <= 53.0) then
      znott = p45*uref**5 + p44*uref**4 + p43*uref**3 + p42*uref**2 + p41*uref + p40
    elseif ( uref > 53.0 .and. uref <= 80.0) then
      znott = p56*uref**6 + p55*uref**5 + p54*uref**4 + p53*uref**3 + p52*uref**2 + p51*uref + p50
    elseif ( uref > 80.0) then
      znott = p60
    else
      print*, 'Wrong input uref value:',uref
    endif

  END SUBROUTINE znot_t_v7

  SUBROUTINE znot_m_v8(uref,znotm)
  IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
! For high winds, try to fit available observational data
! Comparing to znot_t_v6, slightly decrease Cd for higher wind speed
! And this is another variation similar to v7
!
! Bin Liu, NOAA/NCEP/EMC 2018
!
! uref(m/s)   :   wind speed at 10-m height
! znotm(meter):   areodynamical roughness scale over water
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znotm
    REAL             :: p13, p12, p11, p10
    REAL             :: p25, p24, p23, p22, p21, p20
    REAL             :: p35, p34, p33, p32, p31, p30
    REAL             :: p40

    p13 = -1.296521881682694e-02
    p12 =  2.855780863283819e-01
    p11 = -1.597898515251717e+00
    p10 = -8.396975715683501e+00

    p25 =  3.790846746036765e-10
    p24 =  3.281964357650687e-09
    p23 =  1.962282433562894e-07
    p22 = -1.240239171056262e-06
    p21 =  1.739759082358234e-07
    p20 =  2.147264020369413e-05

    p35 =  1.897534489606422e-07
    p34 = -3.019495980684978e-05
    p33 =  1.931392924987349e-03
    p32 = -6.797293095862357e-02
    p31 =  1.346757797103756e+00
    p30 = -1.707846930193362e+01

    p40 =  3.886804744928044e-04

    if (uref >= 0.0 .and.  uref <= 6.5 ) then
      znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
    elseif (uref > 6.5 .and. uref <= 15.7) then
      znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
    elseif (uref > 15.7 .and. uref <= 51.5) then
      znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30 )
    elseif ( uref > 51.5) then
      znotm = p40
    else
      print*, 'Wrong input uref value:',uref
    endif

  END SUBROUTINE znot_m_v8

  SUBROUTINE znot_t_v8(uref,znott)
  IMPLICIT NONE
! Calculate scalar roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
! To be compatible with the slightly decreased Cd for higher wind speed
! And this is another variation similar to v7
!
! Bin Liu, NOAA/NCEP/EMC 2018
!
! uref(m/s)   :   wind speed at 10-m height
! znott(meter):   scalar roughness scale over water
!

    REAL, INTENT(IN) :: uref
    REAL, INTENT(OUT):: znott

    REAL             :: p00
    REAL             :: p15, p14, p13, p12, p11, p10
    REAL             :: p25, p24, p23, p22, p21, p20
    REAL             :: p35, p34, p33, p32, p31, p30
    REAL             :: p45, p44, p43, p42, p41, p40
    REAL             :: p56, p55, p54, p53, p52, p51, p50
    REAL             :: p60

    p00 =  1.100000000000000e-04

    p15 = -9.193764479895316e-10
    p14 =  7.052217518653943e-08
    p13 = -2.163419217747114e-06
    p12 =  3.342963077911962e-05
    p11 = -2.633566691328004e-04
    p10 =  8.644979973037803e-04

    p25 = -9.402722450219142e-12
    p24 =  1.325396583616614e-09
    p23 = -7.299148051141852e-08
    p22 =  1.982901461144764e-06
    p21 = -2.680293455916390e-05
    p20 =  1.484341646128200e-04

    p35 =  7.921446674311864e-12
    p34 = -1.019028029546602e-09
    p33 =  5.251986927351103e-08
    p32 = -1.337841892062716e-06
    p31 =  1.659454106237737e-05
    p30 = -7.558911792344770e-05

    p45 = -2.706461188613193e-10
    p44 =  5.845859022891930e-08
    p43 = -5.027577045502003e-06
    p42 =  2.156326523752734e-04
    p41 = -4.617267288861201e-03
    p40 =  3.951492707214883e-02

    p56 = -1.112896580069263e-13
    p55 =  4.450334755105140e-11
    p54 = -7.375373918500171e-09
    p53 =  6.493685149526543e-07
    p52 = -3.206421106713471e-05
    p51 =  8.407596231678149e-04
    p50 = -9.027924333673693e-03

    p60 =  5.791179079892191e-05

    if (uref >= 0.0 .and. uref < 5.9 ) then
      znott = p00
    elseif (uref >= 5.9 .and. uref <= 15.4) then
      znott = p15*uref**5 + p14*uref**4 + p13*uref**3 + p12*uref**2 + p11*uref + p10
    elseif (uref > 15.4 .and. uref <= 21.6) then
      znott = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
    elseif (uref > 21.6 .and. uref <= 42.6) then
      znott = p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30
    elseif ( uref > 42.6 .and. uref <= 51.5) then
      znott = p45*uref**5 + p44*uref**4 + p43*uref**3 + p42*uref**2 + p41*uref + p40
    elseif ( uref > 51.5 .and. uref <= 80.0) then
      znott = p56*uref**6 + p55*uref**5 + p54*uref**4 + p53*uref**3 + p52*uref**2 + p51*uref + p50
    elseif ( uref > 80.0) then
      znott = p60
    else
      print*, 'Wrong input uref value:',uref
    endif

  END SUBROUTINE znot_t_v8

   SUBROUTINE znot_wind10m(w10m,znott,znotm,icoef_sf)
   IMPLICIT NONE

! w10m(m/s)   :   10-m wind speed
! znott(meter):   Roughness scale for temperature/moisture, zt
! znotm(meter):   Roughness scale for momentum, z0
! Author      :  Weiguo Wang on 02/24/2016
!            convert from icoef=0,1,2 to have 10m level cd, ch match obs
     REAL, INTENT(IN) :: w10m
     INTEGER, INTENT(IN) :: icoef_sf
     REAL, INTENT(OUT):: znott, znotm

     real :: zm,zt,windmks, zlev,z10, tmp, zlevt, aaa, zm1,zt1
        zlev=20.0
        zlevt=10.0
        z10=10.0
            windmks=w10m
            if (windmks > 85.0) windmks=85.0
            if (windmks < 1.0) windmks=1.0
            if ( icoef_sf .EQ. 1) then
              call  znot_m_v1(windmks,zm1)
              call  znot_t_v1(windmks,zt1)

            else if ( icoef_sf .EQ. 0 ) then
              call  znot_m_v0(windmks,zm1)
              call  znot_t_v0(windmks,zt1)

            else  if( icoef_sf .EQ. 2 ) then
              call  znot_m_v1(windmks,zm1)
              call  znot_t_v2(windmks,zt1)

            else  if( icoef_sf .EQ. 3 ) then
              call  znot_m_v1(windmks,zm)
              call  znot_t_v2(windmks,zt)
!! adjust a little to match obs at 10m, cd is reduced
            tmp=0.4*0.4/(alog(zlev/zm))**2   ! cd at zlev
            zm1=z10/exp( sqrt(0.4*0.4/(tmp*0.95-0.0002)) )
!ch
            tmp=0.4*0.4/(alog(zlevt/zm)*alog(zlevt/zt))  ! ch at zlev using old formula
            zt1=z10/exp( 0.4*0.4/( 0.95*tmp*alog(z10/zm1) )  )

            else if( icoef_sf .EQ. 4 ) then

              call  znot_m_v1(windmks,zm)
              call  znot_t_v2(windmks,zt)
!!  for wind<20, cd similar to icoef=2 at 10m, then reduced
             tmp=0.4*0.4/(alog(10.0/zm))**2   ! cd at zlev
             aaa=0.75
            if (windmks < 20) then
              aaa=0.99
            elseif(windmks < 45.0) then
              aaa=0.99+(windmks-20)*(0.75-0.99)/(45.0-20.0)
            endif
            zm1=z10/exp( sqrt(0.4*0.4/(tmp*aaa)) )
!ch
          tmp=0.4*0.4/(alog(zlevt/zm)*alog(zlevt/zt))  ! ch at zlev using old formula
            zt1=z10/exp( 0.4*0.4/( 0.95*tmp*alog(z10/zm1) )  )

            else if( icoef_sf .EQ. 5 ) then

              call  znot_m_v1(windmks,zm)
              call  znot_t_v2(windmks,zt)
!!  for wind<20, cd similar to icoef=2 at 10m, then reduced
             tmp=0.4*0.4/(alog(10.0/zm))**2   ! cd at zlev
             aaa=0.80
            if (windmks < 20) then
              aaa=1.0
            elseif(windmks < 45.0) then
              aaa=1.0+(windmks-20)*(0.80-1.0)/(45.0-20.0)
            endif
            zm1=z10/exp( sqrt(0.4*0.4/(tmp*aaa)) )
!ch
          tmp=0.4*0.4/(alog(zlevt/zm)*alog(zlevt/zt))  ! ch at zlev using old formula
            zt1=z10/exp( 0.4*0.4/( 1.0*tmp*alog(z10/zm1) )  )

            else  if( icoef_sf .EQ. 6 ) then
              call  znot_m_v6(windmks,zm1)
              call  znot_t_v6(windmks,zt1)
            else  if( icoef_sf .EQ. 7 ) then
              call  znot_m_v7(windmks,zm1)
              call  znot_t_v7(windmks,zt1)
            else  if( icoef_sf .EQ. 8 ) then
              call  znot_m_v8(windmks,zm1)
              call  znot_t_v8(windmks,zt1)
           else
             write(0,*)'stop, icoef_sf must be one of 0,1,2,3,4,5,6,7,8'
             stop
          endif
          znott=zt1
          znotm=zm1

  end subroutine znot_wind10m

END MODULE module_sf_exchcoef

max_stars_repo_name: wramstrom/ccpp-physics

max_issues_repo_name: wramstrom/ccpp-physics

max_forks_repo_name: wramstrom/ccpp-physics