Ignore:
Timestamp:
Oct 16, 2012, 2:41:50 PM (12 years ago)
Author:
Laurent Fairhead
Message:

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires


Testing release based on r1668

Location:
LMDZ5/branches/testing
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/calltherm.F90

    r1517 r1669  
    88     &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
    99     &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
    10      &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl)
     10     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl &
     11!!! nrlmd le 10/04/2012
     12     &      ,pbl_tke,pctsrf,omega,airephy &
     13     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
     14     &      ,n2,s2,ale_bl_stat &
     15     &      ,therm_tke_max,env_tke_max &
     16     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
     17     &      ,alp_bl_conv,alp_bl_stat &
     18!!! fin nrlmd le 10/04/2012
     19     &                    )
    1120
    1221      USE dimphy
     
    1625#include "thermcell.h"
    1726#include "iniprint.h"
     27
     28!!! nrlmd le 10/04/2012
     29#include "indicesol.h"
     30!!! fin nrlmd le 10/04/2012
    1831
    1932!IM 140508
     
    7588      !on garde le zmax du pas de temps precedent
    7689      real zmax0(klon), f0(klon)
     90
     91!!! nrlmd le 10/04/2012
     92      real pbl_tke(klon,klev+1,nbsrf)
     93      real pctsrf(klon,nbsrf)
     94      real omega(klon,klev)
     95      real airephy(klon)
     96      real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
     97      real therm_tke_max0(klon),env_tke_max0(klon)
     98      real n2(klon),s2(klon)
     99      real ale_bl_stat(klon)
     100      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
     101      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
     102!!! fin nrlmd le 10/04/2012
     103
    77104!********************************************************
    78105
     
    220247     &      ,Ale,Alp,lalim_conv,wght_th &
    221248     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
    222      &      ,ztla,zthl)
     249     &      ,ztla,zthl &
     250!!! nrlmd le 10/04/2012
     251     &      ,pbl_tke,pctsrf,omega,airephy &
     252     &      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
     253     &      ,n2,s2,ale_bl_stat &
     254     &      ,therm_tke_max,env_tke_max &
     255     &      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
     256     &      ,alp_bl_conv,alp_bl_stat &
     257!!! fin nrlmd le 10/04/2012
     258     &                         )
    223259           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
    224260         else
     
    227263         endif
    228264
    229        flag_bidouille_stratocu=iflag_thermals.eq.14.or.iflag_thermals.eq.16
     265! Attention : les noms sont contre intuitif.
     266! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
     267! Il aurait mieux valu avoir un nobidouille_stratocu
     268! Et pour simplifier :
     269! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
     270! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
     271! fait bien ce qu'on croit.
     272
     273       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16
     274
     275      if (iflag_thermals<=12) then
     276         lmax=1
     277         do k=1,klev-1
     278            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
     279         enddo
     280      endif
    230281
    231282      fact(:)=0.
     
    267318
    268319       DO i=1,klon
    269         if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
    270320            fm_therm(i,klev+1)=0.
    271321            Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
     
    273323            Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
    274324!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
     325        if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
    275326       ENDDO
    276327
  • LMDZ5/branches/testing/libf/phylmd/concvl.F

    r1665 r1669  
    248248         DO i = 1, klon
    249249          cbmf(i) = 0.
    250 !          plcl(i) = 0.
     250          plcl(i) = 0.
    251251          sigd(i) = 0.
    252252         ENDDO
     
    256256      plfc(:)  = 0.
    257257      wbeff(:) = 100.
    258       plcl(:) = 0.
    259258
    260259      DO k = 1, klev+1
     
    369368     $              cape,cin,tvp,
    370369     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
    371      $              asupmaxmin,lalim_conv)
     370     $              asupmaxmin,lalim_conv,
     371!AC!
     372     $              da,phi)
     373!AC!
    372374      endif 
    373375C------------------------------------------------------------------
     
    399401       ENDDO
    400402       endif
     403
     404c!AC!
     405       if (iflag_con.eq.3) then
     406       DO itra = 1,ntra
     407        DO k = 1, klev
     408         DO i = 1, klon
     409            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
     410         ENDDO
     411        ENDDO
     412       ENDDO
     413       endif
     414c!AC!
    401415
    402416      DO k = 1, klev
  • LMDZ5/branches/testing/libf/phylmd/conf_phys.F90

    r1664 r1669  
    110110  integer,SAVE :: iflag_thermals_omp,nsplit_thermals_omp
    111111  real,save :: tau_thermals_omp,alp_bl_k_omp
     112!!! nrlmd le 10/04/2012
     113  integer,SAVE :: iflag_trig_bl_omp,iflag_clos_bl_omp
     114  integer,SAVE :: tau_trig_shallow_omp,tau_trig_deep_omp
     115  real,SAVE    :: s_trig_omp
     116!!! fin nrlmd le 10/04/2012
    112117  real :: alp_offset
    113118  REAL, SAVE :: alp_offset_omp
     
    282287!Config Help = Used in physiq.F
    283288!
     289! - flag_aerosol=0 => no aerosol
    284290! - flag_aerosol=1 => so4 only (defaut)
    285291! - flag_aerosol=2 => bc  only
     
    289295! - flag_aerosol=6 => all aerosol
    290296
    291   flag_aerosol_omp = 1
     297  flag_aerosol_omp = 0
    292298  CALL getin('flag_aerosol',flag_aerosol_omp)
    293299
     
    10831089  alp_bl_k_omp = 1.
    10841090  call getin('alp_bl_k',alp_bl_k_omp)
     1091
     1092!!! nrlmd le 10/04/2012
     1093
     1094!Config Key  = iflag_trig_bl
     1095!Config Desc = 
     1096!Config Def  = 0
     1097!Config Help =
     1098!
     1099  iflag_trig_bl_omp = 0
     1100  call getin('iflag_trig_bl',iflag_trig_bl_omp)
     1101
     1102!Config Key  = s_trig_bl
     1103!Config Desc = 
     1104!Config Def  = 0
     1105!Config Help =
     1106!
     1107  s_trig_omp = 2e7
     1108  call getin('s_trig',s_trig_omp)
     1109
     1110!Config Key  = tau_trig_shallow
     1111!Config Desc = 
     1112!Config Def  = 0
     1113!Config Help =
     1114!
     1115  tau_trig_shallow_omp = 600
     1116  call getin('tau_trig_shallow',tau_trig_shallow_omp)
     1117
     1118!Config Key  = tau_trig_deep
     1119!Config Desc = 
     1120!Config Def  = 0
     1121!Config Help =
     1122!
     1123  tau_trig_deep_omp = 1800
     1124  call getin('tau_trig_deep',tau_trig_deep_omp)
     1125
     1126!Config Key  = iflag_clos_bl
     1127!Config Desc = 
     1128!Config Def  = 0
     1129!Config Help =
     1130!
     1131  iflag_clos_bl_omp = 0
     1132  call getin('iflag_clos_bl',iflag_clos_bl_omp)
     1133
     1134!!! fin nrlmd le 10/04/2012
    10851135
    10861136!
     
    16501700    tau_thermals = tau_thermals_omp
    16511701    alp_bl_k = alp_bl_k_omp
     1702!!! nrlmd le 10/04/2012
     1703    iflag_trig_bl = iflag_trig_bl_omp
     1704    s_trig = s_trig_omp
     1705    tau_trig_shallow = tau_trig_shallow_omp
     1706    tau_trig_deep = tau_trig_deep_omp
     1707    iflag_clos_bl = iflag_clos_bl_omp
     1708!!! fin nrlmd le 10/04/2012
    16521709    iflag_coupl = iflag_coupl_omp
    16531710    iflag_clos = iflag_clos_omp
     
    17101767! il n'est utilisable que lors du couplage avec le SO4 seul
    17111768    IF (ok_ade .OR. ok_aie) THEN
     1769       IF ( flag_aerosol .EQ. 0 ) THEN
     1770          CALL abort_gcm('conf_phys','flag_aerosol=0 not compatible avec ok_ade ou ok_aie=.TRUE.',1)
     1771       END IF
    17121772       IF ( .NOT. new_aod .AND.  flag_aerosol .NE. 1) THEN
    17131773          CALL abort_gcm('conf_phys','new_aod=.FALSE. not compatible avec flag_aerosol=1',1)
     
    18391899  write(lunout,*)' iflag_wake = ', iflag_wake
    18401900  write(lunout,*)' alp_offset = ', alp_offset
     1901!!! nrlmd le 10/04/2012
     1902  write(lunout,*)' iflag_trig_bl = ', iflag_trig_bl
     1903  write(lunout,*)' s_trig = ', s_trig
     1904  write(lunout,*)' tau_trig_shallow = ', tau_trig_shallow
     1905  write(lunout,*)' tau_trig_deep = ', tau_trig_deep
     1906  write(lunout,*)' iflag_clos_bl = ', iflag_clos_bl
     1907!!! fin nrlmd le 10/04/2012
    18411908
    18421909  write(lunout,*)' lonmin lonmax latmin latmax bilKP_ins =',&
  • LMDZ5/branches/testing/libf/phylmd/cv3_routines.F

    r1554 r1669  
    879879 110  continue
    880880
    881       do 121 j=1,ntra
    882 ccccc      do 111 k=1,nl+1
    883       do 111 k=1,nd
    884        nn=0
    885       do 101 i=1,len
    886       if(iflag1(i).eq.0)then
    887        nn=nn+1
    888        tra(nn,k,j)=tra1(i,k,j)
    889       endif
    890  101  continue
    891  111  continue
    892  121  continue
     881!AC!      do 121 j=1,ntra
     882!AC!ccccc      do 111 k=1,nl+1
     883!AC!      do 111 k=1,nd
     884!AC!       nn=0
     885!AC!      do 101 i=1,len
     886!AC!      if(iflag1(i).eq.0)then
     887!AC!       nn=nn+1
     888!AC!       tra(nn,k,j)=tra1(i,k,j)
     889!AC!      endif
     890!AC! 101  continue
     891!AC! 111  continue
     892!AC! 121  continue
    893893
    894894      if (nn.ne.ncum) then
     
    16331633      sij(1:ncum,1:nd,1:nd)=0.0
    16341634     
    1635       do k=1,ntra
    1636        do j=1,nd  ! instead nlp
    1637         do i=1,nd ! instead nlp
    1638          do il=1,ncum
    1639             traent(il,i,j,k)=tra(il,j,k)
    1640          enddo
    1641         enddo
    1642        enddo
    1643       enddo
     1635!AC!      do k=1,ntra
     1636!AC!       do j=1,nd  ! instead nlp
     1637!AC!        do i=1,nd ! instead nlp
     1638!AC!         do il=1,ncum
     1639!AC!            traent(il,i,j,k)=tra(il,j,k)
     1640!AC!         enddo
     1641!AC!        enddo
     1642!AC!       enddo
     1643!AC!      enddo
    16441644      zm(:,:)=0.
    16451645
     
    16971697 710  continue
    16981698
    1699        do k=1,ntra
    1700         do j=minorig,nl
    1701          do il=1,ncum
    1702           if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
    1703      :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
    1704             traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
    1705      :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
    1706           endif
    1707          enddo
    1708         enddo
    1709        enddo
     1699!AC!       do k=1,ntra
     1700!AC!        do j=minorig,nl
     1701!AC!         do il=1,ncum
     1702!AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
     1703!AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
     1704!AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1705!AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1706!AC!          endif
     1707!AC!         enddo
     1708!AC!        enddo
     1709!AC!       enddo
    17101710
    17111711c
     
    17301730 750  continue
    17311731
    1732       do j=1,ntra
    1733        do i=minorig+1,nl
    1734         do il=1,ncum
    1735          if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
    1736           traent(il,i,i,j)=tra(il,nk(il),j)
    1737          endif
    1738         enddo
    1739        enddo
    1740       enddo
     1732!AC!      do j=1,ntra
     1733!AC!       do i=minorig+1,nl
     1734!AC!        do il=1,ncum
     1735!AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
     1736!AC!          traent(il,i,i,j)=tra(il,nk(il),j)
     1737!AC!         endif
     1738!AC!        enddo
     1739!AC!       enddo
     1740!AC!      enddo
    17411741
    17421742      do 100 j=minorig,nl
     
    19041904      enddo ! il
    19051905
    1906       do j=1,ntra
    1907        do il=1,ncum
    1908         if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
    1909      :     .and. csum(il,i).lt.m(il,i) ) then
    1910          traent(il,i,i,j)=tra(il,nk(il),j)
    1911         endif
    1912        enddo
    1913       enddo
     1906!AC!      do j=1,ntra
     1907!AC!       do il=1,ncum
     1908!AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
     1909!AC!     :     .and. csum(il,i).lt.m(il,i) ) then
     1910!AC!         traent(il,i,i,j)=tra(il,nk(il),j)
     1911!AC!        endif
     1912!AC!       enddo
     1913!AC!      enddo
    19141914789   continue
    19151915c     
     
    20142014         enddo
    20152015        enddo
    2016         do k=1,ntra
    2017          do i=1,nd
    2018           do il=1,ncum
    2019            trap(il,i,k)=tra(il,i,k)
    2020           enddo
    2021          enddo
    2022         enddo
     2016!AC!        do k=1,ntra
     2017!AC!         do i=1,nd
     2018!AC!          do il=1,ncum
     2019!AC!           trap(il,i,k)=tra(il,i,k)
     2020!AC!          enddo
     2021!AC!         enddo
     2022!AC!        enddo
    20232023c
    20242024c   ***  check whether ep(inb)=0, if so, skip precipitating    ***
     
    23412341c    ***       find tracer concentrations in precipitating downdraft     ***
    23422342c
    2343       do j=1,ntra
    2344        do il = 1,ncum
    2345        if (i.lt.inb(il) .and. lwork(il)) then
    2346 c
    2347          if(mplus(il))then
    2348           trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
    2349      :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
    2350           trap(il,i,j)=trap(il,i,j)/mp(il,i)
    2351          else ! if (mplus(il))
    2352           if(mp(il,i+1).gt.1.0e-16)then
    2353            trap(il,i,j)=trap(il,i+1,j)
    2354           endif
    2355          endif ! (mplus(il)) else if (.not.mplus(il))
    2356 c
    2357         endif ! (i.lt.inb(il) .and. lwork(il))
    2358        enddo
    2359       end do
     2343!AC!      do j=1,ntra
     2344!AC!       do il = 1,ncum
     2345!AC!       if (i.lt.inb(il) .and. lwork(il)) then
     2346!AC!c
     2347!AC!         if(mplus(il))then
     2348!AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
     2349!AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
     2350!AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
     2351!AC!         else ! if (mplus(il))
     2352!AC!          if(mp(il,i+1).gt.1.0e-16)then
     2353!AC!           trap(il,i,j)=trap(il,i+1,j)
     2354!AC!          endif
     2355!AC!         endif ! (mplus(il)) else if (.not.mplus(il))
     2356!AC!c
     2357!AC!        endif ! (i.lt.inb(il) .and. lwork(il))
     2358!AC!       enddo
     2359!AC!      end do
    23602360
    23612361400   continue
     
    24842484      enddo
    24852485c       print*,'cv3_yield initialisation 2'
    2486       do j=1,ntra
    2487        do i=1,nd
    2488         do il=1,ncum
    2489           ftra(il,i,j)=0.0
    2490         enddo
    2491        enddo
    2492       enddo
     2486!AC!      do j=1,ntra
     2487!AC!       do i=1,nd
     2488!AC!        do il=1,ncum
     2489!AC!          ftra(il,i,j)=0.0
     2490!AC!        enddo
     2491!AC!       enddo
     2492!AC!      enddo
    24932493c       print*,'cv3_yield initialisation 3'
    24942494      do i=1,nl
     
    26492649
    26502650
    2651       do j=1,ntra
    2652        do il=1,ncum
    2653         if (iflag(il) .le. 1) then
    2654         if (cvflag_grav) then
    2655          ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
    2656      :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2657      :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2658         else
    2659          ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
    2660      :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2661      :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2662         endif
    2663         endif  ! iflag
    2664        enddo
    2665       enddo
     2651!AC!     do j=1,ntra
     2652!AC!      do il=1,ncum
     2653!AC!       if (iflag(il) .le. 1) then
     2654!AC!       if (cvflag_grav) then
     2655!AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
     2656!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     2657!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     2658!AC!       else
     2659!AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
     2660!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     2661!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     2662!AC!       endif
     2663!AC!       endif  ! iflag
     2664!AC!      enddo
     2665!AC!     enddo
    26662666
    26672667       do j=2,nl
     
    26872687      enddo
    26882688
    2689       do k=1,ntra
    2690        do j=2,nl
    2691         do il=1,ncum
    2692          if (j.le.inb(il) .and. iflag(il) .le. 1) then
    2693 
    2694           if (cvflag_grav) then
    2695            ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
    2696      :                *(traent(il,j,1,k)-tra(il,1,k))
    2697           else
    2698            ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
    2699      :                *(traent(il,j,1,k)-tra(il,1,k))
    2700           endif
    2701 
    2702          endif
    2703         enddo
    2704        enddo
    2705       enddo
     2689!AC!      do k=1,ntra
     2690!AC!       do j=2,nl
     2691!AC!        do il=1,ncum
     2692!AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
     2693!AC!
     2694!AC!          if (cvflag_grav) then
     2695!AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
     2696!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
     2697!AC!          else
     2698!AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
     2699!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
     2700!AC!          endif
     2701!AC!
     2702!AC!         endif
     2703!AC!        enddo
     2704!AC!       enddo
     2705!AC!      enddo
    27062706c      print*,'cv3_yield apres ft'
    27072707c
     
    286528651350  continue
    28662866
    2867       do k=1,ntra
    2868        do il=1,ncum
    2869         if (i.le.inb(il) .and. iflag(il) .le. 1) then
    2870          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2871          cpinv=1.0/cpn(il,i)
    2872          if (cvflag_grav) then
    2873            ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
    2874      :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    2875      :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    2876          else
    2877            ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
    2878      :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    2879      :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    2880          endif
    2881         endif
    2882        enddo
    2883       enddo
     2867!AC!      do k=1,ntra
     2868!AC!       do il=1,ncum
     2869!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     2870!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2871!AC!         cpinv=1.0/cpn(il,i)
     2872!AC!         if (cvflag_grav) then
     2873!AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
     2874!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     2875!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     2876!AC!         else
     2877!AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
     2878!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     2879!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     2880!AC!         endif
     2881!AC!        endif
     2882!AC!       enddo
     2883!AC!      enddo
    28842884
    28852885      do 480 k=1,i-1
     
    29382938480   continue
    29392939
    2940       do j=1,ntra
    2941        do k=1,i-1
    2942         do il=1,ncum
    2943          if (i.le.inb(il) .and. iflag(il) .le. 1) then
    2944           dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2945           cpinv=1.0/cpn(il,i)
    2946           if (cvflag_grav) then
    2947            ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    2948      :        *(traent(il,k,i,j)-tra(il,i,j))
    2949           else
    2950            ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    2951      :        *(traent(il,k,i,j)-tra(il,i,j))
    2952           endif
    2953          endif
    2954         enddo
    2955        enddo
    2956       enddo
     2940!AC!      do j=1,ntra
     2941!AC!       do k=1,i-1
     2942!AC!        do il=1,ncum
     2943!AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
     2944!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2945!AC!          cpinv=1.0/cpn(il,i)
     2946!AC!          if (cvflag_grav) then
     2947!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     2948!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
     2949!AC!          else
     2950!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     2951!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
     2952!AC!          endif
     2953!AC!         endif
     2954!AC!        enddo
     2955!AC!       enddo
     2956!AC!      enddo
    29572957
    29582958      do 490 k=i,nl+1
     
    30043004490   continue
    30053005
    3006       do j=1,ntra
    3007        do k=i,nl+1
    3008         do il=1,ncum
    3009          if (i.le.inb(il) .and. k.le.inb(il)
    3010      $                .and. iflag(il) .le. 1) then
    3011           dpinv=1.0/(ph(il,i)-ph(il,i+1))
    3012           cpinv=1.0/cpn(il,i)
    3013           if (cvflag_grav) then
    3014            ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    3015      :         *(traent(il,k,i,j)-tra(il,i,j))
    3016           else
    3017            ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    3018      :             *(traent(il,k,i,j)-tra(il,i,j))
    3019           endif
    3020          endif ! i and k
    3021         enddo
    3022        enddo
    3023       enddo
     3006!AC!      do j=1,ntra
     3007!AC!       do k=i,nl+1
     3008!AC!        do il=1,ncum
     3009!AC!         if (i.le.inb(il) .and. k.le.inb(il)
     3010!AC!     $                .and. iflag(il) .le. 1) then
     3011!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3012!AC!          cpinv=1.0/cpn(il,i)
     3013!AC!          if (cvflag_grav) then
     3014!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     3015!AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
     3016!AC!          else
     3017!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     3018!AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
     3019!AC!          endif
     3020!AC!         endif ! i and k
     3021!AC!        enddo
     3022!AC!       enddo
     3023!AC!      enddo
    30243024
    30253025c sb: interface with the cloud parameterization:          ! cld
     
    30523052      enddo
    30533053
    3054       do j=1,ntra
    3055        do il=1,ncum
    3056         if (i.le.inb(il) .and. iflag(il) .le. 1) then
    3057          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    3058          cpinv=1.0/cpn(il,i)
    3059 
    3060          if (cvflag_grav) then
    3061           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
    3062      :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    3063      :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
    3064          else
    3065           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
    3066      :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    3067      :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
    3068          endif
    3069         endif ! i
    3070        enddo
    3071       enddo
     3054!AC!      do j=1,ntra
     3055!AC!       do il=1,ncum
     3056!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3057!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3058!AC!         cpinv=1.0/cpn(il,i)
     3059!AC!
     3060!AC!         if (cvflag_grav) then
     3061!AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
     3062!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     3063!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     3064!AC!         else
     3065!AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
     3066!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     3067!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     3068!AC!         endif
     3069!AC!        endif ! i
     3070!AC!       enddo
     3071!AC!      enddo
    30723072
    30733073
     
    31463146503   continue
    31473147
    3148       do j=1,ntra
    3149        do il=1,ncum
    3150         IF (iflag(il) .le. 1) THEN
    3151         IF (cvflag_grav) then
    3152         ex=0.01*grav*ment(il,inb(il),inb(il))
    3153      :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    3154      :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
    3155         ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
    3156         ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
    3157      :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    3158      :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
    3159         else
    3160         ex=0.1*ment(il,inb(il),inb(il))
    3161      :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    3162      :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
    3163         ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
    3164         ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
    3165      :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    3166      :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
    3167         ENDIF   !cvflag grav
    3168         ENDIF    !iflag
    3169        enddo
    3170       enddo
     3148!AC!      do j=1,ntra
     3149!AC!       do il=1,ncum
     3150!AC!        IF (iflag(il) .le. 1) THEN
     3151!AC!    IF (cvflag_grav) then
     3152!AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
     3153!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3154!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3155!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3156!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3157!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3158!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3159!AC!    else
     3160!AC!        ex=0.1*ment(il,inb(il),inb(il))
     3161!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3162!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3163!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3164!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3165!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3166!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3167!AC!        ENDIF   !cvflag grav
     3168!AC!        ENDIF    !iflag
     3169!AC!       enddo
     3170!AC!      enddo
    31713171
    31723172c
     
    32873287      ENDDO
    32883288      ENDDO
    3289       DO j = 1,ntra
    3290       DO i = 1,nl
    3291        DO il = 1,ncum
    3292         IF (iflag(il) .le. 1) THEN
    3293          ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
    3294         ENDIF
    3295        ENDDO
    3296       ENDDO
    3297       ENDDO
     3289
     3290!AC!      DO j = 1,ntra
     3291!AC!      DO i = 1,nl
     3292!AC!       DO il = 1,ncum
     3293!AC!        IF (iflag(il) .le. 1) THEN
     3294!AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
     3295!AC!        ENDIF
     3296!AC!       ENDDO
     3297!AC!      ENDDO
     3298!AC!      ENDDO
    32983299
    32993300c
     
    35393540        end
    35403541
     3542!AC!
     3543      SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
     3544     &                        ment,sij,da,phi)
     3545        implicit none
     3546c inputs:
     3547        integer ncum, nd, na, nloc,len
     3548        real ment(nloc,na,na),sij(nloc,na,na)
     3549c ouputs:
     3550        real da(nloc,na),phi(nloc,na,na)
     3551c local variables:
     3552        integer i,j,k
     3553c       
     3554        da(:,:)=0.
     3555c
     3556        do j=1,na
     3557          do k=1,na
     3558            do i=1,ncum
     3559            da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
     3560            phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
     3561            end do
     3562          end do
     3563        end do
     3564        return
     3565        end
     3566!AC!
    35413567
    35423568      SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
     
    36093635
    36103636
    3611         do 2100 j=1,ntra
    3612 c oct3         do 2110 k=1,nl
    3613          do 2110 k=1,nd ! oct3
    3614           do 2120 i=1,ncum
    3615             ftra1(idcum(i),k,j)=ftra(i,k,j)
    3616  2120     continue
    3617  2110    continue
    3618  2100   continue
     3637!AC!        do 2100 j=1,ntra
     3638!AC!c oct3         do 2110 k=1,nl
     3639!AC!         do 2110 k=1,nd ! oct3
     3640!AC!          do 2120 i=1,ncum
     3641!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
     3642!AC! 2120     continue
     3643!AC! 2110    continue
     3644!AC! 2100   continue
    36193645        return
    36203646        end
  • LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F

    r1403 r1669  
    116116 110  continue
    117117
    118       do 121 j=1,ntra
    119 ccccc      do 111 k=1,nl+1
    120       do 111 k=1,nd
    121        nn=0
    122       do 101 i=1,len
    123       if(iflag1(i).eq.0)then
    124        nn=nn+1
    125        tra(nn,k,j)=tra1(i,k,j)
    126       endif
    127  101  continue
    128  111  continue
    129  121  continue
     118!AC!      do 121 j=1,ntra
     119!AC!ccccc      do 111 k=1,nl+1
     120!AC!      do 111 k=1,nd
     121!AC!       nn=0
     122!AC!      do 101 i=1,len
     123!AC!      if(iflag1(i).eq.0)then
     124!AC!       nn=nn+1
     125!AC!       tra(nn,k,j)=tra1(i,k,j)
     126!AC!      endif
     127!AC! 101  continue
     128!AC! 111  continue
     129!AC! 121  continue
    130130
    131131      if (nn.ne.ncum) then
  • LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F

    r1518 r1669  
    99     :         ,Plim1,Plim2,asupmax,supmax0
    1010     :         ,asupmaxmin
     11!AC!
     12     :         ,da,phi
     13!AC!
    1114     o         ,iflag1,kbas1,ktop1
    1215     :         ,precip1,cbmf1,plcl1,plfc1,wbeff1,sig1,w01,ptop21
     
    1720     :         ,ftd1,fqd1
    1821     :         ,Plim11,Plim21,asupmax1,supmax01
    19      :         ,asupmaxmin1     )
     22     :         ,asupmaxmin1     
     23!AC!
     24     :         ,da1,phi1  )
     25!AC!
    2026***************************************************************
    2127*                                                             *
     
    5056      real asupmax(nloc,nd),supmax0(nloc)
    5157      real asupmaxmin(nloc)
    52 
     58!AC!
     59      real da(nloc,nd),phi(nloc,nd,nd)
     60!AC!
    5361c outputs:
    5462      integer iflag1(len),kbas1(len),ktop1(len)
     
    6876      real asupmax1(len,nd),supmax01(len)
    6977      real asupmaxmin1(len)
     78!AC!
     79      real da1(nloc,nd),phi1(nloc,nd,nd)
     80!AC!
    7081c
    7182c local variables:
     
    111122            fqd1(idcum(i),k)=fqd(i,k)
    112123            asupmax1(idcum(i),k)=asupmax(i,k)
    113  2010     continue
     124!AC!
     125            da1(idcum(i),k)=da(i,k)
     126!AC!
     127 2010    continue
    114128 2020   continue
    115129
     
    119133
    120134
    121         do 2100 j=1,ntra
    122 c oct3         do 2110 k=1,nl
    123          do 2110 k=1,nd ! oct3
    124           do 2120 i=1,ncum
    125             ftra1(idcum(i),k,j)=ftra(i,k,j)
    126  2120     continue
    127  2110    continue
    128  2100   continue
     135!AC!        do 2100 j=1,ntra
     136!AC!c oct3         do 2110 k=1,nl
     137!AC!         do 2110 k=1,nd ! oct3
     138!AC!          do 2120 i=1,ncum
     139!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
     140!AC! 2120     continue
     141!AC! 2110    continue
     142!AC! 2100   continue
     143
     144!AC!
     145       do j=1,nd
     146         do k=1,nd
     147          do i=1,ncum
     148            phi1(idcum(i),k,j)=phi(i,k,j)
     149          end do
     150         end do
     151        end do
     152!AC!
     153
    129154c
    130155c        do 2220 k2=1,nd
  • LMDZ5/branches/testing/libf/phylmd/cv3p_mixing.F

    r1664 r1669  
    118118            elij(i,k,j)=0.0
    119119            hent(i,k,j)=0.0
    120             ment(i,k,j)=0.0
    121             sij(i,k,j)=0.0
     120!AC!            ment(i,k,j)=0.0
     121!AC!            sij(i,k,j)=0.0
    122122 385      continue
    123123 390    continue
    124124 400  continue
     125
     126!AC!
     127      ment(1:ncum,1:nd,1:nd)=0.0
     128      sij(1:ncum,1:nd,1:nd)=0.0
     129!AC!
    125130
    126131      do k=1,ntra
  • LMDZ5/branches/testing/libf/phylmd/cva_driver.F

    r1518 r1669  
    2020     &                   ftd1,fqd1,
    2121     &                   Plim11,Plim21,asupmax1,supmax01,asupmaxmin1
    22      &                   ,lalim_conv)
     22     &                   ,lalim_conv,
     23!AC!
     24     &                   da1,phi1)
     25!AC!
    2326***************************************************************
    2427*                                                             *
     
    171174      real tvp1(len,nd)
    172175c
     176!AC!
     177      real da1(len,nd),phi1(len,nd,nd)
     178      real da(len,nd),phi(len,nd,nd)
     179!AC!
    173180      real ftd1(len,nd)
    174181      real fqd1(len,nd)
     
    912919      endif
    913920
     921!AC!
     922!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     923! --- passive tracers
     924!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     925
     926      if (iflag_con.eq.3) then
     927       CALL cv3_tracer(nloc,len,ncum,nd,nd,
     928     :                  ment,sij,da,phi)
     929      endif
     930
     931!AC!
     932
    914933!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    915934! --- UNCOMPRESS THE FIELDS
     
    928947     :          ,Plim1,Plim2,asupmax,supmax0
    929948     :          ,asupmaxmin
     949!AC!
     950     :          ,da,phi
     951!AC!
    930952     o          ,iflag1,kbas1,ktop1
    931953     o          ,precip1,cbmf1,plcl1,plfc1,wbeff1,sig1,w01,ptop21
     
    936958     o          ,ftd1,fqd1
    937959     o          ,Plim11,Plim21,asupmax1,supmax01
    938      o          ,asupmaxmin1)
     960     o          ,asupmaxmin1
     961!AC!
     962     o          ,da1,phi1)
     963!AC!
    939964      endif
    940965
  • LMDZ5/branches/testing/libf/phylmd/inistats.F

    r1492 r1669  
    1818      real, dimension(istime) :: lt
    1919      integer :: nvarid
    20       real, dimension(llm) :: pseudoalt
    2120
    2221      write (*,*)
  • LMDZ5/branches/testing/libf/phylmd/phys_output_mod.F90

    r1665 r1669  
    3333  CHARACTER(len=20), dimension(nfiles), private, save   :: type_ecri
    3434  !$OMP THREADPRIVATE(nhorim, nvertm, zoutm,zdtime,type_ecri)
     35 ! swaero_diag : flag indicates if it is necessary to do calculation for some aerosol diagnostics
     36  logical, save                                :: swaero_diag=.FALSE.
     37
    3538
    3639  !   integer, save                     :: nid_hf3d
     
    260263  type(ctrl_out),save :: o_wape         = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'wape')
    261264
     265!!! nrlmd le 10/04/2012
     266
     267!-------Spectre de thermiques de type 2 au LCL
     268  type(ctrl_out),save :: o_n2                = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'n2')
     269  type(ctrl_out),save :: o_s2                = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'s2')
     270                                                                             
     271!-------Déclenchement stochastique                                           
     272  type(ctrl_out),save :: o_proba_notrig      = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'proba_notrig')
     273  type(ctrl_out),save :: o_random_notrig     = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'random_notrig')
     274  type(ctrl_out),save :: o_ale_bl_stat       = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'ale_bl_stat')
     275  type(ctrl_out),save :: o_ale_bl_trig       = ctrl_out((/ 1, 1, 1, 6, 10, 10 /),'ale_bl_trig')
     276
     277!-------Fermeture statistique
     278  type(ctrl_out),save :: o_alp_bl_det        = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_det')
     279  type(ctrl_out),save :: o_alp_bl_fluct_m    = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_fluct_m')
     280  type(ctrl_out),save :: o_alp_bl_fluct_tke  = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_fluct_tke')
     281  type(ctrl_out),save :: o_alp_bl_conv       = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_conv')
     282  type(ctrl_out),save :: o_alp_bl_stat       = ctrl_out((/ 1, 1, 1, 10, 10, 10 /),'alp_bl_stat')
     283
     284!!! fin nrlmd le 10/04/2012
    262285
    263286  ! Champs interpolles sur des niveaux de pression ??? a faire correctement
     
    365388
    366389  type(ctrl_out),save :: o_topswad      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'topswad')
     390  type(ctrl_out),save :: o_topswad0     = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'topswad0')
    367391  type(ctrl_out),save :: o_topswai      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'topswai')
    368392  type(ctrl_out),save :: o_solswad      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'solswad')
     393  type(ctrl_out),save :: o_solswad0     = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'solswad0')
    369394  type(ctrl_out),save :: o_solswai      = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'solswai')
    370395
     
    432457  type(ctrl_out),save :: o_ovap         = ctrl_out((/ 2, 3, 4, 10, 10, 10 /),'ovap')
    433458  type(ctrl_out),save :: o_ovapinit     = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'ovapinit')
     459  type(ctrl_out),save :: o_oliq         = ctrl_out((/ 2, 3, 4, 10, 10, 10 /),'oliq')
    434460  type(ctrl_out),save :: o_wvapp        = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'wvapp')
    435461  type(ctrl_out),save :: o_geop         = ctrl_out((/ 2, 3, 10, 10, 10, 10 /),'geop')
     
    494520  type(ctrl_out),save :: o_dtcon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtcon')
    495521  type(ctrl_out),save :: o_ducon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'ducon')
     522  type(ctrl_out),save :: o_dvcon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dvcon')
    496523  type(ctrl_out),save :: o_dqcon        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqcon')
    497524  type(ctrl_out),save :: o_dtwak        = ctrl_out((/ 4, 5, 10, 10, 10, 10 /),'dtwak')
     
    531558  type(ctrl_out),save :: o_e_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'e_th')
    532559  type(ctrl_out),save :: o_w_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'w_th')
    533   type(ctrl_out),save :: o_lambda_th    = ctrl_out((/ 10, 10, 10, 10, 10, 10 /),'lambda_th')
    534560  type(ctrl_out),save :: o_ftime_th     = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'ftime_th')
    535561  type(ctrl_out),save :: o_q_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'q_th')
     
    537563  type(ctrl_out),save :: o_d_th         = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'d_th')
    538564  type(ctrl_out),save :: o_f0_th        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'f0_th')
    539   type(ctrl_out),save :: o_zmax_th      = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'zmax_th')
     565  type(ctrl_out),save :: o_zmax_th      = ctrl_out((/ 4,  4,  4,  5, 10, 10 /),'zmax_th')
    540566  type(ctrl_out),save :: o_dqthe        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqthe')
    541567  type(ctrl_out),save :: o_dtajs        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtajs')
     
    621647    USE infotrac
    622648    USE ioipsl
     649!    USE phys_cal_mod, only : hour
    623650    USE mod_phys_lmdz_para
    624651    USE aero_mod, only : naero_spc,name_aero
     
    682709    !                 entre [phys_out_lonmin,phys_out_lonmax] et [phys_out_latmin,phys_out_latmax]
    683710
    684     logical, dimension(nfiles), save  :: phys_out_regfkey       = (/ .false., .false., .false., .false., .false., .false. /)
    685     real, dimension(nfiles), save     :: phys_out_lonmin        = (/ -180., -180., -180., -180., -180., -180. /)
    686     real, dimension(nfiles), save     :: phys_out_lonmax        = (/ 180., 180., 180., 180., 180., 180. /)
    687     real, dimension(nfiles), save     :: phys_out_latmin        = (/ -90., -90., -90., -90., -90., -90. /)
    688     real, dimension(nfiles), save     :: phys_out_latmax        = (/ 90., 90., 90., 90., 90., 90. /)
     711    logical, dimension(nfiles), save  :: phys_out_regfkey       = (/ .false., .false., .false.,  .false., .false., .false. /)
     712    real, dimension(nfiles), save     :: phys_out_lonmin        = (/   -180.,   -180.,   -180.,    -180.,   -180.,  -180. /)
     713    real, dimension(nfiles), save     :: phys_out_lonmax        = (/    180.,    180.,    180.,     180.,    180.,    180. /)
     714    real, dimension(nfiles), save     :: phys_out_latmin        = (/    -90.,    -90.,    -90.,     -90.,    -90.,    -90. /)
     715    real, dimension(nfiles), save     :: phys_out_latmax        = (/     90.,     90.,     90.,     90.,     90.,    90. /)
    689716
    690717    write(lunout,*) 'Debut phys_output_mod.F90'
     
    792819    DO iff=1,nfiles
    793820
     821       ! Calculate ecrit_files for all files
     822       if ( chtimestep(iff).eq.'DefFreq' ) then
     823          ! Par defaut ecrit_files = (ecrit_mensuel ecrit_jour ecrit_hf ...)*86400.
     824          ecrit_files(iff)=ecrit_files(iff)*86400.
     825       else
     826          call convers_timesteps(chtimestep(iff),dtime,ecrit_files(iff))
     827       endif
     828       write(lunout,*)'ecrit_files(',iff,')= ',ecrit_files(iff)
     829
     830       zoutm(iff) = ecrit_files(iff) ! Frequence ou l on ecrit en seconde
     831
    794832       IF (clef_files(iff)) THEN
    795833
    796           if ( chtimestep(iff).eq.'DefFreq' ) then
    797              ! Par defaut ecrit_files = (ecrit_mensuel ecrit_jour ecrit_hf ...)*86400.
    798              ecrit_files(iff)=ecrit_files(iff)*86400.
    799           else
    800              call convers_timesteps(chtimestep(iff),dtime,ecrit_files(iff))
    801           endif
    802           write(lunout,*)'ecrit_files(',iff,')= ',ecrit_files(iff)
    803 
    804           zoutm(iff) = ecrit_files(iff) ! Frequence ou l on ecrit en seconde
    805 
    806834          idayref = day_ref
    807           CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
     835          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)       
     836! correction pour l heure initiale                               !jyg
     837!                                                                !jyg
     838!      CALL ymds2ju(annee_ref, 1, idayref, hour, zjulian)         !jyg
    808839
    809840!!!!!!!!!!!!!!!!! Traitement dans le cas ou l'on veut stocker sur un domaine limite !!
     
    10991130                  o_topswad%flag,o_topswad%name, "ADE at TOA", "W/m2")
    11001131             CALL histdef2d(iff,clef_stations(iff), &
     1132                  o_topswad0%flag,o_topswad0%name, "ADE clear-sky at TOA", "W/m2")
     1133             CALL histdef2d(iff,clef_stations(iff), &
    11011134                  o_solswad%flag,o_solswad%name, "ADE at SRF", "W/m2")
     1135             CALL histdef2d(iff,clef_stations(iff), &
     1136                  o_solswad0%flag,o_solswad0%name, "ADE clear-sky at SRF", "W/m2")
    11021137
    11031138             CALL histdef2d(iff,clef_stations(iff), &
     
    12161251                     o_wbeff%flag,o_wbeff%name, "Conv. updraft velocity at LFC (<100)", "m/s")
    12171252             end if
    1218              CALL histdef2d(iff,clef_stations(iff), &
    1219                   o_prw%flag,o_prw%name, "Precipitable water", "kg/m2")
    12201253             IF (.NOT.clef_stations(iff)) THEN
    12211254                !
     
    12531286          ENDIF !iflag_con .GE. 3
    12541287
     1288          CALL histdef2d(iff,clef_stations(iff), &
     1289               o_prw%flag,o_prw%name, "Precipitable water", "kg/m2")
    12551290          CALL histdef2d(iff,clef_stations(iff), &
    12561291               o_s_pblh%flag,o_s_pblh%name, "Boundary Layer Height", "m")
     
    13181353          ! Couplage conv-CL
    13191354          IF (iflag_con.GE.3) THEN
     1355             IF (iflag_coupl>=1) THEN
    13201356                CALL histdef2d(iff,clef_stations(iff), &
    13211357                     o_ale_bl%flag,o_ale_bl%name, "ALE BL", "m2/s2")
    13221358                CALL histdef2d(iff,clef_stations(iff), &
    13231359                     o_alp_bl%flag,o_alp_bl%name, "ALP BL", "m2/s2")
     1360             ENDIF
    13241361          ENDIF !(iflag_con.GE.3)
    13251362
     
    13751412          CALL histdef3d(iff,clef_stations(iff),o_theta%flag,o_theta%name, "Potential air temperature", "K" )
    13761413          CALL histdef3d(iff,clef_stations(iff),o_ovap%flag,o_ovap%name, "Specific humidity", "kg/kg" )
     1414          CALL histdef3d(iff,clef_stations(iff),o_oliq%flag,o_oliq%name, "Condensed water", "kg/kg" )
    13771415          CALL histdef3d(iff,clef_stations(iff), &
    13781416               o_ovapinit%flag,o_ovapinit%name, "Specific humidity (begin of timestep)", "kg/kg" )
     
    14801518               o_ducon%flag,o_ducon%name, "Convection du", "m/s2")
    14811519          CALL histdef3d(iff,clef_stations(iff), &
     1520               o_dvcon%flag,o_dvcon%name, "Convection dv", "m/s2")
     1521          CALL histdef3d(iff,clef_stations(iff), &
    14821522               o_dqcon%flag,o_dqcon%name, "Convection dQ", "(kg/kg)/s")
    14831523
     
    14891529                CALL histdef2d(iff,clef_stations(iff), &
    14901530                     o_alp_wk%flag,o_alp_wk%name, "ALP WK", "m2/s2")
     1531                CALL histdef2d(iff,clef_stations(iff), &
     1532                     o_ale%flag,o_ale%name, "ALE", "m2/s2")
     1533                CALL histdef2d(iff,clef_stations(iff), &
     1534                     o_alp%flag,o_alp%name, "ALP", "W/m2")
     1535                CALL histdef2d(iff,clef_stations(iff),o_cin%flag,o_cin%name, "Convective INhibition", "m2/s2")
     1536                CALL histdef2d(iff,clef_stations(iff),o_wape%flag,o_WAPE%name, "WAPE", "m2/s2")
    14911537                CALL histdef2d(iff,clef_stations(iff),o_wake_h%flag,o_wake_h%name, "wake_h", "-")
    14921538                CALL histdef2d(iff,clef_stations(iff),o_wake_s%flag,o_wake_s%name, "wake_s", "-")
     
    14961542                CALL histdef3d(iff,clef_stations(iff),o_wake_deltaq%flag,o_wake_deltaq%name, "wake_deltaq", " ")
    14971543                CALL histdef3d(iff,clef_stations(iff),o_wake_omg%flag,o_wake_omg%name, "wake_omg", " ")
    1498                 CALL histdef2d(iff,clef_stations(iff),o_wape%flag,o_WAPE%name, "WAPE", "m2/s2")
    14991544             ENDIF
    1500              CALL histdef2d(iff,clef_stations(iff), &
    1501                      o_ale%flag,o_ale%name, "ALE", "m2/s2")
    1502              CALL histdef2d(iff,clef_stations(iff), &
    1503                      o_alp%flag,o_alp%name, "ALP", "W/m2")
    1504              CALL histdef2d(iff,clef_stations(iff),o_cin%flag,o_cin%name, "Convective INhibition", "m2/s2")
    15051545             CALL histdef3d(iff,clef_stations(iff),o_Vprecip%flag,o_Vprecip%name, "precipitation vertical profile", "-")
    15061546             CALL histdef3d(iff,clef_stations(iff),o_ftd%flag,o_ftd%name, "tend temp due aux descentes precip", "-")
    15071547             CALL histdef3d(iff,clef_stations(iff),o_fqd%flag,o_fqd%name,"tend vap eau due aux descentes precip", "-")
    15081548          ENDIF !(iflag_con.EQ.3)
     1549
     1550!!! nrlmd le 10/04/2012
     1551
     1552        IF (iflag_trig_bl>=1) THEN
     1553 CALL histdef2d(iff,clef_stations(iff),o_n2%flag,o_n2%name, "Nombre de panaches de type 2", " ")
     1554 CALL histdef2d(iff,clef_stations(iff),o_s2%flag,o_s2%name, "Surface moyenne des panaches de type 2", "m2")
     1555
     1556 CALL histdef2d(iff,clef_stations(iff),o_proba_notrig%flag,o_proba_notrig%name, "Probabilité de non-déclenchement", " ")
     1557 CALL histdef2d(iff,clef_stations(iff),o_random_notrig%flag,o_random_notrig%name, "Tirage aléatoire de non-déclenchement", " ")
     1558 CALL histdef2d(iff,clef_stations(iff),o_ale_bl_trig%flag,o_ale_bl_trig%name, "ALE_BL_STAT + Condition P>Pseuil", "m2/s2")
     1559 CALL histdef2d(iff,clef_stations(iff),o_ale_bl_stat%flag,o_ale_bl_stat%name, "ALE_BL_STAT", "m2/s2")
     1560       ENDIF  !(iflag_trig_bl>=1)
     1561
     1562        IF (iflag_clos_bl>=1) THEN
     1563 CALL histdef2d(iff,clef_stations(iff),o_alp_bl_det%flag,o_alp_bl_det%name, "ALP_BL_DET", "W/m2")
     1564 CALL histdef2d(iff,clef_stations(iff),o_alp_bl_fluct_m%flag,o_alp_bl_fluct_m%name, "ALP_BL_FLUCT_M", "W/m2")
     1565 CALL histdef2d(iff,clef_stations(iff),o_alp_bl_fluct_tke%flag,o_alp_bl_fluct_tke%name, "ALP_BL_FLUCT_TKE", "W/m2")
     1566 CALL histdef2d(iff,clef_stations(iff),o_alp_bl_conv%flag,o_alp_bl_conv%name, "ALP_BL_CONV", "W/m2")
     1567 CALL histdef2d(iff,clef_stations(iff),o_alp_bl_stat%flag,o_alp_bl_stat%name, "ALP_BL_STAT", "W/m2")
     1568       ENDIF  !(iflag_clos_bl>=1)
     1569
     1570!!! fin nrlmd le 10/04/2012
    15091571
    15101572          CALL histdef3d(iff,clef_stations(iff),o_dtlsc%flag,o_dtlsc%name, "Condensation dT", "K/s")
     
    15191581          CALL histdef3d(iff,clef_stations(iff),o_dtthe%flag,o_dtthe%name, "Thermal dT", "K/s")
    15201582
    1521           if(iflag_thermals.gt.1) THEN
     1583          if(iflag_thermals.ge.1) THEN
    15221584             CALL histdef3d(iff,clef_stations(iff),o_dqlscth%flag,o_dqlscth%name, "dQ therm.", "(kg/kg)/s")
    15231585             CALL histdef3d(iff,clef_stations(iff),o_dqlscst%flag,o_dqlscst%name, "dQ strat.", "(kg/kg)/s")
     
    15311593             CALL histdef3d(iff,clef_stations(iff),o_e_th%flag,o_e_th%name,"Thermal plume entrainment","K/s")
    15321594             CALL histdef3d(iff,clef_stations(iff),o_w_th%flag,o_w_th%name,"Thermal plume vertical velocity","m/s")
    1533              CALL histdef3d(iff,clef_stations(iff), &
    1534                   o_lambda_th%flag,o_lambda_th%name,"Thermal plume vertical velocity","m/s")
    15351595             CALL histdef2d(iff,clef_stations(iff), &
    15361596                  o_ftime_th%flag,o_ftime_th%name,"Fraction of time Shallow convection occurs"," ")
     
    15481608             CALL histdef3d(iff,clef_stations(iff), &
    15491609                  o_dqthe%flag,o_dqthe%name, "Thermal dQ", "(kg/kg)/s")
    1550           endif !iflag_thermals.gt.1
     1610          endif !iflag_thermals.ge.1
    15511611          CALL histdef3d(iff,clef_stations(iff), &
    15521612               o_dtajs%flag,o_dtajs%name, "Dry adjust. dT", "K/s")
     
    17081768
    17091769    ENDDO !  iff
     1770
     1771    ! Updated write frequencies due to phys_out_filetimesteps.
     1772    ! Write frequencies are now in seconds. 
     1773    ecrit_mth = ecrit_files(1)
     1774    ecrit_day = ecrit_files(2)
     1775    ecrit_hf  = ecrit_files(3)
     1776    ecrit_ins = ecrit_files(4)
     1777    ecrit_LES = ecrit_files(5)
     1778    ecrit_ins = ecrit_files(6)
     1779
     1780    write(lunout,*)'swaero_diag=',swaero_diag
    17101781    write(lunout,*)'Fin phys_output_mod.F90'
    17111782  end subroutine phys_output_open
     
    17561827       endif
    17571828    endif
     1829
     1830    ! Set swaero_diag=true if at least one of the concerned variables are defined
     1831    if (nomvar=='topswad' .OR. nomvar=='topswai' .OR. nomvar=='solswad' .OR. nomvar=='solswai' ) THEN
     1832       if  ( flag_var(iff)<=lev_files(iff) ) then
     1833          swaero_diag=.TRUE.
     1834       end if
     1835    end if
    17581836  end subroutine histdef2d
    17591837
  • LMDZ5/branches/testing/libf/phylmd/phys_output_write.h

    r1665 r1669  
    508508        ENDIF
    509509
    510        if (iflag_pbl>1 .and. lev_histday.gt.10 ) then
     510       if (iflag_pbl>1 .and. lev_files(iff).gt.10 ) then
    511511        IF (o_tke_srf(nsrf)%flag(iff)<=lev_files(iff)) THEN
    512512        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     
    631631        end if
    632632
    633         IF (o_prw%flag(iff)<=lev_files(iff)) THEN
    634       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    635      $o_prw%name,itau_w,prw)
    636         ENDIF
    637 
    638633      IF (.NOT.clef_stations(iff)) THEN
    639634      IF (o_cape_max%flag(iff)<=lev_files(iff)) THEN
     
    671666
    672667       IF (o_mc%flag(iff)<=lev_files(iff)) THEN
    673         if(iflag_thermals.gt.1)then
     668        if(iflag_thermals>=1)then
    674669         zx_tmp_fi3d=dnwd+dnwd0+upwd+fm_therm(:,1:klev)
    675670        else
     
    681676     
    682677      ENDIF !iflag_con .GE. 3
     678
     679        IF (o_prw%flag(iff)<=lev_files(iff)) THEN
     680      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     681     $o_prw%name,itau_w,prw)
     682        ENDIF
    683683
    684684        IF (o_s_pblh%flag(iff)<=lev_files(iff)) THEN
     
    801801! Couplage convection-couche limite
    802802      IF (iflag_con.GE.3) THEN
     803      IF (iflag_coupl>=1) THEN
    803804       IF (o_ale_bl%flag(iff)<=lev_files(iff)) THEN
    804805       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     
    809810     $o_alp_bl%name,itau_w,alp_bl)
    810811       ENDIF
     812      ENDIF !iflag_coupl>=1
    811813      ENDIF !(iflag_con.GE.3)
    812814
     
    823825       ENDIF
    824826
     827       IF (o_ale%flag(iff)<=lev_files(iff)) THEN
     828       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     829     $o_ale%name,itau_w,ale)
     830       ENDIF
     831       IF (o_alp%flag(iff)<=lev_files(iff)) THEN
     832       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     833     $o_alp%name,itau_w,alp)
     834       ENDIF
     835       IF (o_cin%flag(iff)<=lev_files(iff)) THEN
     836       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     837     $o_cin%name,itau_w,cin)
     838       ENDIF
    825839       IF (o_wape%flag(iff)<=lev_files(iff)) THEN
    826840       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     
    869883      ENDIF ! iflag_wake>=1
    870884
    871        IF (o_ale%flag(iff)<=lev_files(iff)) THEN
    872        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    873      $o_ale%name,itau_w,ale)
    874        ENDIF
    875        IF (o_alp%flag(iff)<=lev_files(iff)) THEN
    876        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    877      $o_alp%name,itau_w,alp)
    878        ENDIF
    879        IF (o_cin%flag(iff)<=lev_files(iff)) THEN
    880        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    881      $o_cin%name,itau_w,cin)
    882        ENDIF
    883885        IF (o_Vprecip%flag(iff)<=lev_files(iff)) THEN
    884886       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     
    897899      ENDIF !(iflag_con.EQ.3)
    898900 
     901!!! nrlmd le 10/04/2012
     902
     903        IF (iflag_trig_bl>=1) THEN
     904          IF (o_n2%flag(iff)<=lev_files(iff)) THEN
     905        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     906     s                     o_n2%name,itau_w,n2)
     907         ENDIF
     908
     909         IF (o_s2%flag(iff)<=lev_files(iff)) THEN
     910        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     911     s                     o_s2%name,itau_w,s2)
     912         ENDIF
     913
     914          IF (o_proba_notrig%flag(iff)<=lev_files(iff)) THEN
     915        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     916     s                     o_proba_notrig%name,itau_w,proba_notrig)
     917         ENDIF
     918
     919         IF (o_random_notrig%flag(iff)<=lev_files(iff)) THEN
     920        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     921     s                     o_random_notrig%name,itau_w,random_notrig)
     922         ENDIF
     923
     924         IF (o_ale_bl_stat%flag(iff)<=lev_files(iff)) THEN
     925        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     926     s                     o_ale_bl_stat%name,itau_w,ale_bl_stat)
     927         ENDIF
     928
     929         IF (o_ale_bl_trig%flag(iff)<=lev_files(iff)) THEN
     930        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     931     s                     o_ale_bl_trig%name,itau_w,ale_bl_trig)
     932         ENDIF
     933       ENDIF  !(iflag_trig_bl>=1)
     934
     935        IF (iflag_clos_bl>=1) THEN
     936         IF (o_alp_bl_det%flag(iff)<=lev_files(iff)) THEN
     937        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     938     s                     o_alp_bl_det%name,itau_w,alp_bl_det)
     939         ENDIF
     940
     941         IF (o_alp_bl_fluct_m%flag(iff)<=lev_files(iff)) THEN
     942        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     943     s                     o_alp_bl_fluct_m%name,itau_w,alp_bl_fluct_m)
     944         ENDIF
     945
     946         IF (o_alp_bl_fluct_tke%flag(iff)<=lev_files(iff)) THEN
     947        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     948     s                o_alp_bl_fluct_tke%name,itau_w,alp_bl_fluct_tke)
     949         ENDIF
     950
     951         IF (o_alp_bl_conv%flag(iff)<=lev_files(iff)) THEN
     952        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     953     s                     o_alp_bl_conv%name,itau_w,alp_bl_conv)
     954         ENDIF
     955
     956         IF (o_alp_bl_stat%flag(iff)<=lev_files(iff)) THEN
     957        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     958     s                     o_alp_bl_stat%name,itau_w,alp_bl_stat)
     959         ENDIF
     960       ENDIF  !(iflag_clos_bl>=1)
     961
     962!!! fin nrlmd le 10/04/2012
     963
    899964      IF (type_ocean=='slab ') THEN
    900965      IF ( o_slab_bils%flag(iff)<=lev_files(iff))
     
    11931258     $            topswad_aero)
    11941259          ENDIF
     1260          IF (o_topswad0%flag(iff)<=lev_files(iff)) THEN
     1261             CALL histwrite_phy(nid_files(iff),
     1262     $clef_stations(iff),
     1263     $o_topswad0%name,itau_w,
     1264     $            topswad0_aero)
     1265          ENDIF
    11951266          IF (o_solswad%flag(iff)<=lev_files(iff)) THEN
    11961267             CALL histwrite_phy(nid_files(iff),
     
    11981269     $o_solswad%name,itau_w,
    11991270     $            solswad_aero)
     1271          ENDIF
     1272          IF (o_solswad0%flag(iff)<=lev_files(iff)) THEN
     1273             CALL histwrite_phy(nid_files(iff),
     1274     $clef_stations(iff),
     1275     $o_solswad0%name,itau_w,
     1276     $            solswad0_aero)
    12001277          ENDIF
    12011278
     
    14101487      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    14111488     $                   o_ovap%name,itau_w,q_seri)
     1489       ENDIF
     1490                                                               
     1491       IF (o_oliq%flag(iff)<=lev_files(iff)) THEN             
     1492      CALL histwrite_phy(nid_files(iff),clef_stations(iff),   
     1493     $                   o_oliq%name,itau_w,ql_seri)           
    14121494       ENDIF
    14131495
     
    16371719       ENDIF
    16381720
     1721       IF (o_dvcon%flag(iff)<=lev_files(iff)) THEN
     1722      zx_tmp_fi3d(1:klon,1:klev)=d_v_con(1:klon,1:klev)/pdtphys
     1723      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     1724     $o_dvcon%name,itau_w,zx_tmp_fi3d)
     1725       ENDIF
     1726
    16391727       IF (o_dqcon%flag(iff)<=lev_files(iff)) THEN
    16401728      zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys
     
    16801768!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    16811769! Sorties specifiques a la separation thermiques/non thermiques
    1682        if (iflag_thermals>1) then
     1770       if (iflag_thermals>=1) then
    16831771
    16841772       IF (o_dtlscth%flag(iff)<=lev_files(iff)) THEN
     
    17441832       ENDIF
    17451833
    1746       endif ! iflag_thermals>1
     1834      endif ! iflag_thermals>=1
    17471835
    17481836!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    17911879       ENDIF
    17921880
    1793        IF (iflag_thermals.gt.1) THEN
     1881       IF (iflag_thermals>=1) THEN
    17941882        IF (o_ftime_th%flag(iff)<=lev_files(iff)) THEN
    17951883! Pour l instant 0 a y reflichir pour les thermiques
     
    18201908        ENDIF
    18211909
    1822         IF (o_lambda_th%flag(iff)<=lev_files(iff)) THEN
    1823         CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    1824      s                     o_lambda_th%name,itau_w,lambda_th)
    1825         ENDIF
    18261910
    18271911        IF (o_a_th%flag(iff)<=lev_files(iff)) THEN
     
    20512135       ENDIF
    20522136
    2053        IF (o_mcd%flag(iff)<=lev_files(iff)) THEN
    2054       zx_tmp_fi3d(1:klon,1:klev)=-1 * (dnwd(1:klon,1:klev)+
    2055      $                                 dnwd0(1:klon,1:klev))
    2056       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    2057      $o_mcd%name,itau_w,zx_tmp_fi3d)
    2058        ENDIF
    2059 
    2060        IF (o_dmc%flag(iff)<=lev_files(iff)) THEN
    2061       zx_tmp_fi3d(1:klon,1:klev)=upwd(1:klon,1:klev) +
    2062      $  dnwd(1:klon,1:klev)+ dnwd0(1:klon,1:klev)
    2063       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    2064      $o_dmc%name,itau_w,zx_tmp_fi3d)
    2065        ENDIF
     2137       if (iflag_con >= 3) then
     2138          IF (o_mcd%flag(iff)<=lev_files(iff)) THEN
     2139             zx_tmp_fi3d(1:klon,1:klev)=-1 * (dnwd(1:klon,1:klev)+
     2140     $            dnwd0(1:klon,1:klev))
     2141             CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     2142     $            o_mcd%name,itau_w,zx_tmp_fi3d)
     2143          ENDIF
     2144
     2145          IF (o_dmc%flag(iff)<=lev_files(iff)) THEN
     2146             zx_tmp_fi3d(1:klon,1:klev)=upwd(1:klon,1:klev) +
     2147     $            dnwd(1:klon,1:klev)+ dnwd0(1:klon,1:klev)
     2148             CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     2149     $            o_dmc%name,itau_w,zx_tmp_fi3d)
     2150          ENDIF
     2151       else if (iflag_con == 2) then
     2152          IF (o_mcd%flag(iff) <= lev_files(iff)) THEN
     2153             CALL histwrite_phy(nid_files(iff), clef_stations(iff),
     2154     $            o_mcd%name, itau_w, pmfd)
     2155          ENDIF
     2156
     2157          IF (o_dmc%flag(iff) <= lev_files(iff)) THEN
     2158             CALL histwrite_phy(nid_files(iff), clef_stations(iff),
     2159     $            o_dmc%name, itau_w, pmfu + pmfd)
     2160          ENDIF
     2161       end if
    20662162
    20672163       IF (o_ref_liq%flag(iff)<=lev_files(iff)) THEN
  • LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90

    r1539 r1669  
    346346!$OMP THREADPRIVATE(ccm)
    347347
     348!!! nrlmd le 10/04/2012
     349      REAL,SAVE,ALLOCATABLE :: ale_bl_trig(:)
     350!$OMP THREADPRIVATE(ale_bl_trig)
     351!!! fin nrlmd le 10/04/2012
     352
    348353CONTAINS
    349354
     
    496501      ALLOCATE(tau_aero(klon,klev,naero_grp,nbands),piz_aero(klon,klev,naero_grp,nbands),cg_aero(klon,klev,naero_grp,nbands))
    497502      ALLOCATE(ccm(klon,klev,nbands))
     503
     504!!! nrlmd le 10/04/2012
     505      ALLOCATE(ale_bl_trig(klon))
     506!!! fin nrlmd le 10/04/2012
    498507
    499508END SUBROUTINE phys_state_var_init
     
    603612      deallocate(ccm)
    604613       
     614!!! nrlmd le 10/04/2012
     615      deallocate(ale_bl_trig)
     616!!! fin nrlmd le 10/04/2012
     617
    605618END SUBROUTINE phys_state_var_end
    606619
  • LMDZ5/branches/testing/libf/phylmd/physiq.F

    r1665 r1669  
    180180      real facteur,zfratqs1,zfratqs2
    181181
    182       REAL lambda_th(klon,klev),zz,znum,zden
     182      REAL zz,znum,zden
    183183      REAL wmax_th(klon)
    184184      REAL zmax_th(klon)
     
    614614      REAL dd_t(klon,klev),dd_q(klon,klev)
    615615
    616       real, save :: alp_bl_prescr=0.1
    617       real, save :: ale_bl_prescr=4.
     616      real, save :: alp_bl_prescr=0.
     617      real, save :: ale_bl_prescr=0.
    618618
    619619      real, save :: ale_max=1000.
     
    689689      REAL ztla(klon,klev)
    690690      REAL zthl(klon,klev)
     691
     692ccc nrlmd le 10/04/2012
     693
     694c--------Stochastic Boundary Layer Triggering: ALE_BL--------
     695c---Propriétés du thermiques au LCL
     696      real zlcl_th(klon)                                     ! Altitude du LCL calculé continument (pcon dans thermcell_main.F90)
     697      real fraca0(klon)                                      ! Fraction des thermiques au LCL
     698      real w0(klon)                                          ! Vitesse des thermiques au LCL
     699      real w_conv(klon)                                      ! Vitesse verticale de grande échelle au LCL
     700      real therm_tke_max0(klon)                              ! TKE dans les thermiques au LCL
     701      real env_tke_max0(klon)                                ! TKE dans l'environnement au LCL
     702
     703c---Spectre de thermiques de type 2 au LCL
     704      real n2(klon),s2(klon)
     705      real ale_bl_stat(klon)
     706
     707c---Déclenchement stochastique
     708      integer :: tau_trig(klon)
     709      real proba_notrig(klon)
     710      real random_notrig(klon)
     711
     712c--------Statistical Boundary Layer Closure: ALP_BL--------
     713c---Profils de TKE dans et hors du thermique
     714      real pbl_tke_input(klon,klev+1,nbsrf)
     715      real therm_tke_max(klon,klev)                          ! Profil de TKE dans les thermiques
     716      real env_tke_max(klon,klev)                            ! Profil de TKE dans l'environnement
     717
     718c---Fermeture statistique
     719      real alp_bl_det(klon)                                     ! ALP déterministe du thermique unique
     720      real alp_bl_fluct_m(klon)                                 ! ALP liée aux fluctuations de flux de masse sous-nuageux
     721      real alp_bl_fluct_tke(klon)                               ! ALP liée aux fluctuations d'énergie cinétique sous-nuageuse
     722      real alp_bl_conv(klon)                                    ! ALP liée à grande échelle
     723      real alp_bl_stat(klon)                                    ! ALP totale
     724
     725ccc fin nrlmd le 10/04/2012
    691726
    692727c Variables locales pour la couche limite (al1):
     
    12121247      LOGICAL, SAVE :: mskocean_beta
    12131248c$OMP THREADPRIVATE(mskocean_beta)
    1214       REAL, dimension(klon, klev) :: beta       ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
    1215       REAL, dimension(klon, klev) :: cldtaurad  ! epaisseur optique pour radlwsw,COSP
    1216       REAL, dimension(klon, klev) :: cldemirad  ! emissivite pour radlwsw,COSP
     1249      REAL, dimension(klon, klev) :: beta         ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
     1250      REAL, dimension(klon, klev) :: cldtaurad    ! epaisseur optique pour radlwsw,COSP
     1251      REAL, dimension(klon, klev) :: cldtaupirad  ! epaisseur optique pour radlwsw,COSP cas pre-industrial
     1252      REAL, dimension(klon, klev) :: cldemirad    ! emissivite pour radlwsw,COSP
    12171253      INTEGER :: nbtr_tmp ! Number of tracer inside concvl
    12181254      REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
     
    13541390         solswad(:)=0.
    13551391
    1356          lambda_th(:,:)=0.
    13571392         wmax_th(:)=0.
    13581393         tau_overturning_th(:)=0.
     
    14901525cCR:04.12.07: initialisations poches froides
    14911526c Controle de ALE et ALP pour la fermeture convective (jyg)
    1492          CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
     1527          if (iflag_wake>=1) then
     1528            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
    14931529     s                  ,alp_bl_prescr, ale_bl_prescr)
    14941530c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
    14951531c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
     1532          endif
    14961533
    14971534        do i = 1,klon
     
    15161553       print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
    15171554      ENDIF
     1555
    15181556c
    15191557      ALLOCATE(tabCFMIP(nCFMIP))
     
    16241662
    16251663#endif
    1626 
    1627 
    1628          ecrit_hf = ecrit_hf * un_jour
    1629 cIM
    1630          IF(ecrit_day.LE.1.) THEN
    1631           ecrit_day = ecrit_day * un_jour !en secondes
    1632          ENDIF
    1633 cIM
    1634          ecrit_mth = ecrit_mth * un_jour
    1635          ecrit_ins = ecrit_ins * un_jour
    16361664         ecrit_reg = ecrit_reg * un_jour
    16371665         ecrit_tra = ecrit_tra * un_jour
    1638          ecrit_LES = ecrit_LES * un_jour
    1639 c
    1640 
     1666       
    16411667cXXXPB Positionner date0 pour initialisation de ORCHIDEE
    16421668      date0 = jD_ref
     
    17351761!
    17361762      itap   = itap + 1
     1763c
    17371764!
    17381765! Update fraction of the sub-surfaces (pctsrf) and
     
    20422069c
    20432070
    2044       if (iflag_pbl/=0) then 
    2045 
    2046         CALL pbl_surface(
    2047      e       dtime,     date0,     itap,    days_elapsed+1,
    2048      e       debut,     lafin,
    2049      e       rlon,      rlat,      rugoro,  rmu0,     
    2050      e       rain_fall, snow_fall, solsw,   sollw,   
    2051      e       t_seri,    q_seri,    u_seri,  v_seri,   
    2052      e       pplay,     paprs,     pctsrf,           
    2053      +       ftsol,     falb1,     falb2,   u10m,   v10m,
    2054      s       sollwdown, cdragh,    cdragm,  u1,    v1,
    2055      s       albsol1,   albsol2,   sens,    evap, 
    2056      s       zxtsol,    zxfluxlat, zt2m,    qsat2m,
    2057      s       d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
    2058      s       coefh,     coefm,     slab_wfbils,               
    2059      d       qsol,      zq2m,      s_pblh,  s_lcl,
    2060      d       s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
    2061      d       s_therm,   s_trmb1,   s_trmb2, s_trmb3,
    2062      d       zxrugs,    zu10m,     zv10m,   fder,
    2063      d       zxqsurf,   rh2m,      zxfluxu, zxfluxv,
    2064      d       frugs,     agesno,    fsollw,  fsolsw,
    2065      d       d_ts,      fevap,     fluxlat, t2m,
    2066      d       wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
    2067      -       dsens,     devap,     zxsnow,
    2068      -       zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
     2071      if (iflag_pbl/=0) then
     2072
     2073      CALL pbl_surface(
     2074     e     dtime,     date0,     itap,    days_elapsed+1,
     2075     e     debut,     lafin,
     2076     e     rlon,      rlat,      rugoro,  rmu0,     
     2077     e     rain_fall, snow_fall, solsw,   sollw,   
     2078     e     t_seri,    q_seri,    u_seri,  v_seri,   
     2079     e     pplay,     paprs,     pctsrf,           
     2080     +     ftsol,     falb1,     falb2,   u10m,   v10m,
     2081     s     sollwdown, cdragh,    cdragm,  u1,    v1,
     2082     s     albsol1,   albsol2,   sens,    evap, 
     2083     s     zxtsol,    zxfluxlat, zt2m,    qsat2m,
     2084     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
     2085     s     coefh,     coefm,     slab_wfbils,               
     2086     d     qsol,      zq2m,      s_pblh,  s_lcl,
     2087     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
     2088     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
     2089     d     zxrugs,    zu10m,     zv10m,   fder,
     2090     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
     2091     d     frugs,     agesno,    fsollw,  fsolsw,
     2092     d     d_ts,      fevap,     fluxlat, t2m,
     2093     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
     2094     -     dsens,     devap,     zxsnow,
     2095     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
    20692096
    20702097
    20712098!-----------------------------------------------------------------------------------------
    20722099! ajout des tendances de la diffusion turbulente
    2073         CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
     2100      CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
    20742101!-----------------------------------------------------------------------------------------
    20752102
    2076         if (mydebug) then
    2077           call writefield_phy('u_seri',u_seri,llm)
    2078           call writefield_phy('v_seri',v_seri,llm)
    2079           call writefield_phy('t_seri',t_seri,llm)
    2080           call writefield_phy('q_seri',q_seri,llm)
    2081         endif
    2082 
    2083 
    2084         IF (ip_ebil_phy.ge.2) THEN
    2085           ztit='after surface_main'
    2086           CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
     2103      if (mydebug) then
     2104        call writefield_phy('u_seri',u_seri,llm)
     2105        call writefield_phy('v_seri',v_seri,llm)
     2106        call writefield_phy('t_seri',t_seri,llm)
     2107        call writefield_phy('q_seri',q_seri,llm)
     2108      endif
     2109
     2110
     2111      IF (ip_ebil_phy.ge.2) THEN
     2112        ztit='after surface_main'
     2113        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
    20872114     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
    20882115     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2089           call diagphy(airephy,ztit,ip_ebil_phy
     2116         call diagphy(airephy,ztit,ip_ebil_phy
    20902117     e      , zero_v, zero_v, zero_v, zero_v, sens
    20912118     e      , evap  , zero_v, zero_v, ztsol
    20922119     e      , d_h_vcol, d_qt, d_ec
    20932120     s      , fs_bound, fq_bound )
    2094         END IF
     2121      END IF
    20952122
    20962123      ENDIF
    2097 
    20982124c =================================================================== c
    20992125c   Calcul de Qsat
     
    22442270cdans le thermique sinon
    22452271       if (iflag_coupl.eq.0) then
    2246           if (debut.and.prt_level.gt.9)WRITE(lunout,*) 'ALE&ALP imposes'
    2247           Ale_bl(1:klon) = ale_bl_prescr
    2248           Alp_bl(1:klon) = alp_bl_prescr
     2272          if (debut.and.prt_level.gt.9)
     2273     $                     WRITE(lunout,*)'ALE et ALP imposes'
     2274          do i = 1,klon
     2275con ne couple que ale
     2276c           ALE(i) = max(ale_wake(i),Ale_bl(i))
     2277            ALE(i) = max(ale_wake(i),ale_bl_prescr)
     2278con ne couple que alp
     2279c           ALP(i) = alp_wake(i) + Alp_bl(i)
     2280            ALP(i) = alp_wake(i) + alp_bl_prescr
     2281          enddo
    22492282       else
    22502283         IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
    2251        endif
     2284!         do i = 1,klon
     2285!             ALE(i) = max(ale_wake(i),Ale_bl(i))
     2286! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
     2287!             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
     2288!         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
     2289!         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
     2290!         enddo
    22522291
    22532292!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    22562295! w si <0
    22572296!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2258 
    22592297       do i = 1,klon
    22602298          ALE(i) = max(ale_wake(i),Ale_bl(i))
     2299ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
     2300          if (iflag_trig_bl.ge.1) then
     2301             ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
     2302          endif
     2303ccc fin nrlmd le 10/04/2012
    22612304          if (alp_offset>=0.) then
    22622305            ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
     
    22692312          endif
    22702313       enddo
    2271 
    22722314!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    22732315
     2316       endif
    22742317       do i=1,klon
    22752318          if (alp(i)>alp_max) then
     
    25862629
    25872630
    2588          if (iflag_thermals.gt.1) then
     2631ccc nrlmd le 10/04/2012
     2632         DO k=1,klev+1
     2633           DO i=1,klon
     2634           pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
     2635           pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
     2636           pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
     2637           pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
     2638           ENDDO
     2639         ENDDO
     2640ccc fin nrlmd le 10/04/2012
     2641
     2642         if (iflag_thermals>=1) then
    25892643         call calltherm(pdtphys
    25902644     s      ,pplay,paprs,pphi,weak_inversion
     
    25962650con rajoute ale et alp, et les caracteristiques de la couche alim
    25972651     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
    2598      s      ,ztv,zpspsk,ztla,zthl)
     2652     s      ,ztv,zpspsk,ztla,zthl
     2653ccc nrlmd le 10/04/2012
     2654     e      ,pbl_tke_input,pctsrf,omega,airephy
     2655     s      ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0
     2656     s      ,n2,s2,ale_bl_stat
     2657     s      ,therm_tke_max,env_tke_max
     2658     s      ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke
     2659     s      ,alp_bl_conv,alp_bl_stat
     2660ccc fin nrlmd le 10/04/2012
     2661     s                 )
     2662
     2663ccc nrlmd le 10/04/2012
     2664c-----------Stochastic triggering-----------
     2665      if (iflag_trig_bl.ge.1) then
     2666c
     2667        IF (prt_level .GE. 10) THEN
     2668         print *,'cin, ale_bl_stat, alp_bl_stat ',
     2669     $            cin, ale_bl_stat, alp_bl_stat
     2670        ENDIF
     2671
     2672c----Initialisations
     2673      do i=1,klon
     2674      proba_notrig(i)=1.
     2675      random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
     2676        if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
     2677        tau_trig(i)=tau_trig_shallow
     2678        else
     2679        tau_trig(i)=tau_trig_deep
     2680        endif
     2681      enddo
     2682c
     2683        IF (prt_level .GE. 10) THEN
     2684         print *,'random_notrig, tau_trig ',
     2685     $            random_notrig, tau_trig
     2686          print *,'s_trig,s2,n2 ',
     2687     $             s_trig,s2,n2
     2688        ENDIF
     2689
     2690c----Tirage aléatoire et calcul de ale_bl_trig
     2691      do i=1,klon
     2692        if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
     2693        proba_notrig(i)=(1.-exp(-s_trig/s2(i)))**
     2694     $                  (n2(i)*dtime/tau_trig(i))
     2695c        print *, 'proba_notrig(i) ',proba_notrig(i)
     2696          if (random_notrig(i) .ge. proba_notrig(i)) then
     2697          ale_bl_trig(i)=ale_bl_stat(i)
     2698          else
     2699          ale_bl_trig(i)=0.
     2700          endif
     2701        else
     2702        proba_notrig(i)=1.
     2703        random_notrig(i)=0.
     2704        ale_bl_trig(i)=0.
     2705        endif
     2706      enddo
     2707c
     2708        IF (prt_level .GE. 10) THEN
     2709         print *,'proba_notrig, ale_bl_trig ',
     2710     $            proba_notrig, ale_bl_trig
     2711        ENDIF
     2712
     2713      endif !(iflag_trig_bl)
     2714
     2715c-----------Statistical closure-----------
     2716      if (iflag_clos_bl.ge.1) then
     2717
     2718        do i=1,klon
     2719        alp_bl(i)=alp_bl_stat(i)
     2720        enddo
     2721
     2722      else
     2723
     2724      alp_bl_stat(:)=0.
     2725
     2726      endif !(iflag_clos_bl)
     2727
     2728        IF (prt_level .GE. 10) THEN
     2729         print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
     2730        ENDIF
     2731
     2732ccc fin nrlmd le 10/04/2012
    25992733
    26002734! ----------------------------------------------------------------------
     
    26272761c  ==============
    26282762
    2629 ! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
     2763! Dans le cas où on active les thermiques, on fait partir l'ajustement
    26302764! a partir du sommet des thermiques.
    26312765! Dans le cas contraire, on demarre au niveau 1.
     
    28142948! FH 22/09/2009
    28152949! La ligne ci-dessous faisait osciller le modele et donnait une solution
    2816 ! asymptotique bidon et d\'ependant fortement du pas de temps.
     2950! assymptotique bidon et dépendant fortement du pas de temps.
    28172951!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
    28182952!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    28422976c Appeler le processus de condensation a grande echelle
    28432977c et le processus de precipitation
     2978c-------------------------------------------------------------------------
     2979      IF (prt_level .GE.10) THEN
     2980       print *,' ->fisrtilp '
     2981      ENDIF
    28442982c-------------------------------------------------------------------------
    28452983      CALL fisrtilp(dtime,paprs,pplay,
     
    29623100cjq - introduce the aerosol direct and first indirect radiative forcings
    29633101cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
    2964       IF (ok_ade.OR.ok_aie) THEN
     3102      IF (flag_aerosol .gt. 0) THEN
    29653103         IF (.NOT. aerosol_couple)
    29663104     &        CALL readaerosol_optic(
     
    32473385cIM betaCRF
    32483386c
    3249       cldtaurad = cldtau
    3250       cldemirad = cldemi
     3387      cldtaurad   = cldtau
     3388      cldtaupirad = cldtaupi
     3389      cldemirad   = cldemi
    32513390c
    32523391      if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND.
     
    32653404         beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
    32663405        endif
    3267         cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
    3268         cldemirad(i,k) = cldemi(i,k) * beta(i,k)
     3406        cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
     3407        cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
     3408        cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
    32693409       ENDDO
    32703410       ENDDO
     
    32873427          beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
    32883428         endif
    3289         cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
    3290         cldemirad(i,k) = cldemi(i,k) * beta(i,k)
     3429        cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
     3430        cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
     3431        cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
    32913432        endif
    32923433c
     
    33373478     s        topsw_aero, topsw0_aero,
    33383479     s        solsw_aero, solsw0_aero,
    3339      e        cldtaupi,
     3480     e        cldtaupirad,
    33403481     s        topswai_aero, solswai_aero)
    33413482           
     
    33513492       RCFC12 = RCFC12_act
    33523493c
     3494      IF (prt_level .GE.10) THEN
     3495       print *,' ->radlwsw, number 1 '
     3496      ENDIF
     3497c
    33533498         CALL radlwsw
    33543499     e        (dist, rmu0, fract,
     
    33563501     e        t_seri,q_seri,wo,
    33573502     e        cldfra, cldemirad, cldtaurad,
    3358      e        ok_ade, ok_aie,
     3503     e        ok_ade, ok_aie, flag_aerosol,
    33593504     e        tau_aero, piz_aero, cg_aero,
    3360      e        cldtaupi,new_aod,
     3505     e        cldtaupirad,new_aod,
    33613506     e        zqsat, flwc, fiwc,
    33623507     s        heat,heat0,cool,cool0,radsol,albpla,
     
    33883533       RCFC12 = RCFC12_per
    33893534c
     3535      IF (prt_level .GE.10) THEN
     3536       print *,' ->radlwsw, number 2 '
     3537      ENDIF
     3538c
    33903539         CALL radlwsw
    33913540     e        (dist, rmu0, fract,
     
    33933542     e        t_seri,q_seri,wo,
    33943543     e        cldfra, cldemi, cldtau,
    3395      e        ok_ade, ok_aie,
     3544     e        ok_ade, ok_aie, flag_aerosol,
    33963545     e        tau_aero, piz_aero, cg_aero,
    33973546     e        cldtaupi,new_aod,
     
    34793628c Appeler le programme de parametrisation de l'orographie
    34803629c a l'echelle sous-maille:
     3630c
     3631      IF (prt_level .GE.10) THEN
     3632       print *,' call orography ? ', ok_orodr
     3633      ENDIF
    34813634c
    34823635      IF (ok_orodr) THEN
     
    35693722
    35703723       IF (ok_hines) then
     3724
    35713725         CALL hines_gwd(klon,klev,dtime,paprs,pplay,
    35723726     i                  rlat,t_seri,u_seri,v_seri,
     
    35763730c  ajout des tendances
    35773731        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
     3732
    35783733      ENDIF
    3579 
     3734c
     3735
     3736c
     3737cIM cf. FLott BEG
    35803738C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
    35813739
     
    36023760cIM calcul composantes axiales du moment angulaire et couple des montagnes
    36033761c
    3604       IF (is_sequential .and. ok_orodr) THEN
    3605      
     3762      IF (is_sequential .and. ok_orodr) THEN
    36063763        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
    36073764     C                 ra,rg,romega,
     
    38984055c Convertir les incrementations en tendances
    38994056c
     4057      IF (prt_level .GE.10) THEN
     4058        print *,'Convertir les incrementations en tendances '
     4059      ENDIF
     4060c
    39004061      if (mydebug) then
    39014062        call writefield_phy('u_seri',u_seri,llm)
     
    40164177c=============================================================
    40174178
    4018       if (iflag_thermals>1) then
     4179      if (iflag_thermals>=1) then
    40194180      d_t_lscth=0.
    40204181      d_t_lscst=0.
  • LMDZ5/branches/testing/libf/phylmd/radlwsw.F90

    r1664 r1669  
    1010   t,q,wo,&
    1111   cldfra, cldemi, cldtaupd,&
    12    ok_ade, ok_aie,&
     12   ok_ade, ok_aie, flag_aerosol,&
    1313   tau_aero, piz_aero, cg_aero,&
    1414   cldtaupi, new_aod, &
     
    5656  ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
    5757  ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
     58  ! flag_aerosol-input-I- aerosol flag from 0 to 6
    5859  ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
    5960  ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
     
    119120
    120121  LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
     122  INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
    121123  REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
    122124  REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,9,2)                        ! aerosol optical properties (see aeropt.F)
     
    354356               zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
    355357               ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
    356                tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:),&
     358               tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
    357359               PTAUA, POMEGAA,&
    358360               ztopswadaero,zsolswadaero,&
    359361               ztopswaiaero,zsolswaiaero,&
    360                ok_ade, ok_aie)
     362               ok_ade, ok_aie, flag_aerosol)
    361363         
    362364       ELSE ! new_aod=T         
     
    377379               zsolsw_aero,zsolsw0_aero,&
    378380               ztopswcf_aero,zsolswcf_aero, &
    379                ok_ade, ok_aie)
    380          
     381               ok_ade, ok_aie, flag_aerosol)
    381382       ENDIF
    382383
  • LMDZ5/branches/testing/libf/phylmd/regr_lat_time_climoz_m.F90

    r1279 r1669  
    224224    ! Get the  number of months:
    225225    call nf95_inq_dimid(ncid_in, "time", dimid)
    226     call nf95_inquire_dimension(ncid_in, dimid, len=n_month)
     226    call nf95_inquire_dimension(ncid_in, dimid, nclen=n_month)
    227227
    228228    allocate(o3_in(n_lat, n_plev, n_month, read_climoz))
  • LMDZ5/branches/testing/libf/phylmd/sw_aeroAR4.F90

    r1664 r1669  
    1818     PSOLSWAERO,PSOLSW0AERO,&
    1919     PTOPSWCFAERO,PSOLSWCFAERO,&
    20      ok_ade, ok_aie )
     20     ok_ade, ok_aie, flag_aerosol )
    2121
    2222  USE dimphy
    23 
     23  USE phys_output_mod, ONLY : swaero_diag
    2424  IMPLICIT NONE
    2525
     
    5656  !     --------------
    5757  !        ORIGINAL : 89-07-14
    58   !        95-01-01   J.-J. MORCRETTE  Direct/Diffuse Albedo
    59   !        03-11-27   J. QUAAS Introduce aerosol forcings (based on BOUCHER)
    60   !        09-04      A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
     58  !        1995-01-01  J.-J. MORCRETTE  Direct/Diffuse Albedo
     59  !        2003-11-27  J. QUAAS Introduce aerosol forcings (based on BOUCHER)
     60  !        2009-04     A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
     61  !        2012-09     O. BOUCHER - reorganise aerosol cases with ok_ade, ok_aie, flag_aerosol
    6162  !     ------------------------------------------------------------------
    6263  !
     
    8283
    8384  REAL(KIND=8) PCLDSW(KDLON,KFLEV)    ! CLOUD FRACTION
    84   REAL(KIND=8) PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS
     85  REAL(KIND=8) PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
    8586  REAL(KIND=8) PCG(KDLON,2,KFLEV)     ! ASYMETRY FACTOR
    8687  REAL(KIND=8) POMEGA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
     
    132133  !$OMP THREADPRIVATE(initialized)
    133134
    134   !jq-Introduced for aerosol forcings
     135  !jq-local flag introduced for aerosol forcings
    135136  REAL(KIND=8), SAVE :: flag_aer
    136137  !$OMP THREADPRIVATE(flag_aer)
    137138
    138139  LOGICAL ok_ade, ok_aie    ! use aerosol forcings or not?
     140  INTEGER flag_aerosol      ! global flag for aerosol 0 (no aerosol) or 1-5 (aerosols)
    139141  REAL(KIND=8) tauaero(kdlon,kflev,9,2)  ! aerosol optical properties
    140142  REAL(KIND=8) pizaero(kdlon,kflev,9,2)  ! (see aeropt.F)
    141143  REAL(KIND=8) cgaero(kdlon,kflev,9,2)   ! -"-
    142   REAL(KIND=8) PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
     144  REAL(KIND=8) PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (present-day value)
    143145  REAL(KIND=8) POMEGAA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
    144146  REAL(KIND=8) PTOPSWADAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
    145147  REAL(KIND=8) PSOLSWADAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
    146   REAL(KIND=8) PTOPSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
    147   REAL(KIND=8) PSOLSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
     148  REAL(KIND=8) PTOPSWAD0AERO(KDLON)    ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
     149  REAL(KIND=8) PSOLSWAD0AERO(KDLON)    ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
    148150  REAL(KIND=8) PTOPSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
    149151  REAL(KIND=8) PSOLSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
    150   REAL(KIND=8) PTOPSWAERO(KDLON,9)      ! SW TOA AS DRF nat & ant
    151   REAL(KIND=8) PTOPSW0AERO(KDLON,9)      ! SW SRF AS DRF nat & ant
    152   REAL(KIND=8) PSOLSWAERO(KDLON,9)      ! SW TOA CS DRF nat & ant
    153   REAL(KIND=8) PSOLSW0AERO(KDLON,9)      ! SW SRF CS DRF nat & ant
     152  REAL(KIND=8) PTOPSWAERO(KDLON,9)    ! SW TOA AS DRF nat & ant
     153  REAL(KIND=8) PTOPSW0AERO(KDLON,9)    ! SW SRF AS DRF nat & ant
     154  REAL(KIND=8) PSOLSWAERO(KDLON,9)    ! SW TOA CS DRF nat & ant
     155  REAL(KIND=8) PSOLSW0AERO(KDLON,9)    ! SW SRF CS DRF nat & ant
    154156  REAL(KIND=8) PTOPSWCFAERO(KDLON,3)   !  SW TOA AS cloudRF nat & ant
    155157  REAL(KIND=8) PSOLSWCFAERO(KDLON,3)   !  SW SRF AS cloudRF nat & ant
     
    179181
    180182! Key to define the aerosol effect acting on climate
    181 ! 0: aerosol feedback active according to ok_ade, ok_aie  DEFAULT
    182 ! 1: no feedback , zero aerosol fluxes are used for climate, diagnostics according to ok_ade_ok_aie
    183 ! 2: feedback according to total aerosol direct effect used for climate, diagnostics according to ok_ade, ok_aie
    184 ! 3: feedback according to natural aerosol direct effect used for climate, diagnostics according to ok_ade_ok_aie
    185 
    186   INTEGER,SAVE :: AEROSOLFEEDBACK_ACTIVE = 0
     183! OB: AEROSOLFEEDBACK_ACTIVE is now a LOGICAL
     184! TRUE: fluxes use natural and/or anthropogenic aerosols according to ok_ade and ok_aie, DEFAULT
     185! FALSE: fluxes use no aerosols (case 1)
     186
     187  LOGICAL,SAVE :: AEROSOLFEEDBACK_ACTIVE = .TRUE.
    187188!$OMP THREADPRIVATE(AEROSOLFEEDBACK_ACTIVE) 
    188189
    189190      CHARACTER (LEN=20) :: modname='sw_aeroAR4'
    190191      CHARACTER (LEN=80) :: abort_message
    191 
    192   IF ((.not. ok_ade) .and. (AEROSOLFEEDBACK_ACTIVE .ge. 2)) THEN
    193      abort_message ='Error: direct effect is not activated but assumed to be active - see sw_aeroAR4.F90'
    194      CALL abort_gcm (modname,abort_message,1)
    195   ENDIF
    196   AEROSOLFEEDBACK_ACTIVE=MIN(MAX(AEROSOLFEEDBACK_ACTIVE,0),3)
    197   IF  (AEROSOLFEEDBACK_ACTIVE .gt. 3) THEN
    198      abort_message ='Error: AEROSOLFEEDBACK_ACTIVE options go only until 3'
    199      CALL abort_gcm (modname,abort_message,1)
    200   ENDIF
    201192
    202193  IF(.NOT.initialized) THEN
     
    209200     ALLOCATE(ZFSUPAI_AERO(KDLON,KFLEV+1))
    210201     ALLOCATE(ZFSDNAI_AERO(KDLON,KFLEV+1))
    211      ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,9))
    212      ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,9))
    213      ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,9))
    214      ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,9))
     202!-OB decrease size of these arrays to what is needed
     203!                | direct effect
     204!ind effect      | no aerosol   natural  total
     205!natural (PTAU)  |   1            3       2     --ZFSUP/ZFSDN
     206!total (PTAUA)   |                5       4     --ZFSUP/ZFSDN
     207!no cloud        |   1            3       2     --ZFSUP0/ZFSDN0
     208! so we need which case when ?
     209! ok_ade and ok_aie = 4-5, 4-2 and 2
     210! ok_ade and not ok_aie = 2-3 and 2
     211! not ok_ade and ok_aie = 5-3 and 5
     212! not ok_ade and not ok_aie = 3
     213! therefore the cases have the folliwng switches
     214! 3 = not ok_ade or not ok_aie
     215! 4 = ok_ade and ok_aie
     216! 2 = ok_ade
     217! 5 = ok_aie
     218     ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,5))
     219     ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,5))
     220     ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,3))
     221     ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,3))
     222! end OB modif
    215223     ZFSUPAD_AERO(:,:)=0.
    216224     ZFSDNAD_AERO(:,:)=0.
     
    226234
    227235  IF (appel1er) THEN
    228      WRITE(lunout,*) 'SW calling frequency : ', swpas
     236     WRITE(lunout,*)'SW calling frequency : ', swpas
    229237     WRITE(lunout,*) "   In general, it should be 1"
    230238     appel1er = .FALSE.
     
    241249     ENDDO
    242250
    243 ! clear sky is either computed IF no direct effect is asked for, or for extended diag)
    244      IF (( lev_histmth .ge. 4 ) .or. ( .not. ok_ade )) THEN   
     251! clear sky with no aerosols at all is computed IF ACTIVEFEEDBACK_ACTIVE is false or for extended diag
     252     IF ( swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE .OR. flag_aerosol .EQ. 0 ) THEN   
    245253
    246254     ! clear-sky: zero aerosol effect
     
    268276        ENDDO
    269277     ENDDO
    270      ENDIF
    271 
    272 ! cloudy sky is either computed IF no indirect effect is asked for, or for extended diag)
    273      IF (( lev_histmth .ge. 4 ) .or. ( .not. ok_aie )) THEN   
     278     ENDIF ! swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE
     279
     280! cloudy sky with no aerosols at all is either computed IF no indirect effect is asked for, or for extended diag
     281     IF ( swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE .OR. flag_aerosol .EQ. 0 ) THEN   
    274282     ! cloudy-sky: zero aerosol effect
    275283     flag_aer=0.0
     
    297305        ENDDO
    298306     ENDDO
    299      ENDIF
    300 
     307     ENDIF ! swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE
     308
     309     IF (flag_aerosol .GT. 0 ) THEN
     310
     311     IF (ok_ade.and.swaero_diag .or. .not. ok_ade) THEN
     312
     313        ! clear sky direct effect natural aerosol
     314        ! CAS AER (3)
     315        flag_aer=1.0
     316        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
     317             PRMU0,PFRAC,PTAVE,PWV,&
     318             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
     319        INU = 1
     320        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
     321             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
     322             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
     323             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
     324             ZFD, ZFU)
     325        INU = 2
     326        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
     327             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
     328             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
     329             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
     330             PWV, PQS,&
     331             ZFDOWN, ZFUP)
     332
     333        DO JK = 1 , KFLEV+1
     334           DO JL = 1, KDLON
     335              ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
     336              ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
     337           ENDDO
     338        ENDDO
     339     ENDIF !--end not swaero_diag or not ok_ade
    301340
    302341     IF (ok_ade) THEN
    303342
    304         ! clear sky (Anne Cozic 03/07/2007) direct effect of total aerosol
     343        ! clear sky direct effect of total aerosol
    305344        ! CAS AER (2)
    306345        flag_aer=1.0
     
    329368        ENDDO
    330369
    331 ! cloudy sky is either computed IF no indirect effect is asked for, or for extended diag)
    332         IF (( lev_histmth .ge. 2 ) .or. (.not. ok_aie)) THEN 
    333         ! cloudy-sky aerosol direct effect of total aerosol
     370        ! cloudy-sky with natural aerosols for indirect effect
     371        ! but total aerosols for direct effect
     372        ! PTAU
     373        ! CAS AER (2)
    334374        flag_aer=1.0
    335375        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
     
    356396           ENDDO
    357397        ENDDO
    358         ENDIF
    359 
    360 ! natural aeroosl clear sky is  computed  for extended diag)
    361         IF ( lev_histmth .ge. 4 ) THEN           
    362         ! clear sky direct effect natural aerosol
    363         flag_aer=1.0
    364         CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
    365              PRMU0,PFRAC,PTAVE,PWV,&
    366              ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
    367         INU = 1
    368         CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
    369              tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
    370              PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
    371              ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
    372              ZFD, ZFU)
    373         INU = 2
    374         CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
    375              tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
    376              ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
    377              ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
    378              PWV, PQS,&
    379              ZFDOWN, ZFUP)
    380 
    381         DO JK = 1 , KFLEV+1
    382            DO JL = 1, KDLON
    383               ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
    384               ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
    385            ENDDO
    386         ENDDO
    387         ENDIF
    388 
    389 ! cloud sky natural is for extended diagnostics
    390         IF ( lev_histmth .ge. 2 ) THEN
     398
     399     ENDIF !-end ok_ade
     400
     401     IF ( .not. ok_ade .or. .not. ok_aie ) THEN
     402
     403        ! cloudy-sky with natural aerosols for indirect effect
     404        ! and natural aerosols for direct effect
     405        ! PTAU
     406        ! CAS AER (3)
    391407        ! cloudy-sky direct effect natural aerosol
    392408        flag_aer=1.0
     
    414430           ENDDO
    415431        ENDDO
    416         ENDIF
    417 
    418      ENDIF ! ok_ade
    419 
    420 ! cloudy sky needs to be computed in all cases IF ok_aie is activated
    421      IF (ok_aie) THEN
    422         !jq   cloudy-sky + aerosol direct + aerosol indirect of total aerosol
     432
     433     ENDIF  !--true/false or false/true
     434
     435     IF (ok_ade .and. ok_aie) THEN
     436
     437        ! cloudy-sky with total aerosols for indirect effect
     438        ! and total aerosols for direct effect
     439        ! PTAUA
     440        ! CAS AER (2)
    423441        flag_aer=1.0
    424442        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
     
    438456             PWV, PQS,&
    439457             ZFDOWN, ZFUP)
     458
    440459        DO JK = 1 , KFLEV+1
    441460           DO JL = 1, KDLON
     
    444463           ENDDO
    445464        ENDDO
     465 
     466      ENDIF ! ok_ade .and. ok_aie
     467
     468     IF (ok_aie) THEN
     469        ! cloudy-sky with total aerosols for indirect effect
     470        ! and natural aerosols for direct effect
     471        ! PTAUA
     472        ! CAS AER (3)
     473        flag_aer=1.0
     474        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
     475             PRMU0,PFRAC,PTAVE,PWV,&
     476             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
     477        INU = 1
     478        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
     479             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
     480             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
     481             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
     482             ZFD, ZFU)
     483        INU = 2
     484        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
     485             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
     486             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
     487             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
     488             PWV, PQS,&
     489             ZFDOWN, ZFUP)
     490 
     491        DO JK = 1 , KFLEV+1
     492           DO JL = 1, KDLON
     493              ZFSUP_AERO(JL,JK,5) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
     494              ZFSDN_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
     495           ENDDO
     496        ENDDO
     497
    446498     ENDIF ! ok_aie     
    447499
     500     ENDIF !--if flag_aerosol GT 0
     501
    448502     itapsw = 0
    449503  ENDIF
    450504  itapsw = itapsw + 1
    451505
    452   IF  ( AEROSOLFEEDBACK_ACTIVE .eq. 0) THEN
     506  IF  ( AEROSOLFEEDBACK_ACTIVE .AND. flag_aerosol .GT. 0 ) THEN
    453507  IF ( ok_ade .and. ok_aie  ) THEN
    454508    ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
     
    457511    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,2)
    458512  ENDIF
     513
    459514  IF ( ok_ade .and. (.not. ok_aie) )  THEN
    460515    ZFSUP(:,:) =    ZFSUP_AERO(:,:,2)
     
    465520
    466521  IF ( (.not. ok_ade) .and. ok_aie  )  THEN
    467     print*,'Warning: indirect effect in cloudy regions includes direct aerosol effect'
    468     ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
    469     ZFSDN(:,:) =    ZFSDN_AERO(:,:,4)
    470     ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
    471     ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
    472   ENDIF
     522    ZFSUP(:,:) =    ZFSUP_AERO(:,:,5)
     523    ZFSDN(:,:) =    ZFSDN_AERO(:,:,5)
     524    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,3)
     525    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,3)
     526  ENDIF
     527
    473528  IF ((.not. ok_ade) .and. (.not. ok_aie)) THEN
     529    ZFSUP(:,:) =    ZFSUP_AERO(:,:,3)
     530    ZFSDN(:,:) =    ZFSDN_AERO(:,:,3)
     531    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,3)
     532    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,3)
     533  ENDIF
     534
     535! MS the following allows to compute the forcing diagostics without
     536! letting the aerosol forcing act on the meteorology
     537! SEE logic above
     538  ELSE
    474539    ZFSUP(:,:) =    ZFSUP_AERO(:,:,1)
    475540    ZFSDN(:,:) =    ZFSDN_AERO(:,:,1)
     
    478543  ENDIF
    479544
    480 ! MS the following allows to compute the forcing diagostics without
    481 ! letting the aerosol forcing act on the meteorology
    482 ! SEE logic above
    483   ELSEIF  ( AEROSOLFEEDBACK_ACTIVE .gt. 0) THEN
    484     ZFSUP(:,:) =    ZFSUP_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
    485     ZFSDN(:,:) =    ZFSDN_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
    486     ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
    487     ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
    488   ENDIF
    489  
    490 
     545! Now computes heating rates
    491546  DO k = 1, KFLEV
    492547     kpl1 = k+1
     
    511566     PTOPSW(i) = ZFSDN(i,KFLEV+1) - ZFSUP(i,KFLEV+1)
    512567
    513 
    514568! net anthropogenic forcing direct and 1st indirect effect diagnostics
    515569! requires a natural aerosol field read and used
    516570! Difference of net fluxes from double call to radiation
    517571
    518 
    519572IF (ok_ade) THEN
    520573
    521574! indices 1: natural; 2 anthropogenic
     575
    522576! TOA/SRF all sky natural forcing
    523577     PSOLSWAERO(i,1) = (ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))-(ZFSDN_AERO(i,1,1) - ZFSUP_AERO(i,1,1))
    524578     PTOPSWAERO(i,1) = (ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))- (ZFSDN_AERO(i,KFLEV+1,1) - ZFSUP_AERO(i,KFLEV+1,1))
    525579
     580! TOA/SRF clear sky natural forcing
     581     PSOLSW0AERO(i,1) = (ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))-(ZFSDN0_AERO(i,1,1) - ZFSUP0_AERO(i,1,1))
     582     PTOPSW0AERO(i,1) = (ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))-(ZFSDN0_AERO(i,KFLEV+1,1) - ZFSUP0_AERO(i,KFLEV+1,1))
     583
     584   IF (ok_aie) THEN
     585
     586! TOA/SRF all sky anthropogenic forcing
     587     PSOLSWAERO(i,2) = (ZFSDN_AERO(i,1,4) - ZFSUP_AERO(i,1,4))-(ZFSDN_AERO(i,1,5) - ZFSUP_AERO(i,1,5))
     588     PTOPSWAERO(i,2) = (ZFSDN_AERO(i,KFLEV+1,4) - ZFSUP_AERO(i,KFLEV+1,4))- (ZFSDN_AERO(i,KFLEV+1,5) - ZFSUP_AERO(i,KFLEV+1,5))
     589
     590   ELSE
     591
    526592! TOA/SRF all sky anthropogenic forcing
    527593     PSOLSWAERO(i,2) = (ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))-(ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))
    528594     PTOPSWAERO(i,2) = (ZFSDN_AERO(i,KFLEV+1,2) - ZFSUP_AERO(i,KFLEV+1,2))- (ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))
    529595
    530 ! TOA/SRF clear sky natural forcing
    531      PSOLSW0AERO(i,1) = (ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))-(ZFSDN0_AERO(i,1,1) - ZFSUP0_AERO(i,1,1))
    532      PTOPSW0AERO(i,1) = (ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))-(ZFSDN0_AERO(i,KFLEV+1,1) - ZFSUP0_AERO(i,KFLEV+1,1))
     596   ENDIF
    533597
    534598! TOA/SRF clear sky anthropogenic forcing
     
    536600     PTOPSW0AERO(i,2) = (ZFSDN0_AERO(i,KFLEV+1,2) - ZFSUP0_AERO(i,KFLEV+1,2))-(ZFSDN0_AERO(i,KFLEV+1,3) - ZFSUP0_AERO(i,KFLEV+1,3))
    537601
     602! direct anthropogenic forcing , as in old LMDzT, however differences of net fluxes
     603     PSOLSWADAERO(i) = PSOLSWAERO(i,2)
     604     PTOPSWADAERO(i) = PTOPSWAERO(i,2)
     605     PSOLSWAD0AERO(i) = PSOLSW0AERO(i,2)
     606     PTOPSWAD0AERO(i) = PTOPSW0AERO(i,2)
     607
     608! OB: these diagnostics may not always work but who need them
    538609! Cloud forcing indices 1: natural; 2 anthropogenic; 3: zero aerosol direct effect
    539610! Instantaneously computed cloudy sky direct aerosol effect, cloud forcing due to aerosols above clouds
     
    552623     PTOPSWCFAERO(i,3) = (ZFSDN_AERO(i,KFLEV+1,1) - ZFSUP_AERO(i,KFLEV+1,1))- (ZFSDN0_AERO(i,KFLEV+1,1) - ZFSUP0_AERO(i,KFLEV+1,1))
    553624
    554 ! direct anthropogenic forcing , as in old LMDzT, however differences of net fluxes
    555      PSOLSWADAERO(i) = PSOLSWAERO(i,2)
    556      PTOPSWADAERO(i) = PTOPSWAERO(i,2)
    557      PSOLSWAD0AERO(i) = PSOLSW0AERO(i,2)
    558      PTOPSWAD0AERO(i) = PTOPSW0AERO(i,2)
    559 
    560625ENDIF
    561626
    562 
    563627IF (ok_aie) THEN
     628   IF (ok_ade) THEN
    564629     PSOLSWAIAERO(i) = (ZFSDN_AERO(i,1,4) - ZFSUP_AERO(i,1,4))-(ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))
    565630     PTOPSWAIAERO(i) = (ZFSDN_AERO(i,KFLEV+1,4) - ZFSUP_AERO(i,KFLEV+1,4))-(ZFSDN_AERO(i,KFLEV+1,2) - ZFSUP_AERO(i,KFLEV+1,2))
     631   ELSE
     632     PSOLSWAIAERO(i) = (ZFSDN_AERO(i,1,5) - ZFSUP_AERO(i,1,5))-(ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))
     633     PTOPSWAIAERO(i) = (ZFSDN_AERO(i,KFLEV+1,5) - ZFSUP_AERO(i,KFLEV+1,5))-(ZFSDN_AERO(i,KFLEV+1,3) - ZFSUP_AERO(i,KFLEV+1,3))
     634   ENDIF
    566635ENDIF
    567636
    568   ENDDO
     637ENDDO
     638
    569639END SUBROUTINE SW_AEROAR4
  • LMDZ5/branches/testing/libf/phylmd/thermcell.h

    r1496 r1669  
    11      integer            :: iflag_thermals,nsplit_thermals
     2
     3!!! nrlmd le 10/04/2012
     4      integer            :: iflag_trig_bl,iflag_clos_bl
     5      integer            :: tau_trig_shallow,tau_trig_deep
     6      real               :: s_trig
     7!!! fin nrlmd le 10/04/2012
     8
    29      real,parameter     :: r_aspect_thermals=2.,l_mix_thermals=30.
    310      real               :: alp_bl_k
    411      real               :: tau_thermals
    5       integer,parameter  :: w2di_thermals=1
     12      integer,parameter  :: w2di_thermals=0
    613      integer            :: isplit
    714
     
    1421      common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
    1522
     23!!! nrlmd le 10/04/2012
     24      common/ctherm6/iflag_trig_bl,iflag_clos_bl
     25      common/ctherm7/tau_trig_shallow,tau_trig_deep
     26      common/ctherm8/s_trig
     27!!! fin nrlmd le 10/04/2012
     28
    1629!$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
     30!$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
  • LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90

    r1525 r1669  
    1010     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
    1111     &                  ,zmax0, f0,zw2,fraca,ztv &
    12      &                  ,zpspsk,ztla,zthl)
     12     &                  ,zpspsk,ztla,zthl &
     13!!! nrlmd le 10/04/2012
     14     &                  ,pbl_tke,pctsrf,omega,airephy &
     15     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
     16     &                  ,n2,s2,ale_bl_stat &
     17     &                  ,therm_tke_max,env_tke_max &
     18     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
     19     &                  ,alp_bl_conv,alp_bl_stat &
     20!!! fin nrlmd le 10/04/2012
     21     &                         )
    1322
    1423      USE dimphy
     
    4756#include "iniprint.h"
    4857#include "thermcell.h"
     58!!! nrlmd le 10/04/2012
     59#include "indicesol.h"
     60!!! fin nrlmd le 10/04/2012
    4961
    5062!   arguments:
     
    7789      integer,save :: lev_out=10
    7890!$OMP THREADPRIVATE(lev_out)
     91
     92      REAL susqr2pi, Reuler
    7993
    8094      INTEGER ig,k,l,ll,ierr
     
    155169       real seuil
    156170      real csc(klon,klev)
     171
     172!!! nrlmd le 10/04/2012
     173
     174!------Entrées
     175      real pbl_tke(klon,klev+1,nbsrf)
     176      real pctsrf(klon,nbsrf)
     177      real omega(klon,klev)
     178      real airephy(klon)
     179!------Sorties
     180      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
     181      real therm_tke_max0(klon),env_tke_max0(klon)
     182      real n2(klon),s2(klon)
     183      real ale_bl_stat(klon)
     184      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
     185      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
     186!------Local
     187      integer nsrf
     188      real rhobarz0(klon)                    ! Densité au LCL
     189      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
     190      integer klcl(klon)                     ! Niveau du LCL
     191      real interp(klon)                      ! Coef d'interpolation pour le LCL
     192!--Triggering
     193      real Su                                ! Surface unité: celle d'un updraft élémentaire
     194      parameter(Su=4e4)
     195      real hcoef                             ! Coefficient directeur pour le calcul de s2
     196      parameter(hcoef=1)
     197      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
     198      parameter(hmincoef=0.3)
     199      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
     200      parameter(eps1=0.3)
     201      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
     202      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
     203      real zmax_moy_coef
     204      parameter(zmax_moy_coef=0.33)
     205      real depth(klon)                       ! Epaisseur moyenne du cumulus
     206      real w_max(klon)                       ! Vitesse max statistique
     207      real s_max(klon)
     208!--Closure
     209      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne
     210      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
     211      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
     212      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
     213      parameter(coef_m=1.)
     214      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
     215      parameter(coef_tke=1.)
     216
     217!!! fin nrlmd le 10/04/2012
    157218
    158219!
     
    679740      enddo
    680741!
     742
     743!!! nrlmd le 10/04/2012
     744
     745!------------Test sur le LCL des thermiques
     746    do ig=1,ngrid
     747      ok_lcl(ig)=.false.
     748      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
     749    enddo
     750
     751!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
     752    do l=1,nlay-1
     753      do ig=1,ngrid
     754        if (ok_lcl(ig)) then
     755          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
     756          klcl(ig)=l
     757          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
     758          endif
     759        endif
     760      enddo
     761    enddo
     762
     763!------------Hauteur des thermiques
     764!!jyg le 27/04/2012
     765!!    do ig =1,ngrid
     766!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
     767!! &               -rhobarz(ig,klcl(ig)))*interp(ig)
     768!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
     769!!    zmax(ig)=pphi(ig,lmax(ig))/rg
     770!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
     771!!    enddo
     772    do ig =1,ngrid
     773     zmax(ig)=pphi(ig,lmax(ig))/rg
     774     if (ok_lcl(ig)) then
     775      rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
     776 &               -rhobarz(ig,klcl(ig)))*interp(ig)
     777      zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
     778      zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
     779     else
     780      rhobarz0(ig)=0.
     781      zlcl(ig)=zmax(ig)
     782     endif
     783    enddo
     784!!jyg fin
     785
     786!------------Calcul des propriétés du thermique au LCL
     787  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
     788
     789  !-----Initialisation de la TKE moyenne
     790   do l=1,nlay
     791    do ig=1,ngrid
     792     pbl_tke_max(ig,l)=0.
     793    enddo
     794   enddo
     795
     796!-----Calcul de la TKE moyenne
     797   do nsrf=1,nbsrf
     798    do l=1,nlay
     799     do ig=1,ngrid
     800     pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
     801     enddo
     802    enddo
     803   enddo
     804
     805!-----Initialisations des TKE dans et hors des thermiques
     806   do l=1,nlay
     807    do ig=1,ngrid
     808    therm_tke_max(ig,l)=pbl_tke_max(ig,l)
     809    env_tke_max(ig,l)=pbl_tke_max(ig,l)
     810    enddo
     811   enddo
     812
     813!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
     814   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
     815  &           rg,pplev,therm_tke_max)
     816!   print *,' thermcell_tke_transport -> '   !!jyg
     817
     818!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
     819   do l=1,nlay
     820    do ig=1,ngrid
     821     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
     822     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
     823     w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
     824    enddo
     825   enddo
     826!    print *,' apres w_ls = '   !!jyg
     827
     828  do ig=1,ngrid
     829   if (ok_lcl(ig)) then
     830     fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
     831 &             -fraca(ig,klcl(ig)))*interp(ig)
     832     w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
     833 &         -zw2(ig,klcl(ig)))*interp(ig)
     834     w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
     835 &             -w_ls(ig,klcl(ig)))*interp(ig)
     836     therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
     837 &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
     838     env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
     839 &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
     840     pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
     841 &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
     842     if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
     843     if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
     844     if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
     845   else
     846     fraca0(ig)=0.
     847     w0(ig)=0.
     848!!jyg le 27/04/2012
     849!!     zlcl(ig)=0.
     850!!
     851   endif
     852  enddo
     853
     854  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
     855!  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
     856
     857!------------Triggering------------------
     858  IF (iflag_trig_bl.ge.1) THEN
     859
     860!-----Initialisations
     861   depth(:)=0.
     862   n2(:)=0.
     863   s2(:)=0.
     864   s_max(:)=0.
     865
     866!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
     867   do ig=1,ngrid
     868     zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
     869     depth(ig)=zmax_moy(ig)-zlcl(ig)
     870     hmin(ig)=hmincoef*zlcl(ig)
     871     if (depth(ig).ge.10.) then
     872       s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
     873       n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
     874!!
     875!!jyg le 27/04/2012
     876!!       s_max(ig)=s2(ig)*log(n2(ig))
     877!!       if (n2(ig) .lt. 1) s_max(ig)=0.
     878       s_max(ig)=s2(ig)*log(max(n2(ig),1.))
     879!!fin jyg
     880     else
     881       s2(ig)=0.
     882       n2(ig)=0.
     883       s_max(ig)=0.
     884     endif
     885   enddo
     886!   print *,'avant Calcul de Wmax '    !!jyg
     887
     888!-----Calcul de Wmax et ALE_BL_STAT associée
     889!!jyg le 30/04/2012
     890!!   do ig=1,ngrid
     891!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
     892!!     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))))
     893!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
     894!!     else
     895!!     w_max(ig)=0.
     896!!     ale_bl_stat(ig)=0.
     897!!     endif
     898!!   enddo
     899   susqr2pi=su*sqrt(2.*Rpi)
     900   Reuler=exp(1.)
     901   do ig=1,ngrid
     902     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
     903      w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
     904      ale_bl_stat(ig)=0.5*w_max(ig)**2
     905     else
     906      w_max(ig)=0.
     907      ale_bl_stat(ig)=0.
     908     endif
     909   enddo
     910
     911  ENDIF ! iflag_trig_bl
     912!  print *,'ENDIF  iflag_trig_bl'    !!jyg
     913
     914!------------Closure------------------
     915
     916  IF (iflag_clos_bl.ge.1) THEN
     917
     918!-----Calcul de ALP_BL_STAT
     919  do ig=1,ngrid
     920  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
     921  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
     922 &                   (w0(ig)**2)
     923  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
     924 &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
     925    if (iflag_clos_bl.ge.2) then
     926    alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
     927 &                   (w0(ig)**2)
     928    else
     929    alp_bl_conv(ig)=0.
     930    endif
     931  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
     932  enddo
     933
     934!-----Sécurité ALP infinie
     935  do ig=1,ngrid
     936   if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
     937  enddo
     938
     939  ENDIF ! (iflag_clos_bl.ge.1)
     940
     941!!! fin nrlmd le 10/04/2012
     942
    681943      if (prt_level.ge.10) then
    682944         ig=igout
     
    8581120      end
    8591121
     1122!!! nrlmd le 10/04/2012                          Transport de la TKE par le thermique moyen pour la fermeture en ALP
     1123!                                                         On transporte pbl_tke pour donner therm_tke
     1124!                                          Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
     1125      subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
     1126     &           rg,pplev,therm_tke_max)
     1127      implicit none
     1128
     1129#include "iniprint.h"
     1130!=======================================================================
     1131!
     1132!   Calcul du transport verticale dans la couche limite en presence
     1133!   de "thermiques" explicitement representes
     1134!   calcul du dq/dt une fois qu'on connait les ascendances
     1135!
     1136!=======================================================================
     1137
     1138      integer ngrid,nlay,nsrf
     1139
     1140      real ptimestep
     1141      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
     1142      real entr0(ngrid,nlay),rg
     1143      real therm_tke_max(ngrid,nlay)
     1144      real detr0(ngrid,nlay)
     1145
     1146
     1147      real masse(ngrid,nlay),fm(ngrid,nlay+1)
     1148      real entr(ngrid,nlay)
     1149      real q(ngrid,nlay)
     1150      integer lev_out                           ! niveau pour les print
     1151
     1152      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
     1153
     1154      real zzm
     1155
     1156      integer ig,k
     1157      integer isrf
     1158
     1159
     1160      lev_out=0
     1161
     1162
     1163      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
     1164
     1165!   calcul du detrainement
     1166      do k=1,nlay
     1167         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
     1168         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
     1169      enddo
     1170
     1171
     1172! Decalage vertical des entrainements et detrainements.
     1173      masse(:,1)=0.5*masse0(:,1)
     1174      entr(:,1)=0.5*entr0(:,1)
     1175      detr(:,1)=0.5*detr0(:,1)
     1176      fm(:,1)=0.
     1177      do k=1,nlay-1
     1178         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
     1179         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
     1180         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
     1181         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
     1182      enddo
     1183      fm(:,nlay+1)=0.
     1184
     1185!!! nrlmd le 16/09/2010
     1186!   calcul de la valeur dans les ascendances
     1187!       do ig=1,ngrid
     1188!          qa(ig,1)=q(ig,1)
     1189!       enddo
     1190!!!
     1191
     1192!do isrf=1,nsrf
     1193
     1194!   q(:,:)=therm_tke(:,:,isrf)
     1195   q(:,:)=therm_tke_max(:,:)
     1196!!! nrlmd le 16/09/2010
     1197      do ig=1,ngrid
     1198         qa(ig,1)=q(ig,1)
     1199      enddo
     1200!!!
     1201
     1202    if (1==1) then
     1203      do k=2,nlay
     1204         do ig=1,ngrid
     1205            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     1206     &         1.e-5*masse(ig,k)) then
     1207         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     1208     &         /(fm(ig,k+1)+detr(ig,k))
     1209            else
     1210               qa(ig,k)=q(ig,k)
     1211            endif
     1212            if (qa(ig,k).lt.0.) then
     1213!               print*,'qa<0!!!'
     1214            endif
     1215            if (q(ig,k).lt.0.) then
     1216!               print*,'q<0!!!'
     1217            endif
     1218         enddo
     1219      enddo
     1220
     1221! Calcul du flux subsident
     1222
     1223      do k=2,nlay
     1224         do ig=1,ngrid
     1225            wqd(ig,k)=fm(ig,k)*q(ig,k)
     1226            if (wqd(ig,k).lt.0.) then
     1227!               print*,'wqd<0!!!'
     1228            endif
     1229         enddo
     1230      enddo
     1231      do ig=1,ngrid
     1232         wqd(ig,1)=0.
     1233         wqd(ig,nlay+1)=0.
     1234      enddo
     1235
     1236! Calcul des tendances
     1237      do k=1,nlay
     1238         do ig=1,ngrid
     1239            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
     1240     &               -wqd(ig,k)+wqd(ig,k+1))  &
     1241     &               *ptimestep/masse(ig,k)
     1242         enddo
     1243      enddo
     1244
     1245 endif
     1246
     1247   therm_tke_max(:,:)=q(:,:)
     1248
     1249      return
     1250!!! fin nrlmd le 10/04/2012
     1251     end
     1252
  • LMDZ5/branches/testing/libf/phylmd/wake.F

    r1403 r1669  
    561561      ENDDO
    562562
    563 c      On evite kupper = 1
     563c      On evite kupper = 1 et kupper = klev
    564564      DO i=1,klon
    565565        kupper(i) = max(kupper(i),2)
     566        kupper(i) = min(kupper(i),klev-1)
    566567      ENDDO
    567568
Note: See TracChangeset for help on using the changeset viewer.