Changeset 2384


Ignore:
Timestamp:
Nov 5, 2015, 6:49:55 PM (9 years ago)
Author:
fhourdin
Message:

Séparation de la partie de calcul de ALP dans thermcell_main
Spliting thermcell_main.F90 in 2, isolating the ALP computation.

Location:
LMDZ5/trunk/libf/phylmd
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/thermcell_main.F90

    r2351 r2384  
    747747      enddo
    748748!
    749 
     749! $Id$
     750!
     751      print*,'THERM_ALP sorti'
     752      CALL thermcell_alp(ngrid,nlay,ptimestep  &
     753     &                  ,pplay,pplev  &
     754     &                  ,fm0,entr0,lmax  &
     755     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
     756     &                  ,zw2,fraca &
     757!!! necessire en plus
     758     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
    750759!!! nrlmd le 10/04/2012
    751 
    752 !------------Test sur le LCL des thermiques
    753     do ig=1,ngrid
    754       ok_lcl(ig)=.false.
    755       if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
    756     enddo
    757 
    758 !------------Localisation des niveaux entourant le LCL et du coef d'interpolation
    759     do l=1,nlay-1
    760       do ig=1,ngrid
    761         if (ok_lcl(ig)) then
    762 !ATTENTION,zw2 calcule en pplev
    763 !          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
    764 !          klcl(ig)=l
    765 !          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
    766 !          endif
    767           if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
    768           klcl(ig)=l
    769           interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
    770           endif
    771         endif
    772       enddo
    773     enddo
    774 
    775 !------------Hauteur des thermiques
    776 !!jyg le 27/04/2012
    777 !!    do ig =1,ngrid
    778 !!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
    779 !! &               -rhobarz(ig,klcl(ig)))*interp(ig)
    780 !!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
    781 !!    zmax(ig)=pphi(ig,lmax(ig))/rg
    782 !!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
    783 !!    enddo
    784     do ig =1,ngrid
    785 !CR:REHABILITATION ZMAX CONTINU
    786 !     zmax(ig)=pphi(ig,lmax(ig))/rg
    787      if (ok_lcl(ig)) then
    788       rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
    789  &               -rhobarz(ig,klcl(ig)))*interp(ig)
    790       zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
    791       zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
    792      else
    793       rhobarz0(ig)=0.
    794       zlcl(ig)=zmax(ig)
    795      endif
    796     enddo
    797 !!jyg fin
    798 
    799 !------------Calcul des propriétés du thermique au LCL
    800   IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
    801 
    802   !-----Initialisation de la TKE moyenne
    803    do l=1,nlay
    804     do ig=1,ngrid
    805      pbl_tke_max(ig,l)=0.
    806     enddo
    807    enddo
    808 
    809 !-----Calcul de la TKE moyenne
    810    do nsrf=1,nbsrf
    811     do l=1,nlay
    812      do ig=1,ngrid
    813      pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
    814      enddo
    815     enddo
    816    enddo
    817 
    818 !-----Initialisations des TKE dans et hors des thermiques
    819    do l=1,nlay
    820     do ig=1,ngrid
    821     therm_tke_max(ig,l)=pbl_tke_max(ig,l)
    822     env_tke_max(ig,l)=pbl_tke_max(ig,l)
    823     enddo
    824    enddo
    825 
    826 !-----Calcul de la TKE transportée par les thermiques : therm_tke_max
    827    call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
    828   &           rg,pplev,therm_tke_max)
    829 !   print *,' thermcell_tke_transport -> '   !!jyg
    830 
    831 !-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
    832    do l=1,nlay
    833     do ig=1,ngrid
    834      pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l)         !  Recalcul de TKE moyenne aprés transport de TKE_TH
    835      env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l))       !  Recalcul de TKE dans  l'environnement aprés transport de TKE_TH
    836      w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
    837     enddo
    838    enddo
    839 !    print *,' apres w_ls = '   !!jyg
    840 
    841   do ig=1,ngrid
    842    if (ok_lcl(ig)) then
    843      fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
    844  &             -fraca(ig,klcl(ig)))*interp(ig)
    845      w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
    846  &         -zw2(ig,klcl(ig)))*interp(ig)
    847      w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
    848  &             -w_ls(ig,klcl(ig)))*interp(ig)
    849      therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
    850  &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
    851      env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
    852  &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
    853      pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
    854  &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
    855      if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
    856      if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
    857      if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
    858    else
    859      fraca0(ig)=0.
    860      w0(ig)=0.
    861 !!jyg le 27/04/2012
    862 !!     zlcl(ig)=0.
    863 !!
    864    endif
    865   enddo
    866 
    867   ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
    868 !  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
    869 
    870 !------------Triggering------------------
    871   IF (iflag_trig_bl.ge.1) THEN
    872 
    873 !-----Initialisations
    874    depth(:)=0.
    875    n2(:)=0.
    876    s2(:)=100. ! some low value, arbitrary
    877    s_max(:)=0.
    878 
    879 !-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
    880    do ig=1,ngrid
    881      zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
    882      depth(ig)=zmax_moy(ig)-zlcl(ig)
    883      hmin(ig)=hmincoef*zlcl(ig)
    884      if (depth(ig).ge.10.) then
    885        s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
    886        n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
    887 !!
    888 !!jyg le 27/04/2012
    889 !!       s_max(ig)=s2(ig)*log(n2(ig))
    890 !!       if (n2(ig) .lt. 1) s_max(ig)=0.
    891        s_max(ig)=s2(ig)*log(max(n2(ig),1.))
    892 !!fin jyg
    893      else
    894        n2(ig)=0.
    895        s_max(ig)=0.
    896      endif
    897    enddo
    898 !   print *,'avant Calcul de Wmax '    !!jyg
    899 
    900 !-----Calcul de Wmax et ALE_BL_STAT associée
    901 !!jyg le 30/04/2012
    902 !!   do ig=1,ngrid
    903 !!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
    904 !!     w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14))))
    905 !!     ale_bl_stat(ig)=0.5*w_max(ig)**2
    906 !!     else
    907 !!     w_max(ig)=0.
    908 !!     ale_bl_stat(ig)=0.
    909 !!     endif
    910 !!   enddo
    911    susqr2pi=su*sqrt(2.*Rpi)
    912    Reuler=exp(1.)
    913    do ig=1,ngrid
    914      if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
    915       w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
    916       ale_bl_stat(ig)=0.5*w_max(ig)**2
    917      else
    918       w_max(ig)=0.
    919       ale_bl_stat(ig)=0.
    920      endif
    921    enddo
    922 
    923   ENDIF ! iflag_trig_bl
    924 !  print *,'ENDIF  iflag_trig_bl'    !!jyg
    925 
    926 !------------Closure------------------
    927 
    928   IF (iflag_clos_bl.ge.2) THEN
    929 
    930 !-----Calcul de ALP_BL_STAT
    931   do ig=1,ngrid
    932   alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
    933   alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
    934  &                   (w0(ig)**2)
    935   alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
    936  &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
    937     if (iflag_clos_bl.ge.2) then
    938     alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
    939  &                   (w0(ig)**2)
    940     else
    941     alp_bl_conv(ig)=0.
    942     endif
    943   alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
    944   enddo
    945 
    946 !-----Sécurité ALP infinie
    947   do ig=1,ngrid
    948    if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
    949   enddo
    950 
    951   ENDIF ! (iflag_clos_bl.ge.2)
    952 
     760     &                  ,pbl_tke,pctsrf,omega,airephy &
     761     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
     762     &                  ,n2,s2,ale_bl_stat &
     763     &                  ,therm_tke_max,env_tke_max &
     764     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
     765     &                  ,alp_bl_conv,alp_bl_stat &
    953766!!! fin nrlmd le 10/04/2012
    954 
    955       if (prt_level.ge.10) then
    956          ig=igout
    957          do l=1,nlay
    958             print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
    959             print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
    960          enddo
    961       endif
    962 
    963 !      print*,'avant calcul ale et alp'
    964 !calcul de ALE et ALP pour la convection
    965       Alp_bl(:)=0.
    966       Ale_bl(:)=0.
    967 !          print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
    968       do l=1,nlay
    969       do ig=1,ngrid
    970            Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
    971            Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2)
    972 !          print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
    973       enddo
    974       enddo
    975 
    976 ! Ale sec (max de wmax/2 sous la zone d'inhibition) dans
    977 ! le cas iflag_trig_bl=3
    978       IF (iflag_trig_bl==3) Ale_bl(:)=0.5*wmax_sec(:)**2
    979 
    980 !test:calcul de la ponderation des couches pour KE
    981 !initialisations
    982 
    983       fm_tot(:)=0.
    984       wght_th(:,:)=1.
    985       lalim_conv(:)=lalim(:)
    986 
    987       do k=1,klev
    988          do ig=1,ngrid
    989             if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
    990          enddo
    991       enddo
    992 
    993 ! assez bizarre car, si on est dans la couche d'alim et que alim_star et
    994 ! plus petit que 1.e-10, on prend wght_th=1.
    995       do k=1,klev
    996          do ig=1,ngrid
    997             if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
    998                wght_th(ig,k)=alim_star(ig,k)
    999             endif
    1000          enddo
    1001       enddo
    1002 
    1003 !      print*,'apres wght_th'
    1004 !test pour prolonger la convection
    1005       do ig=1,ngrid
    1006 !v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
    1007       if ((alim_star(ig,1).lt.1.e-10)) then
    1008       lalim_conv(ig)=1
    1009       wght_th(ig,1)=1.
    1010 !      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
    1011       endif
    1012       enddo
    1013 
    1014 !------------------------------------------------------------------------
    1015 ! Modif CR/FH 20110310 : Alp integree sur la verticale.
    1016 ! Integrale verticale de ALP.
    1017 ! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
    1018 ! couches
    1019 !------------------------------------------------------------------------
    1020 
    1021       alp_int(:)=0.
    1022       dp_int(:)=0.
    1023       do l=2,nlay
    1024         do ig=1,ngrid
    1025            if(l.LE.lmax(ig)) THEN
    1026            zdp=pplay(ig,l-1)-pplay(ig,l)
    1027            alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
    1028            dp_int(ig)=dp_int(ig)+zdp
    1029            endif
    1030         enddo
    1031       enddo
    1032 
    1033       if (iflag_coupl>=3 .and. iflag_coupl<=5) then
    1034       do ig=1,ngrid
    1035 !valeur integree de alp_bl * 0.5:
    1036         if (dp_int(ig)>0.) then
    1037         Alp_bl(ig)=alp_int(ig)/dp_int(ig)
    1038         endif
    1039       enddo!
    1040       endif
    1041 
    1042 
    1043 ! Facteur multiplicatif sur Alp_bl
    1044       Alp_bl(:)=alp_bl_k*Alp_bl(:)
    1045 
    1046 !------------------------------------------------------------------------
     767     &                   )
     768
    1047769
    1048770
Note: See TracChangeset for help on using the changeset viewer.