!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&      
      program constDChs
!***********************************************************************
!                                                                      *
!         Define commons /dynam/ and /ctlblk/ for dry core version     *
!         /held/ block for held-suarez experiment                      *
!***********************************************************************
      implicit none
      include "param_o.h"
      include "const.h"
!-----------------------------------------------------------------------
!
! *** Universal constants.
!
      real, parameter :: g=9.80616, r=287.04, cp=1004.6, rcp = 1./cp,  &
                         kappa=r*rcp, gamma=0.0065, prf0=101325.,      &
                         om=0.00007292, t0=288.,re=6.37122e6,          &
                         kf=1./(24.*3600.),sigmab=0.7,ka=1/40.*kf,     &
                         ks=0.25*kf,delty=60,delthz=10,p0=1e5
!-----------------------------------------------------------------------
! For common /dynam/ 
!-----------------------------------------------------------------------
      real :: fadv, fadt, rd, f4d, ef4t, fkin, fcp, qmin, rtdf, fadv2  
      real, dimension (lm) :: deta, rdeta, aeta, daeta,f4q2
      real, dimension (lp1) :: eta, dfl
      real, dimension (0:im+1,0:jm+1,nm) :: wpdar, f11, f12, f21, f22, &
                                    p11, p12, p21, p22, fdiv, fddmp
!-----------------------------------------------------------------------
      real, dimension(:,:,:), allocatable :: sqv, sqh, rsqv, rsqh, &
             dxvdx, dxvdy, dyvdx, dyvdy, dzvdx, dzvdy, q11, q12, q22, &
             vcosp, vsinp, vlam, vphi, hcosp, hsinp, hlam, hphi,fvdiff &
             ,hbmsk
!-----------------------------------------------------------------------
      integer, dimension(im,jm):: icorner
!-----------------------------------------------------------------------
!!!      character(len=34)::  Cgrid
      character(len=150)::  Cgrid
      integer::list10,n,i,j,l,list21,idtadv,INDAY,IOUTNUM,IIOUT,IHRST
      real::dst,rdst,cmi,gorg,alpt,omdt,fpgf,fcor05,rq,omg,a110,a120,a210,&
           a220,rgog
!
! held suarez constant
!
      real,dimension(0:im+1,0:jm+1,nm,lm)::ktdt,teq
      real,dimension(lm)::kvdt
      real::sgmsgb,pres

!-----------------------------------------------------------------------
      allocate(sqv     (0:im+1,0:jm+1,nm))
      allocate(sqh     (0:im+1,0:jm+1,nm))
      allocate(rsqv    (0:im+1,0:jm+1,nm))
      allocate(rsqh    (0:im+1,0:jm+1,nm))
      allocate(dxvdx   (0:im+1,0:jm+1,nm))
      allocate(dxvdy   (0:im+1,0:jm+1,nm))
      allocate(dyvdx   (0:im+1,0:jm+1,nm))
      allocate(dyvdy   (0:im+1,0:jm+1,nm))
      allocate(dzvdx   (0:im+1,0:jm+1,nm))
      allocate(dzvdy   (0:im+1,0:jm+1,nm))
      allocate(q11     (0:im+1,0:jm+1,nm))
      allocate(q12     (0:im+1,0:jm+1,nm))
      allocate(q22     (0:im+1,0:jm+1,nm))
      allocate(vcosp   (0:im+1,0:jm+1,nm))
      allocate(vsinp   (0:im+1,0:jm+1,nm))
      allocate(vlam    (0:im+1,0:jm+1,nm))
      allocate(vphi    (0:im+1,0:jm+1,nm))
      allocate(hcosp   (0:im+1,0:jm+1,nm))
      allocate(hsinp   (0:im+1,0:jm+1,nm))
      allocate(hlam    (0:im+1,0:jm+1,nm))
      allocate(hphi    (0:im+1,0:jm+1,nm))
      allocate(fvdiff  (0:im+1,0:jm+1,nm))
      allocate(hbmsk   (0:im+1,0:jm+1,nm))

!GSM
      NAMELIST /CNSTDATA/INDAY,IHRST,IDAT,IOUTNUM,IIOUT

      open(19,file='cnstdata',status='old')
      READ(19,CNSTDATA)
