Changeset 1326


Ignore:
Timestamp:
Mar 15, 2010, 5:37:02 PM (15 years ago)
Author:
idelkadi
Message:

Optimisation de thermcell_dv2.F90

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dv2.F90

    r1146 r1326  
    1010!   de "thermiques" explicitement representes
    1111!   calcul du dq/dt une fois qu'on connait les ascendances
     12!
     13! Vectorisation, FH : 2010/03/08
    1214!
    1315!=======================================================================
     
    3133      real qa(ngrid,nlay),detr(ngrid,nlay),zf,zf2
    3234      real wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
    33       real gamma0,gamma(ngrid,nlay+1)
     35      real gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
    3436      real ue(ngrid,nlay),ve(ngrid,nlay)
    35       real dua,dva
     37      LOGICAL ltherm(ngrid,nlay)
     38      real dua(ngrid,nlay),dva(ngrid,nlay)
    3639      integer iter
    3740
    38       integer ig,k
     41      integer ig,k,nlarga0
     42
     43!-------------------------------------------------------------------------
    3944
    4045!   calcul du detrainement
     46!---------------------------
     47
     48      print*,'THERMCELL DV2 OPTIMISE 3'
     49
     50      nlarga0=0.
    4151
    4252      do k=1,nlay
     
    5969      do k=2,nlay
    6070         do ig=1,ngrid
    61             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
    62      &         1.e-5*masse(ig,k)) then
    63 !   On itère sur la valeur du coeff de freinage.
    64 !              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
    65 !IM 060508 beg
    66 !             if(0.5*(fraca(ig,k+1)+fraca(ig,k)).LT.0.) THEN
    67 !              print*,'th_dv2 ig k fraca(:,k) fraca(:k+1)', &
    68 !    &         ig,k,fraca(ig,k),fraca(ig,k+1)
    69 !             endif
    70 !             if(larga(ig).EQ.0.) THEN
    71 !              print*,'th_dv2 ig larga=0.',ig
    72 !             endif
    73               if(larga(ig).GT.0.) THEN
    74 !IM 060508 end
    75                gamma0=masse(ig,k)  &
     71            ltherm(ig,k)=(fm(ig,k+1)+detr(ig,k))*ptimestep > 1.e-5*masse(ig,k)
     72            if(ltherm(ig,k).and.larga(ig)>0.) then
     73               gamma0(ig,k)=masse(ig,k)  &
    7674     &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
    7775     &         *0.5/larga(ig)  &
    7876     &         *1.
    79 !IM 060508 beg
    80               else
    81                if(prt_level.GE.10) print*,'WARNING cas ELSE on initialise gamma0=0.'
    82                gamma0=0.
    83               endif !(larga(ig).GT.0.) THEN
    84 !IM 060508 end
    85 !    s         *0.5
    86 !              gamma0=0.
    87                zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
     77            else
     78               gamma0(ig,k)=0.
     79            endif
     80            if (ltherm(ig,k).and.larga(ig)<=0.) nlarga0=nlarga0+1
     81         enddo
     82      enddo
     83
     84      gamma(:,:)=0.
     85
     86      do k=2,nlay
     87
     88         do ig=1,ngrid
     89            if (ltherm(ig,k)) then
     90               dua(ig,k)=ua(ig,k-1)-u(ig,k-1)
     91               dva(ig,k)=va(ig,k-1)-v(ig,k-1)
     92            else
     93               ua(ig,k)=u(ig,k)
     94               va(ig,k)=v(ig,k)
     95               ue(ig,k)=u(ig,k)
     96               ve(ig,k)=v(ig,k)
     97            endif
     98         enddo
     99
     100
     101! Debut des iterations
     102!----------------------
     103do iter=1,5
     104         do ig=1,ngrid
     105! Pour memoire : calcul prenant en compte la fraction reelle
     106!              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
     107!              zf2=1./(1.-zf)
     108! Calcul avec fraction infiniement petite
    88109               zf=0.
    89                zf2=1./(1.-zf)
    90 !   la première fois on multiplie le coefficient de freinage
    91  par le module du vent dans la couche en dessous.
    92                dua=ua(ig,k-1)-u(ig,k-1)
    93                dva=va(ig,k-1)-v(ig,k-1)
    94                do iter=1,5
     110               zf2=1.
     111
     112la première fois on multiplie le coefficient de freinage
     113!  par le module du vent dans la couche en dessous.
     114!  Mais pourquoi donc ???
     115               if (ltherm(ig,k)) then
    95116!   On choisit une relaxation lineaire.
    96                   gamma(ig,k)=gamma0
     117!                 gamma(ig,k)=gamma0(ig,k)
    97118!   On choisit une relaxation quadratique.
    98                   gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
     119                  gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
    99120                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)  &
    100121     &               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))  &
     
    105126     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
    106127     &                 +gamma(ig,k))
    107 !                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
    108                   dua=ua(ig,k)-u(ig,k)
    109                   dva=va(ig,k)-v(ig,k)
     128!                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua(ig,k),dva(ig,k)
     129                  dua(ig,k)=ua(ig,k)-u(ig,k)
     130                  dva(ig,k)=va(ig,k)-v(ig,k)
    110131                  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
    111132                  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
    112                enddo
    113             else
    114                ua(ig,k)=u(ig,k)
    115                va(ig,k)=v(ig,k)
    116                ue(ig,k)=u(ig,k)
    117                ve(ig,k)=v(ig,k)
    118                gamma(ig,k)=0.
    119             endif
     133               endif
    120134         enddo
    121       enddo
     135! Fin des iterations
     136!--------------------
     137enddo
    122138
     139      enddo ! k=2,nlay
     140
     141
     142! Calcul du flux vertical de moment dans l'environnement.
     143!---------------------------------------------------------
    123144      do k=2,nlay
    124145         do ig=1,ngrid
     
    134155      enddo
    135156
     157! calcul des tendances.
     158!-----------------------
    136159      do k=1,nlay
    137160         do ig=1,ngrid
    138 !IM
    139          if(prt_level.GE.10) print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
    140      &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
    141      &   masse(ig,k)
    142 !
    143161            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
    144162     &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
     
    152170      enddo
    153171
     172
     173! Sorties eventuelles.
     174!----------------------
     175
     176   if(prt_level.GE.10) then
     177      do k=1,nlay
     178         do ig=1,ngrid
     179           print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
     180     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
     181     &   masse(ig,k)
     182         enddo
     183      enddo
     184   endif
     185!
     186     if (nlarga0>0) then
     187          print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
     188          print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
     189          print*,'Il faudrait decortiquer ces points'
     190     endif
     191
    154192      return
    155193      end
Note: See TracChangeset for help on using the changeset viewer.