Changeset 4230


Ignore:
Timestamp:
Aug 23, 2022, 8:06:04 PM (21 months ago)
Author:
fhourdin
Message:

New code for the wakes. Ready for new adventures.
Jean-Yves, Lamine, Frederic

Location:
LMDZ6/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90

    r3577 r4230  
    162162     
    163163       IF (prt_level > 9) THEN
    164           WRITE(lunout,*) 'USTAR = ',yustar
     164          WRITE(lunout,*) 'USTAR = ',(yustar(i),i=1,knon)
    165165       ENDIF
    166166         
  • LMDZ6/trunk/libf/phylmd/wake.F90

    r4089 r4230  
    1111                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
    1212                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
    13                 d_deltat_gw, &
    14                 d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)     ! tendencies
     13                d_deltat_gw, &                                                      ! tendencies
     14                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)             ! tendencies
    1515
    1616
     
    2727  USE wake_ini_mod , ONLY : wake_ini
    2828  USE wake_ini_mod , ONLY : prt_level,epsim1,RG,RD
    29   USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, wdens_ref, stark, coefgw, alpk
     29  USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, pupperbyphs
    3030  USE wake_ini_mod , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
    3131  USE wake_ini_mod , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin
     
    119119  ! wdens_ref: initial number of wakes per unit area (3D) or per
    120120  ! unit length (2D), at the beginning of each time step
    121   ! Tgw    : 1 sur la période de onde de gravité
    122   ! Cgw    : vitesse de propagation de onde de gravité
     121  ! Tgw    : 1 sur la periode de onde de gravite
     122  ! Cgw    : vitesse de propagation de onde de gravite
    123123  ! LL     : distance entre 2 poches
    124124
    125125  ! -------------------------------------------------------------------------
    126   ! Déclaration de variables
     126  ! Declaration de variables
    127127  ! -------------------------------------------------------------------------
    128128
     
    163163  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
    164164  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
    165   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
    166165  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
    167166  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
    168   ! Tendencies of state variables
     167  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     168  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     169  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
    169170  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
    170171  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
     
    173174  ! -------------------
    174175
    175   ! Variables à fixer
     176  ! Variables a fixer
    176177
    177178  REAL                                                  :: delta_t_min
     
    193194  REAL, DIMENSION(klon)                                 :: f_shear
    194195  REAL, DIMENSION(klon)                                 :: drdt
    195   REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
     196!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
    196197  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
    197198  LOGICAL, DIMENSION (klon)                             :: kill_wake
    198199  REAL                                                  :: drdt_pos
    199200  REAL                                                  :: tau_wk_inv_min
     201  ! Some components of the tendencies of state variables 
     202  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_bnd2
     203  REAL, DIMENSION (klon)                                :: d_sig_spread2
     204  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
    200205
    201206  ! Variables pour les GW
     
    219224  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
    220225  REAL, DIMENSION (klon, klev)                          :: d_tenv, d_qe
    221   REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
    222   REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
     226  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens, d_sigmaw
     227  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     228  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     229  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
    223230  REAL, DIMENSION (klon)                                :: q0_min, q1_min
    224231  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
     
    227234  INTEGER                                               ::isubstep, k, i, igout
    228235
     236  REAL                                                  :: sigmaw_targ
    229237  REAL                                                  :: wdens_targ
    230   REAL                                                  :: sigmaw_targ
     238  REAL                                                  :: d_sigmaw_targ
     239  REAL                                                  :: d_wdens_targ
    231240
    232241  REAL, DIMENSION (klon)                                :: sum_thu, sum_tu, sum_qu, sum_thvu
     
    289298  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
    290299
    291   ! coefgw : Coefficient pour les ondes de gravité
     300  ! coefgw : Coefficient pour les ondes de gravite
    292301  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    293   ! wdens : Densité surfacique de poche froide
     302  ! wdens : Densite surfacique de poche froide
    294303  ! -------------------------------------------------------------------------
    295304
     
    380389   sigmaw_in(:) = sigmaw(:)
    381390!>jyg
    382 
    383   ! sigmaw1=sigmaw
    384   ! IF (sigd_con.GT.sigmaw1) THEN
    385   ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
    386   ! ENDIF
    387   IF (iflag_wk_pop_dyn >=1) THEN
    388     DO i = 1, klon
    389       wdens_targ = max(wdens(i),wdensmin)
    390       d_wdens2(i) = wdens_targ - wdens(i)
    391       wdens(i) = wdens_targ
    392     END DO
    393   ELSE
    394     DO i = 1, klon
    395       d_awdens2(i) = 0.
    396       d_wdens2(i) = 0.
    397     END DO
    398   ENDIF  ! (iflag_wk_pop_dyn >=1)
    399 !
    400   DO i = 1, klon
    401     ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    402 !jyg<
    403 !!    sigmaw(i) = amax1(sigmaw(i), sigmad)
    404 !!    sigmaw(i) = amin1(sigmaw(i), 0.99)
    405     sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
    406     d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
    407     sigmaw(i) = sigmaw_targ
    408 !>jyg
    409   END DO
    410 
    411391!
    412392  IF (iflag_wk_pop_dyn >= 1) THEN
     
    417397!!  ELSE
    418398  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     399
     400
     401  ! sigmaw1=sigmaw
     402  ! IF (sigd_con.GT.sigmaw1) THEN
     403  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     404  ! ENDIF
     405  IF (iflag_wk_pop_dyn >=1) THEN
     406    DO i = 1, klon
     407      d_dens_gen2(i)   = 0.
     408      d_dens_death2(i) = 0.
     409      d_dens_col2(i)   = 0.
     410      d_awdens2(i) = 0.
     411      wdens_targ = max(wdens(i),wdensmin)
     412      d_dens_bnd2(i) = wdens_targ - wdens(i)
     413      d_wdens2(i) = wdens_targ - wdens(i)
     414      wdens(i) = wdens_targ
     415    END DO
     416  ELSE
     417    DO i = 1, klon
     418      d_awdens2(i) = 0.
     419      d_wdens2(i) = 0.
     420    END DO
     421  ENDIF  ! (iflag_wk_pop_dyn >=1)
     422!
     423  DO i = 1, klon
     424    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
     425!jyg<
     426!!    sigmaw(i) = amax1(sigmaw(i), sigmad)
     427!!    sigmaw(i) = amin1(sigmaw(i), 0.99)
     428    d_sig_gen2(i)   = 0.
     429    d_sig_death2(i) = 0.
     430    d_sig_col2(i)   = 0.
     431    d_sig_spread2(i)= 0.
     432    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
     433    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
     434    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
     435    sigmaw(i) = sigmaw_targ
     436!>jyg
     437  END DO
    419438
    420439  wape(:) = 0.
     
    572591  END DO
    573592
    574 
     593 
    575594  ! Choose an integration bound well above wake top
    576595  ! -----------------------------------------------------------------
    577 
    578   ! Pupper = 50000.  ! melting level
    579   ! Pupper = 60000.
    580   ! Pupper = 80000.  ! essais pour case_e
    581   DO i = 1, klon
    582     pupper(i) = 0.6*ph(i, 1)
    583     pupper(i) = max(pupper(i), 45000.)
    584     ! cc        Pupper(i) = 60000.
    585   END DO
    586 
    587596
    588597  ! Determine Wake top pressure (Ptop) from buoyancy integral
     
    657666  ! -5/ Determination de ktop et kupper
    658667
     668  CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
     669 
    659670  DO k = klev, 1, -1
    660671    DO i = 1, klon
    661672      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
    662       IF (ph(i,k+1)<pupper(i)) kupper(i) = k
    663     END DO
    664   END DO
    665 
    666   ! On evite kupper = 1 et kupper = klev
    667   DO i = 1, klon
    668     kupper(i) = max(kupper(i), 2)
    669     kupper(i) = min(kupper(i), klev-1)
    670   END DO
     673    END DO
     674  END DO
     675  print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
     676 
    671677
    672678
     
    808814!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
    809815      sigmaw_targ = max(sigmad, sigd_con(i))
     816      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    810817      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    811818      sigmaw(i) = sigmaw_targ
     
    840847    ok_qx_qw(i) = q0_min(i) >= 0.
    841848    alpha(i) = 1.
     849    alpha_tot(i) = 1.
    842850  END DO
    843851
     
    855863  dtimesub = dtime/nsub
    856864
     865
     866 
    857867  ! ------------------------------------------------------------
    858868  DO isubstep = 1, nsub
    859869    ! ------------------------------------------------------------
     870  CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
     871 
     872    print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
    860873
    861874    ! wk_adv is the logical flag enabling wake evolution in the time advance
     
    867880      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
    868881                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
     882     
    869883    ENDIF
    870884
    871885    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
    872     ! négatif de ktop à kupper --------
    873     ! cc           On calcule pour cela une densité wdens0 pour laquelle on
     886    ! negatif de ktop a kupper --------
     887    ! cc           On calcule pour cela une densite wdens0 pour laquelle on
    874888    ! aurait un entrainement nul ---
    875889    !jyg<
     
    878892    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
    879893    ! convection profonde cree directement de nouvelles poches, sans passer
    880     ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
     894    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.
    881895
    882896    DO i = 1, klon
     
    894908        IF (wdens(i)<=wdens0*1.1) THEN
    895909          IF (iflag_wk_pop_dyn >= 1) THEN
     910             d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
    896911             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
    897912          ENDIF
     
    908923!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
    909924        sigmaw_targ = min(sigmaw(i), sigmaw_max)
     925        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    910926        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    911927        sigmaw(i) = sigmaw_targ
     
    914930    END DO
    915931
    916     IF (iflag_wk_pop_dyn >= 1) THEN
     932    IF (iflag_wk_pop_dyn == 1) THEN
     933 
     934     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     935                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     936                  d_awdens, d_wdens, d_sigmaw, &
     937                  iflag_wk_act, wk_adv, cin, wape, &
     938                  drdt, &
     939                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     940                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     941                  d_wdens_targ, d_sigmaw_targ)
     942                     
    917943!    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
    918944!    Here, it has to be set to zero.
    919945      death_rate(:) = 0.
    920 
    921       IF (iflag_wk_act ==2) THEN
     946   
     947    ELSEIF (iflag_wk_pop_dyn == 2) THEN
     948     CALL wake_popdyn_2 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     949                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     950                  d_awdens, d_wdens, d_sigmaw, &
     951                  iflag_wk_act, wk_adv, cin, wape, &
     952                  drdt, &
     953                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     954                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     955                  d_wdens_targ, d_sigmaw_targ)
     956     death_rate(:) = 0.
     957   
     958    ELSEIF (iflag_wk_pop_dyn == 0.) THEN
     959   
     960    ! cc nrlmd
     961
    922962      DO i = 1, klon
    923963        IF (wk_adv(i)) THEN
    924           wape1_act(i) = abs(cin(i))
    925           wape2_act(i) = 2.*wape1_act(i) + 1.
    926           act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
    927         ENDIF  ! (wk_adv(i))
    928       ENDDO
    929       ENDIF  ! (iflag_wk_act ==2)
    930 
    931 
    932       DO i = 1, klon
    933         IF (wk_adv(i)) THEN
    934 !!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
    935           tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
    936           tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
    937           drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
    938                     (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
    939 !!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
    940           drdt_pos=max(drdt(i),0.)
    941 
    942 !!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
    943 !!                     - wdens(i)*tau_wk_inv_min &
    944 !!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
    945           d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
    946           d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
    947                          2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
    948           d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
    949 
    950 !!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
    951 !!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
    952 !!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
    953           d_sig_gen(i) = wgen(i)*aa0
    954           d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
    955 !!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
    956           d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
    957           d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
    958           d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
    959         ENDIF
    960       ENDDO
    961 
    962       IF (prt_level >= 10) THEN
    963         print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
    964                        cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
    965         print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
    966                        wdens(1), awdens(1), act(1), d_awdens(1)
    967         print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
    968                        wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
    969         print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
    970                        d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
    971       ENDIF
    972    
    973     ELSE  !  (iflag_wk_pop_dyn >= 1)
    974 
    975     ! cc nrlmd
    976 
    977       DO i = 1, klon
    978         IF (wk_adv(i)) THEN
    979           ! cc nrlmd          Introduction du taux de mortalité des poches et
     964          ! cc nrlmd          Introduction du taux de mortalite des poches et
    980965          ! test sur sigmaw_max=0.4
    981966          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     
    1008993    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
    1009994    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
    1010     ! IM 060208 au niveau k=1..?
     995    ! IM 060208 au niveau k=1...
    1011996    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
    1012997    DO k = 1, klev
     
    10901075
    10911076    ! c      DO i=1,klon
    1092     ! c        print*,'Pente entre 0 et kupper (référence)'
     1077    ! c        print*,'Pente entre 0 et kupper (reference)'
    10931078    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
    10941079    ! c        print*,'Pente entre ktop et kupper'
     
    11661151        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
    11671152          dth(i, k) = deltatw(i, k)/ppi(i, k)
     1153  print *, 'VVVVwake k, the(i,k), dth(i,k), sigmaw(i) ', k, the(i,k), dth(i,k), sigmaw(i)
    11681154          th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
    11691155          th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
     
    12141200    END DO
    12151201
    1216     IF (prt_level>=10) THEN
    1217       PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
    1218       PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
    1219       PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
    1220       PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
     1202!!    IF (prt_level>=10) THEN
     1203    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     1204      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
     1205      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
     1206      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
     1207      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
    12211208    ENDIF
    12221209
     
    13101297
    13111298
    1312           ! Coefficient de répartition
     1299          ! Coefficient de repartition
    13131300
    13141301          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
     
    13401327!
    13411328
    1342           ! cc nrlmd          Prise en compte du taux de mortalité
    1343           ! cc               Définitions de entr, detr
     1329          ! cc nrlmd          Prise en compte du taux de mortalite
     1330          ! cc               Definitions de entr, detr
    13441331!jyg<
    13451332!!            detr(i, k) = 0.
     
    13581345
    13591346
    1360           ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
     1347          ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
    13611348          ! Jingmei
    13621349
     
    14141401    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
    14151402      d_deltaqw, sigmaw, d_sigmaw, alpha)
     1403    !
     1404    ! Alpha_tot = Product of all the alpha's
     1405    DO i = 1, klon
     1406      IF (wk_adv(i)) THEN
     1407        alpha_tot(i) = alpha_tot(i)*alpha(i)   
     1408      END IF
     1409    END DO
    14161410
    14171411    ! cc nrlmd
     
    14781472      DO i = 1, klon
    14791473        IF (wk_adv(i)) THEN
     1474          d_sig_gen2(i)   = d_sig_gen2(i)   + d_sig_gen(i)
     1475          d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
     1476          d_sig_col2(i)   = d_sig_col2(i)   + d_sig_col(i)
     1477          d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
     1478          d_sig_bnd2(i)   = d_sig_bnd2(i)   + d_sig_bnd(i)
     1479        END IF
     1480      END DO
     1481      DO i = 1, klon
     1482        IF (wk_adv(i)) THEN
    14801483          awdens(i) = awdens(i) + d_awdens(i)
    14811484          wdens(i) = wdens(i) + d_wdens(i)
    14821485          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
    14831486          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     1487          d_dens_gen2(i)   = d_dens_gen2(i)   + d_dens_gen(i)
     1488          d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
     1489          d_dens_col2(i)   = d_dens_col2(i)   + d_dens_col(i)
     1490          d_dens_bnd2(i)   = d_dens_bnd2(i)   + d_dens_bnd(i)
    14841491        END IF
    14851492      END DO
     
    14871494        IF (wk_adv(i)) THEN
    14881495          wdens_targ = max(wdens(i),wdensmin)
     1496          d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
    14891497          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
    14901498          wdens(i) = wdens_targ
     
    14981506        IF (wk_adv(i)) THEN
    14991507          sigmaw_targ = max(sigmaw(i),sigmad)
     1508          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    15001509          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    15011510          sigmaw(i) = sigmaw_targ
     
    17441753!!          sigmaw(i) = max(sigmad, sigd_con(i))
    17451754          sigmaw_targ = max(sigmad, sigd_con(i))
     1755          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    17461756          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    17471757          sigmaw(i) = sigmaw_targ
     
    19611971!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
    19621972      sigmaw_targ = max(sigmad, sigd_con(i))
     1973      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
    19631974      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
    19641975      sigmaw(i) = sigmaw_targ
     
    20582069!!      sigmaw(i) = 0.                        !jyg
    20592070      sigmaw_targ = 0.
    2060       d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2071      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     2072!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2073      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
    20612074      sigmaw(i) = sigmaw_targ
    20622075      IF (iflag_wk_pop_dyn >= 1) THEN
     
    20642077!!        wdens(i) = 0.
    20652078        wdens_targ = 0.
     2079        d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
    20662080        d_wdens2(i) = wdens_targ - wdens(i)
    20672081        wdens(i) = wdens_targ
     
    21062120    END DO
    21072121  END DO
    2108 
    2109   DO i = 1, klon
     2122!jyg<
     2123  IF (iflag_wk_pop_dyn >= 1) THEN
     2124  DO i = 1, klon
     2125      IF (ok_qx_qw(i)) THEN
     2126    d_sig_gen2(i) = d_sig_gen2(i)/dtime
     2127    d_sig_death2(i) = d_sig_death2(i)/dtime
     2128    d_sig_col2(i) = d_sig_col2(i)/dtime
     2129    d_sig_spread2(i) = d_sig_spread2(i)/dtime
     2130    d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
    21102131    d_sigmaw2(i) = d_sigmaw2(i)/dtime
     2132!
     2133    d_dens_gen2(i) = d_dens_gen2(i)/dtime
     2134    d_dens_death2(i) = d_dens_death2(i)/dtime
     2135    d_dens_col2(i) = d_dens_col2(i)/dtime
     2136    d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
    21112137    d_awdens2(i) = d_awdens2(i)/dtime
    21122138    d_wdens2(i) = d_wdens2(i)/dtime
     2139      ENDIF
    21132140  ENDDO
     2141  ENDIF
     2142!>jyg
    21142143
    21152144 RETURN
     
    21612190        discrim = b*b - 4.*a*c
    21622191        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
    2163         IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
     2192        IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
    21642193          alpha1(i) = 1.
    21652194        ELSE
     
    21862215  RETURN
    21872216END SUBROUTINE wake_vec_modulation
     2217
     2218
     2219
     2220SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper)
     2221
     2222USE wake_ini_mod , ONLY : pupperbyphs
     2223IMPLICIT NONE
     2224
     2225INTEGER,  INTENT(IN)                              :: klon,klev
     2226REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: ph
     2227REAL,     INTENT(IN),   DIMENSION (klon)          :: ptop
     2228REAL,     INTENT(OUT),  DIMENSION (klon)          :: pupper
     2229INTEGER,  INTENT(OUT),  DIMENSION (klon)          :: kupper
     2230INTEGER                                           :: i,k
     2231
     2232
     2233 kupper = 0
     2234 
     2235IF (pupperbyphs<1.) THEN
     2236 ! Choose an integration bound well above wake top
     2237  ! -----------------------------------------------------------------
     2238
     2239  ! Pupper = 50000.  ! melting level
     2240  ! Pupper = 60000.
     2241  ! Pupper = 80000.  ! essais pour case_e
     2242  DO i = 1, klon
     2243  !  pupper(i) = 0.6*ph(i, 1)
     2244    pupper(i) = pupperbyphs*ph(i, 1)
     2245    pupper(i) = max(pupper(i), 45000.)
     2246    ! cc        Pupper(i) = 60000.
     2247  END DO
     2248
     2249ELSE
     2250 
     2251  DO i=1, klon
     2252     ! pupper(i) = pupperbyphs*ptop(i)+(1.-pupperbyphs)*ph(i, 1)
     2253     pupper(i) = min( pupperbyphs*ptop(i)+(1.-pupperbyphs)*ph(i, 1) , ptop(i)-5000.)
     2254  END DO
     2255END IF
     2256 
     2257  ! -5/ Determination de kupper
     2258
     2259  DO k = klev, 1, -1
     2260    DO i = 1, klon
     2261      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
     2262    END DO
     2263  END DO
     2264
     2265  ! On evite kupper = 1 et kupper = klev
     2266  DO i = 1, klon
     2267    kupper(i) = max(kupper(i), 2)
     2268    kupper(i) = min(kupper(i), klev-1)
     2269  END DO
     2270    RETURN
     2271END SUBROUTINE pkupper
     2272
     2273
     2274SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     2275                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     2276                  d_awdens, d_wdens, d_sigmaw, &
     2277                  iflag_wk_act, wk_adv, cin, wape, &
     2278                  drdt, &
     2279                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     2280                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     2281                  d_wdens_targ, d_sigmaw_targ)
     2282               
     2283
     2284  USE wake_ini_mod , ONLY : wake_ini
     2285  USE wake_ini_mod , ONLY : prt_level,RG
     2286  USE wake_ini_mod , ONLY : stark, wdens_ref
     2287  USE wake_ini_mod , ONLY : tau_cv, rzero, aa0
     2288  USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin
     2289  USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max
     2290 
     2291IMPLICIT NONE
     2292
     2293  INTEGER, INTENT(IN)                                   :: klon,klev
     2294  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
     2295  REAL,                             INTENT(IN)          :: dtime
     2296  REAL,                             INTENT(IN)          :: dtimesub
     2297  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
     2298  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
     2299  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
     2300  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
     2301  REAL, DIMENSION (klon),           INTENT(IN)          :: gfl, cstar
     2302  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
     2303  REAL, DIMENSION (klon),           INTENT(IN)          :: rad_wk
     2304  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
     2305  INTEGER,                          INTENT(IN)          :: iflag_wk_act
     2306
     2307 
     2308  !
     2309 
     2310  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     2311  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     2312  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
     2313  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
     2314  ! Some components of the tendencies of state variables 
     2315  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
     2316  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_spread
     2317  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     2318  REAL,                             INTENT(OUT)         :: d_wdens_targ, d_sigmaw_targ
     2319 
     2320 
     2321  REAL                                                  :: delta_t_min
     2322  INTEGER                                               :: nsub
     2323  INTEGER                                               :: i, k
     2324  REAL                                                  :: wdens0
     2325  ! IM 080208
     2326  LOGICAL, DIMENSION (klon)                             :: gwake
     2327 
     2328   ! Variables liees a la dynamique de population
     2329  REAL, DIMENSION(klon)                                 :: act
     2330  REAL, DIMENSION(klon)                                 :: tau_wk_inv
     2331  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     2332  LOGICAL, DIMENSION (klon)                             :: kill_wake
     2333  REAL                                                  :: drdt_pos
     2334  REAL                                                  :: tau_wk_inv_min
     2335 
     2336     
     2337
     2338      IF (iflag_wk_act == 0) THEN
     2339        act(:) = 0.
     2340      ELSEIF (iflag_wk_act == 1) THEN
     2341        act(:) = 1.
     2342      ELSEIF (iflag_wk_act ==2) THEN
     2343      DO i = 1, klon
     2344        IF (wk_adv(i)) THEN
     2345          wape1_act(i) = abs(cin(i))
     2346          wape2_act(i) = 2.*wape1_act(i) + 1.
     2347          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
     2348        ENDIF  ! (wk_adv(i))
     2349      ENDDO
     2350      ENDIF  ! (iflag_wk_act ==2)
     2351
     2352
     2353      DO i = 1, klon
     2354       ! print*, 'XXX wk_adv(i)', wk_adv(i)
     2355        IF (wk_adv(i)) THEN
     2356!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     2357          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     2358          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     2359          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
     2360                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
     2361!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
     2362          drdt_pos=max(drdt(i),0.)
     2363
     2364!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
     2365!!                     - wdens(i)*tau_wk_inv_min &
     2366!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
     2367!jyg+mlt<
     2368          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
     2369          d_dens_gen(i) = wgen(i)
     2370          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
     2371          d_dens_col(i) =  -2.*wdens(i)*gfl(i)*drdt_pos
     2372          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
     2373          d_dens_death(i) = d_dens_death(i)*dtimesub
     2374          d_dens_col(i) =  d_dens_col(i)*dtimesub
     2375
     2376          d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
     2377!!          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
     2378!!                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
     2379!>jyg+mlt
     2380!
     2381!jyg<
     2382          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
     2383!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     2384          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
     2385          d_wdens(i) = d_wdens_targ
     2386!!          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
     2387!>jyg
     2388
     2389!jyg+mlt<
     2390!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
     2391!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
     2392!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
     2393          d_sig_gen(i) = wgen(i)*aa0
     2394          print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', &
     2395                   sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min
     2396          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
     2397!!       
     2398         
     2399          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
     2400          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
     2401          d_sig_spread(i) = gfl(i)*cstar(i)
     2402          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
     2403          d_sig_death(i) = d_sig_death(i)*dtimesub
     2404          d_sig_col(i) =  d_sig_col(i)*dtimesub
     2405          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
     2406          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
     2407!>jyg+mlt
     2408!
     2409!jyg<
     2410          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
     2411!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     2412!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
     2413          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
     2414          d_sigmaw(i) = d_sigmaw_targ
     2415!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     2416!>jyg
     2417        ENDIF
     2418      ENDDO
     2419
     2420      IF (prt_level >= 10) THEN
     2421        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
     2422                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
     2423        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
     2424                       wdens(1), awdens(1), act(1), d_awdens(1)
     2425        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
     2426                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
     2427        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     2428                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     2429      ENDIF
     2430   
     2431    RETURN
     2432    END SUBROUTINE wake_popdyn_1
     2433   
     2434    SUBROUTINE wake_popdyn_2(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     2435                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     2436                  d_awdens, d_wdens, d_sigmaw, &
     2437                  iflag_wk_act, wk_adv, cin, wape, &
     2438                  drdt, &
     2439                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     2440                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     2441                  d_wdens_targ, d_sigmaw_targ)
     2442               
     2443
     2444  USE wake_ini_mod , ONLY : wake_ini
     2445  USE wake_ini_mod , ONLY : prt_level,RG
     2446  USE wake_ini_mod , ONLY : stark, wdens_ref
     2447  USE wake_ini_mod , ONLY : tau_cv, rzero, aa0
     2448  USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin
     2449  USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max
     2450 
     2451IMPLICIT NONE
     2452
     2453  INTEGER, INTENT(IN)                                   :: klon,klev
     2454  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
     2455  REAL,                             INTENT(IN)          :: dtime
     2456  REAL,                             INTENT(IN)          :: dtimesub
     2457  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
     2458  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
     2459  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
     2460  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
     2461  REAL, DIMENSION (klon),           INTENT(IN)          :: gfl, cstar
     2462  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
     2463  REAL, DIMENSION (klon),           INTENT(IN)          :: rad_wk
     2464  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
     2465  INTEGER,                          INTENT(IN)          :: iflag_wk_act
     2466
     2467 
     2468  !
     2469 
     2470  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     2471  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     2472  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
     2473  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
     2474  ! Some components of the tendencies of state variables 
     2475  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
     2476  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_spread
     2477  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     2478  REAL,                             INTENT(OUT)         :: d_wdens_targ, d_sigmaw_targ
     2479 
     2480 
     2481  REAL                                                  :: delta_t_min
     2482  INTEGER                                               :: nsub
     2483  INTEGER                                               :: i, k
     2484  REAL                                                  :: wdens0
     2485  ! IM 080208
     2486  LOGICAL, DIMENSION (klon)                             :: gwake
     2487 
     2488   ! Variables liees a la dynamique de population
     2489  REAL, DIMENSION(klon)                                 :: act
     2490  REAL, DIMENSION(klon)                                 :: tau_wk_inv
     2491  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     2492  LOGICAL, DIMENSION (klon)                             :: kill_wake
     2493  REAL                                                  :: drdt_pos
     2494  REAL                                                  :: tau_wk_inv_min
     2495 
     2496     
     2497
     2498      IF (iflag_wk_act == 0) THEN
     2499        act(:) = 0.
     2500      ELSEIF (iflag_wk_act == 1) THEN
     2501        act(:) = 1.
     2502      ELSEIF (iflag_wk_act ==2) THEN
     2503      DO i = 1, klon
     2504        IF (wk_adv(i)) THEN
     2505          wape1_act(i) = abs(cin(i))
     2506          wape2_act(i) = 2.*wape1_act(i) + 1.
     2507          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
     2508        ENDIF  ! (wk_adv(i))
     2509      ENDDO
     2510      ENDIF  ! (iflag_wk_act ==2)
     2511
     2512
     2513      DO i = 1, klon
     2514       ! print*, 'XXX wk_adv(i)', wk_adv(i)
     2515        IF (wk_adv(i)) THEN
     2516!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     2517          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     2518          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     2519          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
     2520                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
     2521!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
     2522          drdt_pos=max(drdt(i),0.)
     2523
     2524!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
     2525!!                     - wdens(i)*tau_wk_inv_min &
     2526!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
     2527!jyg+mlt<
     2528          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
     2529          d_dens_gen(i) = wgen(i)
     2530          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
     2531          d_dens_col(i) =  -2.*wdens(i)*gfl(i)*drdt_pos
     2532          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
     2533          d_dens_death(i) = d_dens_death(i)*dtimesub
     2534          d_dens_col(i) =  d_dens_col(i)*dtimesub
     2535
     2536          d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
     2537!!          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
     2538!!                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
     2539!>jyg+mlt
     2540!
     2541!jyg<
     2542          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
     2543!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
     2544          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
     2545          d_wdens(i) = d_wdens_targ
     2546!!          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
     2547!>jyg
     2548
     2549!jyg+mlt<
     2550!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
     2551!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
     2552!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
     2553          d_sig_gen(i) = wgen(i)*aa0
     2554          print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', &
     2555                   sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min
     2556          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
     2557!!       
     2558         
     2559          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
     2560          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
     2561          d_sig_spread(i) = gfl(i)*cstar(i)
     2562          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
     2563          d_sig_death(i) = d_sig_death(i)*dtimesub
     2564          d_sig_col(i) =  d_sig_col(i)*dtimesub
     2565          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
     2566          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
     2567!>jyg+mlt
     2568!
     2569!jyg<
     2570          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
     2571!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
     2572!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
     2573          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
     2574          d_sigmaw(i) = d_sigmaw_targ
     2575!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     2576!>jyg
     2577        ENDIF
     2578      ENDDO
     2579
     2580      IF (prt_level >= 10) THEN
     2581        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
     2582                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
     2583        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
     2584                       wdens(1), awdens(1), act(1), d_awdens(1)
     2585        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
     2586                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
     2587        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     2588                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     2589      ENDIF
     2590   
     2591    RETURN
     2592    END SUBROUTINE wake_popdyn_2 
     2593 
     2594 
  • LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90

    r4085 r4230  
    1717
    1818  ! -------------------------------------------------------------------------
    19   ! Déclaration de variables
    20   ! -------------------------------------------------------------------------
    21 
    22   ! Variables à fixer
     19  ! Declaration de variables
     20  ! -------------------------------------------------------------------------
     21
     22  ! Variables a fixer
    2323!jyg<
    2424!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
    2525  INTEGER, SAVE, PROTECTED                                         :: prt_level
    2626  REAL, SAVE, PROTECTED, DIMENSION(2)                              :: wdens_ref
    27   REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk
     27  REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk, pupperbyphs
    2828!>jyg
    2929  REAL, SAVE, PROTECTED                                            :: crep_upper, crep_sol 
    30   !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
     30  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, pupperbyphs, crep_upper, crep_sol)
    3131
    3232  REAL, SAVE, PROTECTED                                            :: tau_cv
     
    111111  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
    112112
    113   ! coefgw : Coefficient pour les ondes de gravité
     113  ! coefgw : Coefficient pour les ondes de gravite
    114114  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    115   ! wdens : Densité surfacique de poche froide
     115  ! wdens : Densite surfacique de poche froide
    116116  ! -------------------------------------------------------------------------
    117117
     
    149149  alpk=0.25
    150150  CALL getin_p('alpk',alpk)
     151 
     152  pupperbyphs=0.6
     153  CALL getin_p('pupperbyphs',pupperbyphs)
     154
     155
    151156!jyg<
    152157!!  wdens_ref=8.E-12
     
    184189  WRITE(*,*) 'stark=', stark
    185190  WRITE(*,*) 'alpk=', alpk
     191  WRITE(*,*) 'pupperbyphs=', pupperbyphs
    186192!jyg<
    187193!!  WRITE(*,*) 'wdens_ref=', wdens_ref
Note: See TracChangeset for help on using the changeset viewer.