!GSM

!-----------------------------------------------------------------------
      print *,iiout,outend,outstr,ioutnum
!-----------------------------------------------------------------------

!
! Read grid elements grid metrcis
!
!ctest
      Cgrid='gridinit.dat'
      print *,'Cgrid=', Cgrid
!ctest

      list10=10
      Open(unit=list10,file=Cgrid,form="unformatted")
      read (list10) sqv, sqh, q11, q12, q22, & 
                    dxvdx, dxvdy, dyvdx, dyvdy, dzvdx, dzvdy, &
                    rsqv, vcosp, vsinp, vlam, vphi, &
                    rsqh, hcosp, hsinp, hlam, hphi
      close(list10)
!      print *,hlam(123,41,5)*180/3.1415926,hphi(123,41,5)*180/3.1415926
!      print *,hlam(124,41,5)*180/3.1415926,hphi(124,41,5)*180/3.1415926
!      print *,hlam(123,42,5)*180/3.1415926,hphi(123,42,5)*180/3.1415926
!      print *,hlam(124,42,5)*180/3.1415926,hphi(124,42,5)*180/3.1415926
!      stop
!-----------------------------------------------------------------------
!
! *** Distance
!
      dst = 2.*Re    
!!!!      print *,"lm=",lm,Re,dst
     
      rdst = 1./dst

      cmi = dt/dst
!
! *** Top of models atmosphere
!
      rd = r

!
! *** Advection constants
!
      idtadv = idtad
      
      fadv = cmi * idtadv /(3.*4.*4.)
      fadv2 = -cmi * idtadv /(3.)
      fadt = - cmi * idtadv /(3.*2.*4.)
      fkin = - cmi /(2.*4.)
      fcp = rcp * 0.125 * cmi * (1./3.)
      f4d = - idtadv /(2.)
      ef4t = rcp /(2.)

      if(nm.eq.6)then
!        rtdf=1.81e-12*(im-1)**4
!        rtdf=0.4e-12*(im-1)**4
         rtdf=1.9e-6*16
!         rtdf=1.9e-6*16*3
      else if(nm.eq.14) then
!        rtdf=8.4e-12*(im-1)**4
!        rtdf=2.1e-12*(im-1)**4
!        rtdf=2.1e-13*(im-1)**4
         rtdf=1.9e-6*16
!         rtdf=1.9e-6*16*3.
!         rtdf=1.9e-6*16*30.
      endif

!
! *** Vertical discretization constants
!
      gorg = g/(r*gamma)
      rgog = 1./gorg

        call getDeta( deta )
        rdeta = 1./deta
        f4q2=-0.25*dt*idtadv*rdeta

          eta(1) = 0.
        do l=1,lm
          eta(l+1) = eta(l) + deta(l)
!          print *,l,eta(l),eta(l+1)
        end do
          eta(lp1) = 1.   
  
        do l=1,lm
          aeta(l) = 0.5*(eta(l) + eta(l+1))
        end do 

!!        print *,(aeta(l),l=lm,1,-1)

      do l=1,lp1
         dfl(l)=t0*(1.-((pt+eta(l)*(prf0-pt))/prf0)**rgog)/gamma
      end do

        daeta(1) = kappa * 0.5 * deta(1)
      do l=2,lm
        daeta(l) = kappa * 0.5 * (deta(l) + deta(l-1))
      end do

      call sfcEta(eta,aeta,dfl,im,jm,lm,nm,flat)
!
! *** Hozizontal discretization constants
!
        icorner(:,:)=1
      
        icorner(1 ,1 )=0
        icorner(1,jm )=0
        icorner(im,jm)=0
        icorner(im,1 )=0


        do n=1,nm
        do i=1,im
        do j=1,jm
           wpdar(i,j,n) =-0.5 *cmi*cmi*weig*rsqh(i,j,n)*7.84e4
           fdiv(i,j,n) = 0.25 *cmi * rsqh(i,j,n)
