Ignore:
Timestamp:
Nov 5, 2018, 3:24:59 PM (6 years ago)
Author:
Laurent Fairhead
Message:

Undoing merge with trunk (r3356) to properly register Yann's latest modifications

Location:
LMDZ6/branches/DYNAMICO-conv
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/DYNAMICO-conv

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/wake.F90

    r3356 r3411  
    44SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
    55                te0, qe0, omgb, &
    6                 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
    7                 sigd_con, Cin, &
    8                 deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
     6                dtdwn, dqdwn, amdwn, amup, dta, dqa, &
     7                sigd_con, &
     8                deltatw, deltaqw, sigmaw, wdens, &                          ! state variables
    99                dth, hw, wape, fip, gfl, &
    1010                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
    1111                dtke, dqke, omg, dp_deltomg, spread, cstar, &
    1212                d_deltat_gw, &
    13                 d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)     ! tendencies
     13                d_deltatw2, d_deltaqw2, d_sigmaw2, d_wdens2)                ! tendencies
    1414
    1515
    … …  
    4848  ! dtls : large scale temperature tendency due to wake
    4949  ! dqls : large scale humidity tendency due to wake
    50   ! hw   : wake top hight (given by hw*deltatw(1)/2=wape)
     50  ! hw   : hauteur de la poche
    5151  ! dp_omgb : vertical gradient of large scale omega
    52   ! awdens  : densite de poches actives
    5352  ! wdens   : densite de poches
    5453  ! omgbdth: flux of Delta_Theta transported by LS omega
    … …  
    7372  ! dta  : source de chaleur due courants satures et detrain  (K/s)
    7473  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
    75   ! wgen : number of wakes generated per unit area and per sec (/m^2/s)
    7674  ! amdwn: flux de masse total des descentes, par unite de
    77   !        surface de la maille (kg/m2/s)
     75  ! surface de la maille (kg/m2/s)
    7876  ! amup : flux de masse total des ascendances, par unite de
    79   !        surface de la maille (kg/m2/s)
    80   ! sigd_con:
    81   ! Cin  : convective inhibition
     77  ! surface de la maille (kg/m2/s)
    8278  ! p    : pressions aux milieux des couches (Pa)
    8379  ! ph   : pressions aux interfaces (Pa)
    … …  
    109105  ! deltatw0   : deltatw initial
    110106  ! deltaqw0   : deltaqw initial
    111   ! hw0    : wake top hight (defined as the altitude at which deltatw=0)
     107  ! hw0    : hw initial
     108  ! sigmaw0: sigmaw initial
    112109  ! amflux : horizontal mass flux through wake boundary
    113110  ! wdens_ref: initial number of wakes per unit area (3D) or per
    … …  
    136133  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
    137134  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
    138   REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
    139135  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
    140   REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
    141136
    142137  !
    … …  
    145140  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
    146141  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
    147   REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
    148142  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    149143
    … …  
    155149  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
    156150  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
    157   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
     151  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread
    158152  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
    159153  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
    … …  
    163157  ! Tendencies of state variables
    164158  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
    165   REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
     159  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_wdens2
    166160
    167161  ! Variables internes
    … …  
    171165  INTEGER, SAVE                                         :: igout
    172166  !$OMP THREADPRIVATE(igout)
     167  REAL                                                  :: alon
    173168  LOGICAL, SAVE                                         :: first = .TRUE.
    174169  !$OMP THREADPRIVATE(first)
    … …  
    181176  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
    182177
    183   REAL, SAVE                                            :: tau_cv
    184   !$OMP THREADPRIVATE(tau_cv)
    185   REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
    186   !$OMP THREADPRIVATE(rzero, aa0)
    187 
    188178  LOGICAL, SAVE                                         :: flag_wk_check_trgl
    189179  !$OMP THREADPRIVATE(flag_wk_check_trgl)
    190180  INTEGER, SAVE                                         :: iflag_wk_check_trgl
    191181  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
    192   INTEGER, SAVE                                         :: iflag_wk_pop_dyn
    193   !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
    194182
    195183  REAL                                                  :: delta_t_min
    196184  INTEGER                                               :: nsub
    197185  REAL                                                  :: dtimesub
    198   REAL                                                  :: wdensmin
    199   REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
    200   !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
     186  REAL                                                  :: sigmad, hwmin, wapecut
    201187  REAL                                                  :: sigmaw_max
    202188  REAL                                                  :: dens_rate
    … …  
    209195  REAL, DIMENSION (klon, klev)                          :: deltaqw0
    210196  REAL, DIMENSION (klon, klev)                          :: te, qe
     197  REAL, DIMENSION (klon)                                :: sigmaw0
    211198!!  REAL, DIMENSION (klon)                                :: sigmaw1
    212 
    213   ! Variables liees a la dynamique de population
    214   REAL, DIMENSION(klon)                                 :: act
    215   REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
    216   REAL, DIMENSION(klon)                                 :: f_shear
    217   REAL, DIMENSION(klon)                                 :: drdt
    218   REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
    219   REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
    220   LOGICAL, DIMENSION (klon)                             :: kill_wake
    221   INTEGER, SAVE                                         :: iflag_wk_act
    222   !$OMP THREADPRIVATE(iflag_wk_act)
    223   REAL                                                  :: drdt_pos
    224   REAL                                                  :: tau_wk_inv_min
    225199
    226200  ! Variables pour les GW
    … …  
    230204  REAL, DIMENSION (klon, klev)                          :: tgw
    231205
    232   ! Variables liees au calcul de hw
     206  ! Variables liées au calcul de hw
    233207  REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
    234208  REAL, DIMENSION (klon)                                :: sum_dth
    … …  
    237211  INTEGER, DIMENSION (klon)                             :: ktop, kupper
    238212
    239   ! Variables liees au test de la forme triangulaire du profil de Delta_theta
     213  ! Variables liées au test de la forme triangulaire du profil de Delta_theta
    240214  REAL, DIMENSION (klon)                                :: sum_half_dth
    241215  REAL, DIMENSION (klon)                                :: dz_half
    … …  
    244218  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
    245219  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
    246   REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
    247220  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
    248221  REAL, DIMENSION (klon)                                :: q0_min, q1_min
    … …  
    255228  INTEGER                                               ::isubstep, k, i
    256229
    257   REAL                                                  :: wdens_targ
    258230  REAL                                                  :: sigmaw_targ
    259231
    … …  
    301273  REAL, DIMENSION (klon, klev)                          :: detr
    302274
    303   REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
    304   REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
     275  REAL, DIMENSION(klon)                                 :: sigmaw_in   ! pour les prints
    305276
    306277  ! -------------------------------------------------------------------------
    … …  
    313284  ! -------------------------------------------------------------------------
    314285
    315 !!  DATA wapecut, sigmad, hwmin/5., .02, 10./
    316   DATA wapecut, sigmad, hwmin/1., .02, 10./
    317 !!  DATA wdensmin/1.e-12/
    318   DATA wdensmin/1.e-14/
     286  DATA wapecut, sigmad, hwmin/5., .02, 10./
    319287  ! cc nrlmd
    320288  DATA sigmaw_max/0.4/
    321289  DATA dens_rate/0.1/
    322   DATA rzero /5000./
    323290  ! cc
    324291  ! Longueur de maille (en m)
    … …  
    326293
    327294  ! ALON = 3.e5
    328   ! alon = 1.E6
    329 
    330   !  Provisionnal; to be suppressed when f_shear is parameterized
    331   f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
     295  alon = 1.E6
    332296
    333297
    … …  
    336300  ! coefgw : Coefficient pour les ondes de gravité
    337301  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    338   ! wdens : Densité surfacique de poche froide
     302  ! wdens : Densité de poche froide par maille
    339303  ! -------------------------------------------------------------------------
    340304
    … …  
    357321  crep_sol = 1.0
    358322
    359   aa0 = 3.14*rzero*rzero
    360 
    361   tau_cv = 4000.
    362 
    363323  ! cc nrlmd Lecture du fichier wake_param.data
    364324  stark=0.33
    365325  CALL getin_p('stark',stark)
    366   cstart = stark*sqrt(2.*wapecut)
    367 
    368326  alpk=0.25
    369327  CALL getin_p('alpk',alpk)
    … …  
    376334  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
    377335!>jyg
    378   iflag_wk_pop_dyn = 0
    379   CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed
    380                                                     ! and wdens prognostic
    381   iflag_wk_act = 0
    382   CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
    383                                             ! 1: act(:)=1.
    384                                             ! 2: act(:)=f(Wape)
    385336  coefgw=4.
    386337  CALL getin_p('coefgw',coefgw)
    … …  
    393344  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
    394345!>jyg
    395   WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
    396   WRITE(*,*) 'iflag_wk_act',iflag_wk_act
    397346  WRITE(*,*) 'coefgw=', coefgw
    398347
    … …  
    408357 endif
    409358
    410  IF (iflag_wk_pop_dyn == 0) THEN
    411359  ! Initialisation de toutes des densites a wdens_ref.
    412360  ! Les densites peuvent evoluer si les poches debordent
    413361  ! (voir au tout debut de la boucle sur les substeps)
    414   !jyg<
    415   !!  wdens(:) = wdens_ref
    416    DO i = 1,klon
    417      wdens(i) = wdens_ref(znatsurf(i)+1)
    418    ENDDO
    419   !>jyg
    420  ENDIF  ! (iflag_wk_pop_dyn == 0)
     362!jyg<
     363!!  wdens(:) = wdens_ref
     364  DO i = 1,klon
     365    wdens(i) = wdens_ref(znatsurf(i)+1)
     366  ENDDO
     367!>jyg
    421368
    422369  ! print*,'stark',stark
    … …  
    468415      d_deltatw2(:,:) = 0.
    469416      d_deltaqw2(:,:) = 0.
    470 
    471       IF (iflag_wk_act == 0) THEN
    472         act(:) = 0.
    473       ELSEIF (iflag_wk_act == 1) THEN
    474         act(:) = 1.
    475       ENDIF
    476 
    477417!!  DO i = 1, klon
    478418!!   sigmaw_in(i) = sigmaw(i)
    … …  
    485425  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
    486426  ! ENDIF
    487   IF (iflag_wk_pop_dyn >=1) THEN
    488     DO i = 1, klon
    489       wdens_targ = max(wdens(i),wdensmin)
    490       d_wdens2(i) = wdens_targ - wdens(i)
    491       wdens(i) = wdens_targ
    492     END DO
    493   ELSE
    494     DO i = 1, klon
    495       d_awdens2(i) = 0.
    496       d_wdens2(i) = 0.
    497     END DO
    498   ENDIF  ! (iflag_wk_pop_dyn >=1)
    499 !
    500427  DO i = 1, klon
    501428    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    … …  
    507434    sigmaw(i) = sigmaw_targ
    508435!>jyg
    509   END DO
    510 
    511 !
    512   IF (iflag_wk_pop_dyn >= 1) THEN
    513     awdens_in(:) = awdens(:)
    514     wdens_in(:) = wdens(:)
    515 !!    wdens(:) = wdens(:) + wgen(:)*dtime
    516 !!    d_wdens2(:) = wgen(:)*dtime
    517 !!  ELSE
    518   ENDIF  ! (iflag_wk_pop_dyn >= 1)
    519 
    520   wape(:) = 0.
    521   wape2(:) = 0.
    522   d_sigmaw(:) = 0.
    523   ktopw(:) = 0
     436    sigmaw0(i) = sigmaw(i)
     437    wape(i) = 0.
     438    wape2(i) = 0.
     439    d_sigmaw(i) = 0.
     440    d_wdens2(i) = 0.
     441    ktopw(i) = 0
     442  END DO
    524443!
    525444!<jyg
    … …  
    914833      gwake(i) = .FALSE.
    915834    ELSE
    916       hw(i) = hw0(i)
    917835      cstar(i) = stark*sqrt(2.*wape(i))
    918836      gwake(i) = .TRUE.
    … …  
    973891    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
    974892    ! aurait un entrainement nul ---
    975     !jyg<
    976     ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
    977     ! les poches sont insuffisantes pour accueillir tout le flux de masse
    978     ! des descentes unsaturees. Nous faisons alors l'hypothese que la
    979     ! convection profonde cree directement de nouvelles poches, sans passer
    980     ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
    981 
    982893    DO i = 1, klon
    983894      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
    … …  
    988899        wdens0 = (sigmaw(i)/(4.*3.14))* &
    989900          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
    990         IF (prt_level >= 10) THEN
    991              print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
    992                      omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
    993         ENDIF
    994901        IF (wdens(i)<=wdens0*1.1) THEN
    995           IF (iflag_wk_pop_dyn >= 1) THEN
    996              d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
    997           ENDIF
    998902          wdens(i) = wdens0
    999903        END IF
    1000       END IF
    1001     END DO
     904        ! c        print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
     905        ! c     $     ,ph(i,1)-pupper(i)',
     906        ! c     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
     907        ! c     $     ,ph(i,1)-pupper(i)
     908      END IF
     909    END DO
     910
     911    ! cc nrlmd
    1002912
    1003913    DO i = 1, klon
    1004914      IF (wk_adv(i)) THEN
    1005915        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
    1006         rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
    1007916!jyg<
    1008917!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
    … …  
    1014923    END DO
    1015924
    1016     IF (iflag_wk_pop_dyn >= 1) THEN
    1017 
    1018       IF (iflag_wk_act ==2) THEN
    1019       DO i = 1, klon
    1020         IF (wk_adv(i)) THEN
    1021           wape1_act(i) = abs(cin(i))
    1022           wape2_act(i) = 2.*wape1_act(i) + 1.
    1023           act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
    1024         ENDIF  ! (wk_adv(i))
    1025       ENDDO
    1026       ENDIF  ! (iflag_wk_act ==2)
    1027 
    1028 
    1029       DO i = 1, klon
    1030         IF (wk_adv(i)) THEN
    1031 !!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
    1032           tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
    1033           tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
    1034           drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
    1035                     (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
    1036 !!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
    1037           drdt_pos=max(drdt(i),0.)
    1038 
    1039 !!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
    1040 !!                     - wdens(i)*tau_wk_inv_min &
    1041 !!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
    1042           d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
    1043           d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
    1044                          2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
    1045           d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
    1046 
    1047 !!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
    1048 !!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
    1049 !!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
    1050           d_sig_gen(i) = wgen(i)*aa0
    1051           d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
    1052 !!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
    1053           d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
    1054           d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
    1055           d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
    1056         ENDIF
    1057       ENDDO
    1058 
    1059       IF (prt_level >= 10) THEN
    1060         print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
    1061                        cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
    1062         print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
    1063                        wdens(1), awdens(1), act(1), d_awdens(1)
    1064         print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
    1065                        wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
    1066         print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
    1067                        d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
    1068       ENDIF
    1069    
    1070     ELSE  !  (iflag_wk_pop_dyn >= 1)
    1071 
    1072     ! cc nrlmd
    1073 
    1074       DO i = 1, klon
    1075         IF (wk_adv(i)) THEN
    1076           ! cc nrlmd          Introduction du taux de mortalité des poches et
    1077           ! test sur sigmaw_max=0.4
    1078           ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
    1079           IF (sigmaw(i)>=sigmaw_max) THEN
    1080             death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
    1081           ELSE
    1082             death_rate(i) = 0.
    1083           END IF
    1084    
    1085           d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
    1086             dtimesub
    1087           ! $              - nat_rate(i)*sigmaw(i)*dtimesub
    1088           ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
    1089           ! c     $  death_rate(i),ktop(i),kupper(i)',
    1090           ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
    1091           ! c     $  death_rate(i),ktop(i),kupper(i)
    1092    
    1093           ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
    1094           ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
    1095           ! wdens = wdens0/(10.*sigmaw)
    1096           ! sigmaw =max(sigmaw,sigd_con)
    1097           ! sigmaw =max(sigmaw,sigmad)
    1098         END IF
    1099       END DO
    1100 
    1101     ENDIF   !  (iflag_wk_pop_dyn >= 1)
    1102 
     925    DO i = 1, klon
     926      IF (wk_adv(i)) THEN
     927        ! cc nrlmd          Introduction du taux de mortalité des poches et
     928        ! test sur sigmaw_max=0.4
     929        ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     930        IF (sigmaw(i)>=sigmaw_max) THEN
     931          death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
     932        ELSE
     933          death_rate(i) = 0.
     934        END IF
     935
     936        d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
     937          dtimesub
     938        ! $              - nat_rate(i)*sigmaw(i)*dtimesub
     939        ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     940        ! c     $  death_rate(i),ktop(i),kupper(i)',
     941        ! c     $                d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     942        ! c     $  death_rate(i),ktop(i),kupper(i)
     943
     944        ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     945        ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     946        ! wdens = wdens0/(10.*sigmaw)
     947        ! sigmaw =max(sigmaw,sigd_con)
     948        ! sigmaw =max(sigmaw,sigmad)
     949      END IF
     950    END DO
    1103951
    1104952    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    … …  
    13751223
    13761224    ! Increment state variables
    1377 !jyg<
    1378     IF (iflag_wk_pop_dyn >= 1) THEN
    1379       DO k = 1, klev
    1380         DO i = 1, klon
    1381           IF (wk_adv(i) .AND. k<=kupper(i)) THEN
    1382             detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
    1383             entr(i,k) = d_sig_gen(i)
    1384           ENDIF
    1385         ENDDO
    1386       ENDDO
    1387       ELSE  ! (iflag_wk_pop_dyn >= 1)
    1388       DO k = 1, klev
    1389         DO i = 1, klon
    1390           IF (wk_adv(i) .AND. k<=kupper(i)) THEN
    1391             detr(i, k) = 0.
    1392    
    1393             entr(i, k) = 0.
    1394           ENDIF
    1395         ENDDO
    1396       ENDDO
    1397     ENDIF  ! (iflag_wk_pop_dyn >= 1)
    1398 
    1399    
    14001225
    14011226    DO k = 1, klev
    … …  
    14391264          ! cc nrlmd          Prise en compte du taux de mortalité
    14401265          ! cc               Définitions de entr, detr
    1441 !jyg<
    1442 !!            detr(i, k) = 0.
    1443 !!   
    1444 !!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
    1445 !!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
    1446 !!
    1447             entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
    1448                          sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
    1449 !>jyg
    1450             spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
    1451 
     1266          detr(i, k) = 0.
     1267
     1268          entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
     1269            sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
     1270
     1271          spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
    14521272          ! cc        spread(i,k) =
    14531273          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
    … …  
    15641384      END DO
    15651385    END DO
    1566 !
    15671386    DO i = 1, klon
    15681387      IF (wk_adv(i)) THEN
    15691388        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     1389!jyg<
    15701390        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
    1571       END IF
    1572     END DO
    1573 !jyg<
    1574     IF (iflag_wk_pop_dyn >= 1) THEN
    1575       DO i = 1, klon
    1576         IF (wk_adv(i)) THEN
    1577           awdens(i) = awdens(i) + d_awdens(i)
    1578           wdens(i) = wdens(i) + d_wdens(i)
    1579           d_awdens2(i) = d_awdens2(i) + d_awdens(i)
    1580           d_wdens2(i) = d_wdens2(i) + d_wdens(i)
    1581         END IF
    1582       END DO
    1583       DO i = 1, klon
    1584         IF (wk_adv(i)) THEN
    1585           wdens_targ = max(wdens(i),wdensmin)
    1586           d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
    1587           wdens(i) = wdens_targ
    1588 !
    1589           wdens_targ = min( max(awdens(i),0.), wdens(i) )
    1590           d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
    1591           awdens(i) = wdens_targ
    1592         END IF
    1593       END DO
    1594       DO i = 1, klon
    1595         IF (wk_adv(i)) THEN
    1596           sigmaw_targ = max(sigmaw(i),sigmad)
    1597           d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    1598           sigmaw(i) = sigmaw_targ
    1599         END IF
    1600       END DO
    1601     ENDIF  ! (iflag_wk_pop_dyn >= 1)
    16021391!>jyg
     1392      END IF
     1393    END DO
    16031394
    16041395
    … …  
    21101901  ! ENDDO
    21111902  ! cc
    2112 
    2113   !jyg<
    2114   IF (iflag_wk_pop_dyn >= 1) THEN
    2115     DO i = 1, klon
    2116       kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
    2117           .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
    2118     ENDDO
    2119   ELSE  ! (iflag_wk_pop_dyn >= 1)
    2120     DO i = 1, klon
    2121       kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
    2122           .NOT. ok_qx_qw(i)
    2123     ENDDO
    2124   ENDIF  ! (iflag_wk_pop_dyn >= 1)
    2125   !>jyg
    2126 
    21271903  DO k = 1, klev
    21281904    DO i = 1, klon
    2129 !!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
    2130 !!jyg          .NOT. ok_qx_qw(i)) THEN
    2131       IF (kill_wake(i)) THEN
     1905
     1906      ! cc nrlmd      On maintient désormais constant sigmaw en régime
     1907      ! permanent
     1908      ! cc      IF ((sigmaw(i).GT.sigmaw_max).or.
     1909      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
     1910          .NOT. ok_qx_qw(i)) THEN
    21321911        ! cc
    21331912        dtls(i, k) = 0.
    … …  
    21371916        d_deltatw2(i,k) = -deltatw0(i,k)
    21381917        d_deltaqw2(i,k) = -deltaqw0(i,k)
    2139       END IF  ! (kill_wake(i))
    2140     END DO
    2141   END DO
    2142 
    2143   DO i = 1, klon
    2144 !!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
    2145 !!jyg        .NOT. ok_qx_qw(i)) THEN
    2146       IF (kill_wake(i)) THEN
     1918      END IF
     1919    END DO
     1920  END DO
     1921
     1922  ! cc nrlmd      On maintient désormais constant sigmaw en régime permanent
     1923  DO i = 1, klon
     1924    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
     1925        .NOT. ok_qx_qw(i)) THEN
    21471926      ktopw(i) = 0
    21481927      wape(i) = 0.
    21491928      cstar(i) = 0.
    2150 !!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
     1929!!jyg   Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
    21511930!!      hw(i) = hwmin                       !jyg
    21521931!!      sigmaw(i) = sigmad                  !jyg
    21531932      hw(i) = 0.                            !jyg
     1933      sigmaw(i) = 0.                        !jyg
    21541934      fip(i) = 0.
    2155 !!      sigmaw(i) = 0.                        !jyg
    2156       sigmaw_targ = 0.
    2157       d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    2158       sigmaw(i) = sigmaw_targ
    2159       IF (iflag_wk_pop_dyn >= 1) THEN
    2160 !!        awdens(i) = 0.
    2161 !!        wdens(i) = 0.
    2162         wdens_targ = 0.
    2163         d_wdens2(i) = wdens_targ - wdens(i)
    2164         wdens(i) = wdens_targ
    2165         wdens_targ = 0.
    2166         d_awdens2(i) = wdens_targ - awdens(i)
    2167         awdens(i) = wdens_targ
    2168       ENDIF  ! (iflag_wk_pop_dyn >= 1)
    2169     ELSE  ! (kill_wake(i))
     1935    ELSE
    21701936      wape(i) = wape2(i)
    21711937      cstar(i) = cstar2(i)
    2172     END IF  ! (kill_wake(i))
     1938    END IF
    21731939    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
    21741940    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
    … …  
    22061972  DO i = 1, klon
    22071973    d_sigmaw2(i) = d_sigmaw2(i)/dtime
    2208     d_awdens2(i) = d_awdens2(i)/dtime
    22091974    d_wdens2(i) = d_wdens2(i)/dtime
    22101975  ENDDO
Note: See TracChangeset for help on using the changeset viewer.