cv3_inip.F

Go to the documentation of this file.

        SUBROUTINE cv3_inip()
***************************************************************
*                                                             *
* cv3_inip lecture des choix de lois de probabilité de mélange*
*          et calcul de leurs coefficients normalisés.        *
*                                                             *
* written by   : jean-yves grandpeix, 06/06/2006, 19.39.27    *
* modified by :                                               *
***************************************************************
*
#include "YOMCST2.h"
c
c      INTEGER iflag_mix
      include 'iniprint.h'

      CHARACTER (LEN=20) :: modname='cv3_inip'
      CHARACTER (LEN=80) :: abort_message

c
c --   mixing probability distribution functions
c
      real qcoef1,qcoef2,qff,qfff,qmix,rmix,qmix1,rmix1,qmix2,rmix2,f
      qcoef1(f) = tanh(f/gammas)
      qcoef2(f) = ( tanh(f/gammas) + gammas *
     $            log(cosh((1.- f)/gammas)/cosh(f/gammas)))
      qff(f) = max(min(f,1.),0.)
      qfff(f) = min(qff(f),scut)
      qmix1(f) = ( tanh((qff(f) - fmax)/gammas)+qcoef1max )/
     $           qcoef2max
      rmix1(f) = ( gammas*log(cosh((qff(f)-fmax)/gammas))
     1             +qff(f)*qcoef1max ) / qcoef2max
      qmix2(f) = -log(1.-qfff(f))/scut
      rmix2(f) = (qfff(f)+(1.-qff(f))*log(1.-qfff(f)))/scut
      qmix(f) = qqa1*qmix1(f) + qqa2*qmix2(f)
      rmix(f) = qqa1*rmix1(f) + qqa2*rmix2(f)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c===========================================================================
c               READ in parameters for the mixing distribution
c               and pass these through a COMMON block to SUBROUTINE convect etc.
c(written by v.t.j. phillips, 20-30/jan/99)
c===========================================================================
c
c   line 1:  a flag(0 or 1) to decide whether p(f) = 1 or the general p(f) is to be
c         used, followed by scut, which is the cut-off value of f in convect
c   line 2:  blank
c   line 3:  the coefficients for the linear combination of p(f)s to
c                 make the general p(f)
c   line 4:  blank
c   line 5:  gammas, fmax for the cosh^2 component of p(f)
c   line 6:  blank
c   line 7:  alphas for the 1st irrational p(f)
c   line 8:  blank
c   line 9:  betas  for the 2nd irrational p(f)
c

cc$$$        open(57,file='parameter_mix.data')
cc$$$
cc$$$        read(57,*) iflag_mix, scut
cc$$$        read(57,*)
cc$$$        if(iflag_mix .gt. 0) then
cc$$$                 read(57,*) qqa1, qqa2
cc$$$              read(57,*)
cc$$$              read(57,*) gammas, fmax
cc$$$              read(57,*)
cc$$$              read(57,*) alphas
cc$$$         endif
cc$$$            close(57)

c
      if(iflag_mix .gt. 0) then
c
c--      normalize pdf weights
c
        sumcoef=qqa1+qqa2
        qqa1=qqa1/sumcoef
        qqa2=qqa2/sumcoef
c
        qcoef1max = qcoef1(fmax)
        qcoef2max = qcoef2(fmax)
c
        sigma = 0.
        aire=0.0
        pdf=0.0
        mu=0.0
        df = 0.0001
c
c        do ff = 0.0 + df, 1.0 - 2.*df, df
         ff=df
         dowhile( ff .le. 1.0 - 2.*df )
              pdf = (qmix(ff+df) -  qmix(ff)) * (1.-ff) / df
              aire=aire+(qmix(ff+df) - qmix(ff)) * (1.-ff)
              mu = mu + pdf * ff * df
         IF(prt_level>9)WRITE(lunout,*)                                 &
     &               pdf,  qmix(ff), aire, ff
         ff=ff+df
         enddo
c
c         do ff=0.0+df,1.0 - 2.*df,df
          ff=df
          dowhile( ff .le. 1.0 - 2.*df )
              pdf = (qmix(ff+df)- qmix(ff)) * (1.-ff) / df
              sigma = sigma+pdf*(ff - mu)*(ff - mu)*df
         ff=ff+df
         enddo
         sigma = sqrt(sigma)
c
        if (abs(aire-1.0) .gt. 0.02) then
            write(lunout,*)'WARNING:: AREA OF MIXING PDF IS::', aire
            abort_message = ''
            CALL abort_gcm(modname,abort_message,1)
        else
            print *,'Area, mean & std deviation are ::', aire,mu,sigma
        endif
      endif     !  (iflag_mix .gt. 0)

      RETURN
      END