var/da/da_radiance/da_cloud_sim.inc

   1 subroutine da_cloud_sim(KINDIC,KDIM,PX,PF,PG,IZS,RZS,DZS)
   2
   3 ! Purpose :
   4 ! -------
   5 ! Simulate the cloud as a linear combination of grey clouds at model levels
   6
   7 ! Interface :
   8 ! ---------
   9 ! KINDIC
  10 ! KDIM           : Dimension of cloud fraction variable
  11 ! PX             : Cloud fraction variable             -> Input
  12 ! PF             : Fitness Error
  13 ! PG             : Gradient of cloud fraction variable -> Output
  14
  15 ! Externals :
  16 ! ---------
  17
  18 ! Method :
  19 ! ------
  20
  21 ! Reference :
  22 ! ---------
  23
  24 ! Author :
  25 ! ------
  26 ! 01/08/2005  Thomas Auligne         *ECMWF*
  27
  28 ! Modifications :
  29 ! -------------
  30
  31 ! ---------------------------------------------------------------------------------------------
  32
  33 IMPLICIT NONE
  34
  35 !! Parameters !!
  36 INTEGER,INTENT(IN)      :: KINDIC
  37 INTEGER,INTENT(IN)      :: KDIM
  38 INTEGER,INTENT(IN)      :: IZS(2)
  39 double precision   ,INTENT(INOUT)   :: PX(KDIM) !izs(2)
  40 double precision   ,INTENT(OUT)     :: PF
  41 double precision   ,INTENT(OUT)     :: PG(KDIM)
  42 real               ,INTENT(IN)      :: RZS(kdim*izs(2))      ! Eigenvectors
  43 DOUBLE PRECISION   ,INTENT(IN)      :: DZS(IZS(1)*KDIM) ! AMAT
  44 !! Local arrays !!
  45 INTEGER                 :: JCH, ilev, JLEV, nchan, neignvec
  46 REAL                    :: ZNORM_PG, ZCLR, ZDCLR, eignvec(kdim,izs(2)), eignval(izs(2))
  47 double precision        :: AMAT(IZS(1),KDIM)
  48 double precision        :: alpha, beta
  49 double precision        :: zx(KDIM), zgx(KDIM, KDIM), zx_eof(KDIM)
  50
  51 !IF (KINDIC == 1) RETURN
  52  PF       = 0.0
  53  PG       = 0.0
  54  nchan    = izs(1)
  55  neignvec = izs(2)
  56  !eignvec  = RESHAPE(rzs(1:KDIM*neignvec),(/KDIM,neignvec/))
  57  !eignval  = rzs(KDIM*neignvec+1:(KDIM+1)*neignvec)
  58
  59  AMAT     = RESHAPE(DZS(1:NCHAN*KDIM),(/NCHAN,KDIM/))
  60  PX(KDIM) = 1.0 - SUM(PX(1:kdim-1))
  61 ! where (PX < 0.0) PX = 0.0
  62 ! where (PX > 1.0) PX = 1.0
  63
  64  !ZX_EOF   = MATMUL(eignvec,eignval*PX)
  65 !!! ZX_EOF   = MATMUL(eignvec,MATMUL(TRANSPOSE(eignvec),PX))
  66  zx_eof = PX
  67  zx     = zx_eof
  68
  69  ! Softmax (= multiple-logistic) variable transform
  70  !beta = 100.0
  71
  72  !zx = exp(beta*zx_eof) / SUM(exp(beta*zx_eof))
  73  !do ilev = 1, kdim
  74  !   do jlev = 1, kdim
  75  !      zgx(ilev,jlev) = - beta * zx(ilev) * zx(jlev)
  76  !      if (ilev == jlev) zgx(ilev,jlev) =  zgx(ilev,jlev) + zx(ilev) * beta
  77  !   end do
  78  !end do
  79
  80  DO JCH=1,NCHAN
  81    PF = PF + 0.5 * (SUM(ZX*AMAT(JCH,:)) - 1.0)**2
  82    DO JLEV=1,KDIM
  83       PG(JLEV) = PG(JLEV) + (AMAT(JCH,JLEV)-AMAT(JCH,KDIM)) * (SUM(ZX*AMAT(JCH,:)) - 1.0)
  84    ENDDO
  85  ENDDO
  86  !PG = MATMUL(PG, zgx)
  87
  88  alpha = float(nchan)*100.0
  89  PF = PF + 0.5*alpha*SUM(ZX**2, MASK=ZX<0.0)
  90  WHERE (ZX<0.0) PG = PG + alpha*ZX
  91
  92 !write(*,'(a,2f10.2,50f6.1)') 'ACD_PX',PF,sqrt(sum(pg**2)),sum(px(1:kdim-1))*100.,PX*100.
  93 !write(*,'(a,2f10.5,f10.2,50f7.2)') '888888 ',PF,sqrt(sum(pg**2)),sum(zx(1:kdim-1))*100.,ZX*100.
  94
  95 end subroutine da_cloud_sim