!           print *,i,j,n,fdiv(i,j,n)
        end do
        end do
        end do

        if(sigma)then
          wpdar=0.5*wpdar
        endif

      omdt = om * dt
      fpgf = -cmi/6.
        
      do n=1,nm
        do i=1,im1
        do j=1,jm1
          fcor05=omdt*vsinp(i,j,n)
          rQ = 1./(1.+fcor05**2)
          omg = fcor05 * sqv(i,j,n)
            a110 = (1. + omg*q12(i,j,n)) * rQ
            a120 =       omg*q22(i,j,n)  * rQ
            a210 =     - omg*q11(i,j,n)  * rQ
            a220 =(1.  - omg*q12(i,j,n)) * rQ
            f11(i,j,n) = 2. * a110 - 1.
            f12(i,j,n) = 2. * a120
            f21(i,j,n) = 2. * a210
            f22(i,j,n) = 2. * a220 - 1.
            p11(i,j,n) = a110 *fpgf
            p12(i,j,n) = a120 *fpgf
            p21(i,j,n) = a210 *fpgf
            p22(i,j,n) = a220 *fpgf
        end do
        end do
        end do
!----------------------------------------------------------------------
      fddmp(:,:,:) = 0.

!
!  calculate held suarez constants
!

       do L=1,lm
         sgmsgb=(aeta(l)-sigmab)/(1.-sigmab)
         if(sgmsgb.gt.0)then
           kvdt(l)=kf*dt*sgmsgb
         else
           kvdt(l)=0.
         endif
!!         print *,l,aeta(l),sgmsgb,kvdt(l)

         do n=1,nm
         do i=0,im+1
         do j=0,jm+1
           if(sgmsgb.gt.0)then
             ktdt(i,j,n,l)=(ka+(ks-ka)*sgmsgb*hcosp(i,j,n)**4)*dt
           else
             ktdt(i,j,n,l)=ka*dt
           endif
         enddo
         enddo
         enddo

       enddo

!  velocity diffusion constants

      fvdiff=0  !dt*2.e-5
!      fvdiff=dt*2.0e-4

      do n=1,nm
      do i=1,im
      do j=1,jm
        if (sqv(i,j,n)==0) then
!         write(0,*)'sqv(i,j,n)==0.0', i,j,n,sqv(i,j,n)
        fvdiff(i,j,n)=0.0
        else
        fvdiff(i,j,n)=fvdiff(i,j,n)*sqv(1,1,1)/sqv(i,j,n)
        end if
      enddo
      enddo
      enddo

!  boundary condition mask
    
      call init_hbmsk(hbmsk,im,jm,nm)
      call cucvparm(pt)
      
!----------------------------------------------------------------------
!!!       print *,im,jm,nm
       dfl=dfl*g
       list21 = 21
       open(list21, file='dynam.dat',form ="unformatted")
       write(list21) fadv, fadv2,fadt, rd,  f4d,f4q2, ef4t, fkin, &
                  deta, rdeta, aeta, eta, dfl, daeta, &
                  wpdar, f11, f12, f21, f22, &
                  p11, p12, p21, p22, fcp, fdiv, fddmp,hsinp,hcosp,&
                  delty,delthz,kappa,p0,rtdf,fvdiff,hbmsk
       close(list21)

!!!       print *,fadv,fadt,rd,f4d,f4q2,ef4t,fkin
! 
! define control block values
!
       print *,'total time:',ntstm*dt/(3600*24)
       print *,'date:',idat
       run = .true.
       first = .true.
!!!       restrt = .false.   !commented 04/02/2019 when introduced new restart option for the model, Dragan
       nfcst = 15
       nbc = 20
       list = 6
       ntsd = 0
       nddamp = 0
       nprec = 0
       nboco = 0
       nshde = 0
       ncp = 0
       subpost = .false.
       ntddmp = 0
       open(list21,file='ctlblk.dat',form="unformatted")
       write(list21)run,idat,ihrst,first,restrt,nfcst,nbc,list, &
        ntsd,nddamp,nprec,nboco,nshde,ncp, &
        subpost,ntddmp
       close(list21)

       list21 = 21
       open(list21, file='hs.dat',form ="unformatted")
       write(list21)kvdt,ktdt
       close(list21)

       open(list21, file='dgns.dat',form ="unformatted")
       write(list21)dgnnum,dgnstr,idgns,ickmm
       close(list21)


       print *,'ntstm:',ntstm

       print *, 'FINISH WRITING'
!----------------------------------------------------------------------
       end program constDChs