Ignore:
Timestamp:
Mar 19, 2009, 11:38:04 AM (16 years ago)
Author:
idelkadi
Message:

Corrections sur les wakes et la convection pour surmonter le probleme de l'eau negative

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/wake.F

    r1059 r1127  
    1212     o                ,Cstar,d_deltat_gw
    1313     o                ,d_deltatw2,d_deltaqw2)
     14
    1415
    1516***************************************************************
    … …  
    157158      REAL alpk
    158159      REAL delta_t_min
    159       REAL Pupper
    160160      INTEGER nsub
    161161      REAL dtimesub
    162162      REAL sigmad, hwmin
     163      REAL :: sigmaw_max
    163164cIM 080208
    164165      LOGICAL, dimension(klon) :: gwake
    … …  
    183184      INTEGER, DIMENSION(klon) :: ktop, kupper
    184185
     186c Sub-timestep tendencies and related variables
     187       REAL d_deltatw(klon,klev),d_deltaqw(klon,klev)
     188       REAL d_te(klon,klev),d_qe(klon,klev)
     189       REAL d_sigmaw(klon),alpha(klon)
     190       REAL q0_min(klon),q1_min(klon)
     191       LOGICAL wk_adv(klon), OK_qx_qw(klon)
     192       REAL epsilon
     193       DATA epsilon/1.e-9/
     194
    185195c Autres variables internes
    186196      INTEGER isubstep, k, i
    … …  
    202212      REAL, DIMENSION(klon,klev) :: the, thu
    203213
    204       REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
     214!      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
    205215
    206216      REAL, DIMENSION(klon,klev+1) :: omgbw
     217      REAL, DIMENSION(klon) :: pupper
    207218      REAL, DIMENSION(klon) :: omgtop
    208219      REAL, DIMENSION(klon,klev) :: dp_omgbw
    … …  
    279290        dqls(i,k) = 0.
    280291        d_deltat_gw(i,k)=0.
     292        d_te(i,k) = 0.
     293        d_qe(i,k) = 0.
     294        d_deltatw(i,k) = 0.
     295        d_deltaqw(i,k) = 0.
    281296!IM 060508 beg
    282297        d_deltatw2(i,k)=0.
    … …  
    294309      sigmaw(i) = amin1(sigmaw(i),0.99)
    295310      sigmaw0(i) = sigmaw(i)
     311      wape(i) = 0.
     312      wape2(i) = 0.
     313      d_sigmaw(i) = 0.
     314      ktopw(i) = 0
    296315      ENDDO
    297316C
    … …  
    406425c
    407426C       Pupper = 50000.  ! melting level
    408        Pupper = 60000.
     427c       Pupper = 60000.
    409428c       Pupper = 80000.  ! essais pour case_e
     429       DO i = 1,klon
     430ccc       Pupper(i) = 0.6*ph(i,1)
     431        Pupper(i) = 60000.
     432       ENDDO
     433
    410434C
    411435C    Determine Wake top pressure (Ptop) from buoyancy integral
    … …  
    481505      DO i=1,klon
    482506        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
    483         IF (ph(i,k+1) .lt. pupper) kupper(i)=k
     507        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
    484508      ENDDO
    485509      ENDDO
    … …  
    622646      ENDIF
    623647      ENDDO
    624 c
    625 C
     648
     649c
     650c Check qx and qw positivity
     651c --------------------------
     652      DO i = 1,klon
     653       q0_min(i)=min(  (qe(i,1)-sigmaw(i)*deltaqw(i,1)),
     654     $              (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1))  )
     655      ENDDO
     656      DO k = 2,klev
     657      DO i = 1,klon
     658        q1_min(i)=min(  (qe(i,k)-sigmaw(i)*deltaqw(i,k)),
     659     $              (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k))  )
     660        IF (q1_min(i).le.q0_min(i)) THEN
     661          q0_min(i)=q1_min(i)
     662        ENDIF
     663      ENDDO
     664      ENDDO
     665c
     666      DO i = 1,klon
     667       OK_qx_qw(i) = q0_min(i) .GE. 0.
     668       alpha(i) = 1.
     669      ENDDO
     670c
    626671CC -----------------------------------------------------------------
    627672C    Sub-time-stepping
    … …  
    634679      DO isubstep = 1,nsub
    635680c------------------------------------------------------------
    636       DO i=1,klon
     681c
     682c wk_adv is the logical flag enabling wake evolution in the time advance loop
     683      DO i = 1,klon
     684       wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.
     685      ENDDO
     686c
     687      DO i=1,klon
     688        IF (wk_adv(i)) THEN
    637689        gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))
    638       ENDDO
    639       DO i=1,klon
    640         sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
    641         sigmaw(i) =amin1(sigmaw(i),0.99)     !!!!!!!!
     690        ENDIF
     691      ENDDO
     692      DO i=1,klon
     693        IF (wk_adv(i)) THEN
     694         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     695c        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     696c        sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
    642697c        wdens = wdens0/(10.*sigmaw)
    643698c        sigmaw =max(sigmaw,sigd_con)
    644699c        sigmaw =max(sigmaw,sigmad)
     700        ENDIF
    645701      ENDDO
    646702C
    … …  
    650706cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
    651707cIM 060208 au niveau k=1..?
    652       dp_deltomg(1:klon,1:klev)=0.
     708      DO k= 1,klev
     709      DO i = 1,klon
     710        dp_deltomg(i,k)=0.
     711      ENDDO
     712      ENDDO
    653713      DO k= 1,klev+1
    654714      DO i = 1,klon
    … …  
    658718c
    659719      DO i=1,klon
     720        IF (wk_adv(i)) THEN
    660721        z(i)= 0.
    661722        omg(i,1) = 0.
    662723        dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
     724        ENDIF
    663725      ENDDO
    664726c
    665727      DO k= 2,klev
    666728      DO i = 1,klon
    667        IF( k .LE. ktop(i)) THEN
     729       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
    668730          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
    669731          z(i) = z(i)+dz(i)
    … …  
    675737c
    676738      DO i = 1,klon
     739        IF (wk_adv(i)) THEN
    677740        dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
    678741        ztop(i) = z(i)+dztop(i)
    679742        omgtop(i)=dp_deltomg(i,1)*ztop(i)
     743        ENDIF
    680744      ENDDO
    681745c
    … …  
    685749c
    686750       DO i=1,klon
     751        IF (wk_adv(i)) THEN
    687752        omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
    688753        dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
     754        ENDIF
    689755       ENDDO
    690756c
    691757      DO k= 1,klev
    692758      DO i = 1,klon
    693        IF( k .LE. ktop(i)) THEN
     759       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
    694760          omg(i,k) = - rho(i,k)*rg*omg(i,k)
    695761          dp_deltomg(i,k) = dp_deltomg(i,1)
    … …  
    701767
    702768      DO i=1,klon
    703       IF (kupper(i) .GT. ktop(i)) THEN
     769      IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
    704770        omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
    705771     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
    706772        dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
    707      $                     (ptop(i)-pupper)
     773     $                     (ptop(i)-pupper(i))
    708774      ENDIF
    709775      ENDDO
    … …  
    711777      DO k= 1,klev
    712778      DO i = 1,klon
    713        IF( k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
     779       IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
    714780          dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
    715781          omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
    … …  
    718784      ENDDO
    719785c
     786c
    720787c--    Compute wake average vertical velocity omgbw
    721788c
    … …  
    723790      DO k = 1,klev+1
    724791      DO i=1,klon
     792        IF ( wk_adv(i)) THEN
    725793        omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
     794        ENDIF
    726795      ENDDO
    727796      ENDDO
    … …  
    730799      DO k = 1,klev
    731800      DO i=1,klon
     801        IF ( wk_adv(i)) THEN
    732802        dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     803        ENDIF
    733804      ENDDO
    734805      ENDDO
    … …  
    739810      DO k = 1,klev
    740811      DO i=1,klon
    741        alpha_up(i,k) = 0.
    742        IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
     812        IF ( wk_adv(i)) THEN
     813         alpha_up(i,k) = 0.
     814         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
     815        ENDIF
    743816      ENDDO
    744817      ENDDO
    … …  
    747820
    748821      DO i=1,klon
    749         RRe1(i) = 1.-sigmaw(i)
    750         RRe2(i) = sigmaw(i)
     822        IF ( wk_adv(i)) THEN
     823         RRe1(i) = 1.-sigmaw(i)
     824         RRe2(i) = sigmaw(i)
     825        ENDIF
    751826      ENDDO
    752827      RRd1 = -1.
    … …  
    757832      DO k= 1,klev
    758833      DO i = 1,klon
    759        IF(k .LE. kupper(i)+1) THEN
     834       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
    760835        dth(i,k) = deltatw(i,k)/ppi(i,k)
    761836        Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
    … …  
    778853      DO k= 2,klev
    779854      DO i = 1,klon
    780        IF(k .LE. kupper(i)+1) THEN
     855       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
    781856        D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
    782857        D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
    … …  
    790865
    791866      DO i=1,klon
    792         omgbdth(i,1) = 0.
    793         omgbdq(i,1) = 0.
     867        IF( wk_adv(i)) THEN
     868         omgbdth(i,1) = 0.
     869         omgbdq(i,1) = 0.
     870        ENDIF
    794871      ENDDO
    795872
    796873      DO k= 2,klev
    797874      DO i = 1,klon
    798        IF(k .LE. kupper(i)+1) THEN  !   loop on interfaces
     875       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN  !   loop on interfaces
    799876        omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
    800877        omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
    … …  
    806883      DO k= 1,klev
    807884      DO i = 1,klon
    808        IF(k .LE. kupper(i)-1) THEN
     885       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    809886c-----------------------------------------------------------------
    810887c
    … …  
    829906c   and increment large scale tendencies
    830907c
    831          dtls(i,k) = dtls(i,k) +
    832      $               dtimesub*(
     908
     909c
     910C
     911CC -----------------------------------------------------------------
     912         d_te(i,k) =  dtimesub*(
    833913     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
    834914     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
    … …  
    836916     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
    837917     $                      )*ppi(i,k)
    838 c         print*,'dtls=',dtls(i,k)
    839 c
    840          dqls(i,k) = dqls(i,k) +
    841      $               dtimesub*(
     918c
     919         d_qe(i,k) =  dtimesub*(
    842920     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
    843921     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
    … …  
    845923     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
    846924     $                      )
    847 c         print*,'dqls=',dqls(k)
    848        ENDIF
     925       ENDIF
     926
    849927c-------------------------------------------------------------------
    850928      ENDDO
    … …  
    856934      DO k= 1,klev
    857935      DO i = 1,klon
    858        IF(k .LE. kupper(i)-1) THEN
     936       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    859937c
    860938c Coefficient de répartition
    … …  
    912990
    913991        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
    914           deltatw(i,k) = deltatw(i,k)+dtimesub*
    915      $          (ff(i)+dtKE(i,k)+dtPBL(i,k) 
     992          d_deltatw(i,k) = dtimesub*
     993     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
    916994     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
    917995        ELSE
    918            deltatw(i,k) = deltatw(i,k)+1/Tgw(i,k)*(1-exp(-dtimesub*
     996           d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
    919997     $          Tgw(i,k)))*
    920998     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
    921999     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
    9221000        ENDIF
    923    
     1001
    9241002        dth(i,k) = deltatw(i,k)/ppi(i,k)
    9251003
    9261004        gg(i)=d_deltaqw(i,k)/dtimesub
    9271005
    928        deltaqw(i,k) = deltaqw(i,k) +
    929      $         dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) - spread(i,k)*
    930      $         deltaqw(i,k))
     1006       d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)
     1007     $                            - spread(i,k)*deltaqw(i,k))
    9311008
    9321009       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
    … …  
    9361013      ENDDO
    9371014
    938 C   And update large scale variables
     1015C
     1016C   Scale tendencies so that water vapour remains positive in w and x.
     1017C
     1018      call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw,
     1019     $                d_deltaqw,sigmaw,d_sigmaw,alpha)
     1020c
     1021      DO k = 1,klev
     1022      DO i = 1,klon
     1023       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1024        d_te(i,k)=alpha(i)*d_te(i,k)
     1025        d_qe(i,k)=alpha(i)*d_qe(i,k)
     1026        d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
     1027        d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
     1028        d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
     1029       ENDIF
     1030      ENDDO
     1031      ENDDO
     1032      DO i = 1,klon
     1033       IF( wk_adv(i)) THEN
     1034        d_sigmaw(i)=alpha(i)*d_sigmaw(i)
     1035       ENDIF
     1036      ENDDO
     1037
     1038C   Update large scale variables and wake variables
    9391039cIM 060208 manque DO i + remplace DO k=1,kupper(i)
    9401040cIM 060208     DO k = 1,kupper(i)
    9411041      DO k= 1,klev
    9421042      DO i = 1,klon
    943        IF(k .LE. kupper(i)) THEN
     1043       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1044        dtls(i,k)=dtls(i,k)+d_te(i,k)
     1045        dqls(i,k)=dqls(i,k)+d_qe(i,k)
     1046       ENDIF
     1047      ENDDO
     1048      ENDDO
     1049      DO k= 1,klev
     1050      DO i = 1,klon
     1051       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
    9441052        te(i,k) = te0(i,k) + dtls(i,k)
    9451053        qe(i,k) = qe0(i,k) + dqls(i,k)
    9461054        the(i,k) = te(i,k)/ppi(i,k)
    947        ENDIF
    948       ENDDO
     1055        deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
     1056        deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
     1057        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1058       ENDIF
     1059      ENDDO
     1060      ENDDO
     1061      DO i = 1,klon
     1062       IF( wk_adv(i)) THEN
     1063        sigmaw(i) = sigmaw(i)+d_sigmaw(i)
     1064       ENDIF
    9491065      ENDDO
    9501066c
    … …  
    9561072c
    9571073      DO i=1,klon
    958       Ptop_provis(i)=ph(i,1)
     1074       IF ( wk_adv(i)) THEN
     1075        Ptop_provis(i)=ph(i,1)
     1076       ENDIF
    9591077      ENDDO
    9601078c
    9611079      DO k= 2,klev
    9621080      DO i=1,klon
    963         IF (Ptop_provis(i) .EQ. ph(i,1) .AND.
     1081        IF ( wk_adv(i) .AND.
     1082     $       Ptop_provis(i) .EQ. ph(i,1) .AND.
    9641083     $      dth(i,k) .GT. -delta_t_min .and.
    9651084     $      dth(i,k-1).LT. -delta_t_min) THEN
    … …  
    9811100      DO k = 1,klev
    9821101      DO i=1,klon
     1102       IF ( wk_adv(i)) THEN
    9831103        dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
    9841104        IF (dz(i) .gt. 0) THEN
    … …  
    9871107         dthmin(i) = amin1(dthmin(i),dth(i,k))
    9881108        ENDIF
     1109       ENDIF
    9891110      ENDDO
    9901111      ENDDO
    … …  
    9931114
    9941115      DO i=1,klon
    995        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
    996        hw(i) = amax1(hwmin,hw(i))
     1116       IF ( wk_adv(i)) THEN
     1117         hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
     1118         hw(i) = amax1(hwmin,hw(i))
     1119       ENDIF
    9971120      ENDDO
    9981121c
    … …  
    10061129      DO k = 1,klev
    10071130      DO i=1,klon
     1131       IF ( wk_adv(i)) THEN
    10081132        dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
    10091133        IF (dz(i) .gt. 0) THEN
    … …  
    10121136         ktop(i) = k
    10131137        ENDIF
     1138       ENDIF
    10141139      ENDDO
    10151140      ENDDO
    … …  
    10181143c
    10191144      DO i=1,klon
     1145       IF ( wk_adv(i)) THEN
    10201146        Ptop_new(i)=ptop(i)
     1147       ENDIF
    10211148      ENDDO
    10221149c
    … …  
    10241151      DO i=1,klon
    10251152cIM v3JYG; IF (k .GE. ktop(i)
    1026        IF (k .LE. ktop(i) .AND.
     1153       IF ( wk_adv(i) .AND.
     1154     $      k .LE. ktop(i) .AND.
    10271155     $      ptop_new(i) .EQ. ptop(i) .AND.
    10281156     $      dth(i,k) .GT. -delta_t_min .and.
    … …  
    10371165c
    10381166      DO i=1,klon
    1039       ptop(i) = ptop_new(i)
     1167       IF ( wk_adv(i)) THEN
     1168        ptop(i) = ptop_new(i)
     1169       ENDIF
    10401170      ENDDO
    10411171
    … …  
    10501180      DO k = 1,klev
    10511181      DO i=1,klon
    1052         IF (k .GE. kupper(i)) THEN
     1182        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
    10531183         deltatw(i,k) = 0.
    10541184         deltaqw(i,k) = 0.
    … …  
    10581188c
    10591189C
     1190c-------------Cstar computation---------------------------------
     1191      DO i=1, klon
     1192      sum_thu(i) = 0.
     1193      sum_tu(i) = 0.
     1194      sum_qu(i) = 0.
     1195      sum_thvu(i) = 0.
     1196      sum_dth(i) = 0.
     1197      sum_dq(i) = 0.
     1198      sum_rho(i) = 0.
     1199      sum_dtdwn(i) = 0.
     1200      sum_dqdwn(i) = 0.
     1201
     1202      av_thu(i) = 0.
     1203      av_tu(i) =0.
     1204      av_qu(i) =0.
     1205      av_thvu(i) = 0.
     1206      av_dth(i) = 0.
     1207      av_dq(i) = 0.
     1208      av_rho(i) =0.
     1209      av_dtdwn(i) =0.
     1210      av_dqdwn(i) = 0.
     1211      ENDDO
     1212C
     1213C Integrals (and wake top level number)
     1214C --------------------------------------
     1215C
     1216C Initialize sum_thvu to 1st level virt. pot. temp.
     1217
     1218      DO i=1,klon
     1219      z(i) = 1.
     1220      dz(i) = 1.
     1221      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     1222      sum_dth(i) = 0.
     1223      ENDDO
     1224
     1225      DO k = 1,klev
     1226      DO i=1,klon
     1227        dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     1228        IF (dz(i) .GT. 0) THEN
     1229         z(i) = z(i)+dz(i)
     1230         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     1231         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     1232         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     1233         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     1234         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     1235         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     1236         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     1237         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     1238         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     1239        ENDIF
     1240      ENDDO
     1241      ENDDO
     1242c
     1243      DO i=1,klon
     1244        hw0(i) = z(i)
     1245      ENDDO
     1246c
     1247C
     1248C - WAPE and mean forcing computation
     1249C ---------------------------------------
     1250C
     1251C ---------------------------------------
     1252C
     1253C Means
     1254
     1255      DO i=1,klon
     1256       av_thu(i) = sum_thu(i)/hw0(i)
     1257       av_tu(i) = sum_tu(i)/hw0(i)
     1258       av_qu(i) = sum_qu(i)/hw0(i)
     1259       av_thvu(i) = sum_thvu(i)/hw0(i)
     1260       av_dth(i) = sum_dth(i)/hw0(i)
     1261       av_dq(i) = sum_dq(i)/hw0(i)
     1262       av_rho(i) = sum_rho(i)/hw0(i)
     1263       av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     1264       av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     1265c
     1266       wape(i) = - rg*hw0(i)*(av_dth(i)
     1267     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
     1268     $     av_dq(i) ))/av_thvu(i)
     1269      ENDDO
     1270C
     1271C Filter out bad wakes
     1272
     1273      DO k = 1,klev
     1274       DO i=1,klon
     1275        IF ( wape(i) .LT. 0.) THEN
     1276          deltatw(i,k) = 0.
     1277          deltaqw(i,k) = 0.
     1278          dth(i,k) = 0.
     1279        ENDIF
     1280       ENDDO
     1281      ENDDO
     1282c
     1283      DO i=1,klon
     1284      IF ( wape(i) .LT. 0.) THEN
     1285        wape(i) = 0.
     1286        Cstar(i) = 0.
     1287        hw(i) = hwmin
     1288        sigmaw(i) = max(sigmad,sigd_con(i))
     1289        fip(i) = 0.
     1290        gwake(i) = .FALSE.
     1291      ELSE
     1292        Cstar(i) = stark*sqrt(2.*wape(i))
     1293        gwake(i) = .TRUE.
     1294      ENDIF
     1295      ENDDO
     1296
    10601297       ENDDO      ! end sub-timestep loop
    10611298C
    … …  
    10651302      DO k = 1,klev
    10661303      DO i=1,klon
    1067         IF (k .LE. kupper(i)-1) THEN
     1304        IF ( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    10681305         dtls(i,k) = dtls(i,k)/dtime
    10691306         dqls(i,k) = dqls(i,k)/dtime
    … …  
    11111348      DO k =1,klev
    11121349      DO i=1,klon
     1350       IF ( wk_adv(i)) THEN
    11131351        rho(i,k) = p(i,k)/(rd*te(i,k))
    11141352        IF(k .eq. 1) THEN
    … …  
    11251363        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
    11261364        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1365       ENDIF
    11271366      ENDDO
    11281367      ENDDO
    … …  
    11341373
    11351374      DO i=1,klon
     1375       IF ( wk_adv(i)) THEN
    11361376        z(i) = 1.
    11371377        dz(i) = 1.
    11381378        sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
    11391379        sum_dth(i) = 0.
     1380      ENDIF
    11401381      ENDDO
    11411382
    11421383      DO k = 1,klev
    11431384      DO i=1,klon
     1385       IF ( wk_adv(i)) THEN
    11441386        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
    11451387        IF (dz(i) .GT. 0) THEN
    … …  
    11551397         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
    11561398        ENDIF
    1157       ENDDO
    1158       ENDDO
    1159 c
    1160       DO i=1,klon
     1399       ENDIF
     1400      ENDDO
     1401      ENDDO
     1402c
     1403      DO i=1,klon
     1404       IF ( wk_adv(i)) THEN
    11611405        hw0(i) = z(i)
    1162       ENDDO
    1163 c
    1164 C 2.1 - WAPE and mean forcing computation
     1406       ENDIF
     1407      ENDDO
     1408c
     1409C - WAPE and mean forcing computation
    11651410C-------------------------------------------------------------
    11661411
    … …  
    11681413
    11691414      DO i=1, klon
     1415       IF ( wk_adv(i)) THEN
    11701416        av_thu(i) = sum_thu(i)/hw0(i)
    11711417        av_tu(i) = sum_tu(i)/hw0(i)
    … …  
    11811427     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+
    11821428     $     av_dth(i)*av_dq(i) ))/av_thvu(i)
    1183       ENDDO
    1184 
    1185 C 2.2 Prognostic variable update
     1429       ENDIF
     1430      ENDDO
     1431
     1432C Prognostic variable update
    11861433C ------------------------------------------------------------
    11871434
    … …  
    11901437      DO k = 1,klev
    11911438      DO i=1,klon
    1192         IF ( wape2(i) .LT. 0.) THEN
     1439        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
    11931440          deltatw(i,k) = 0.
    11941441          deltaqw(i,k) = 0.
    … …  
    12001447
    12011448      DO i=1, klon
    1202       IF ( wape2(i) .LT. 0.) THEN
     1449       IF ( wk_adv(i)) THEN
     1450       IF ( wape2(i) .LT. 0.) THEN
    12031451        wape2(i) = 0.
    12041452        Cstar2(i) = 0.
    … …  
    12121460        gwake(i) = .TRUE.
    12131461      ENDIF
     1462      ENDIF
    12141463      ENDDO
    12151464c
    12161465      DO i=1, klon
     1466       IF ( wk_adv(i)) THEN
    12171467        ktopw(i) = ktop(i)
     1468       ENDIF
    12181469      ENDDO
    12191470c
    12201471      DO i=1, klon
    1221       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
     1472       IF ( wk_adv(i)) THEN
     1473       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
    12221474
    12231475Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
    … …  
    12341486         FIP(i) = 0.
    12351487       ENDIF
     1488       ENDIF
    12361489      ENDDO
    12371490c
    … …  
    12411494C              alors il disparait en se mélangeant à la partie undisturbed
    12421495c
     1496      sigmaw_max = 0.9
    12431497      DO k = 1,klev
    12441498       DO i=1, klon
    1245          IF ((sigmaw(i).GT.0.9).or.
     1499c correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
     1500!         print*,'wape wape2 ktopw OK_qx_qw =',
     1501!     $           wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
     1502         IF ((sigmaw(i).GT.sigmaw_max).or.
    12461503     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
    1247      $      (ktopw(i).le.2)) THEN
     1504     $      (ktopw(i).le.2) .OR.
     1505     $     .not. OK_qx_qw(i)) THEN
    12481506cIM cf NR/JYG 251108  $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    12491507ccc      IF (sigmaw(i).GT.0.9) THEN
    … …  
    12571515c
    12581516      DO i=1, klon
    1259          IF ((sigmaw(i).GT.0.9).or.
    1260      $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
     1517         IF ( (sigmaw(i).GT.sigmaw_max).or.
     1518     $      ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
     1519     $      (ktopw(i).le.2) .OR.
     1520     $     .not. OK_qx_qw(i)) THEN
     1521! correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    12611522ccc      IF (sigmaw(i).GT.0.9) THEN
    12621523         wape(i) = 0.
    … …  
    12721533      RETURN
    12731534      END
     1535
     1536      SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe,
     1537     $           deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
     1538c------------------------------------------------------
     1539cDtermination du coefficient alpha tel que les tendances
     1540c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
     1541c a une humidite positive dans la zone (x) et dans la zone (w).
     1542c------------------------------------------------------
     1543c
     1544 
     1545c  Input
     1546      REAL qe(nlon,nl),d_qe(nlon,nl)
     1547      REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
     1548      REAL sigmaw(nlon),d_sigmaw(nlon)
     1549      LOGICAL wk_adv(nlon)
     1550      INTEGER nl,nlon
     1551c  Output
     1552      REAL alpha(nlon)
     1553c  Internal variables
     1554      REAL alpha1(nlon)
     1555      REAL x,a,b,c,discrim,zeta(nlon)
     1556      REAL epsilon
     1557!      DATA epsilon/1.e-9/
     1558c
     1559      DO k=1,nl
     1560      DO i = 1,nlon
     1561       IF (wk_adv(i)) THEN
     1562        IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
     1563         zeta(i)=0.
     1564        ELSE
     1565         zeta(i)=1.
     1566        END IF
     1567       ENDIF
     1568      ENDDO
     1569      DO i = 1,nlon
     1570       IF (wk_adv(i)) THEN
     1571        x = qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
     1572     $   +d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
     1573     $   -d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))
     1574      a=-d_sigmaw(i)*d_deltaqw(i,k)
     1575      b=d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
     1576     $           -deltaqw(i,k)*d_sigmaw(i)
     1577!      c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)-epsilon
     1578       c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
     1579
     1580      discrim=b*b-4.*a*c
     1581!       print*,'ZETA *********************' 
     1582!       print*,'zeta sigmaw ',zeta(:)
     1583!       print*,'SIGMA *********************'
     1584!       print*,'sigmaw ',sigmaw(:)
     1585
     1586!       print*,' x ************************'
     1587!       print*,'x ',x
     1588!       print*,'  a+b ************************'
     1589!       print*,'a+b ',a+b
     1590
     1591!       print*,'a b c delta zeta ',a,b,c,discrim
     1592        IF (a+b .GE. 0.) THEN
     1593         alpha1(i)=1.
     1594        ELSE
     1595         IF (x .GE. 0.) THEN
     1596            alpha1(i)=1.
     1597         ELSE
     1598!              IF (a .GE. 0.) THEN
     1599              IF (a .GT. 0.) THEN
     1600!       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
     1601!       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
     1602                 alpha1(i)=0.9*min(   (2.*c)/(-b+sqrt(discrim)),
     1603     $                        (-b+sqrt(discrim))/(2.*a)   )
     1604              ELSE IF (a.eq.0.) THEN
     1605                 alpha1(i)=0.9*(-c/b)
     1606              ELSE
     1607!       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
     1608!       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
     1609                 alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
     1610     $                        (-b+sqrt(discrim))/(2.*a)   )
     1611              ENDIF
     1612         ENDIF
     1613        ENDIF
     1614       ENDIF
     1615      ENDDO
     1616      ENDDO
     1617c
     1618      DO i = 1,nlon
     1619       IF (wk_adv(i)) THEN
     1620        alpha(i) = min(alpha(i),alpha1(i))
     1621       ENDIF
     1622      ENDDO
     1623c
     1624      return
     1625      end
     1626
    12741627      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
    12751628     :                ,te0,qe0,omgb
    … …  
    23152668C              alors il disparait en se mélangeant à la partie undisturbed
    23162669
     2670! correction NICOLAS     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
    23172671      IF ((sigmaw.GT.0.9).or.
    23182672     .     ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
Note: See TracChangeset for help on using the changeset viewer.