MODULE lmdz_thermcell_main

! $Id: lmdz_thermcell_main.F90 4590 2023-06-29 01:03:15Z fhourdin $

!

CONTAINS

     &                  ,pplay,pplev,pphi,debut  &

     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &

     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &

     &                  ,ratqscth,ratqsdiff,zqsatth  &

     &                  ,zpspsk,ztla,zthl,ztva &

#ifdef ISO

     &      ,xtpo,xtpdoadj &

#endif

     &   )

      USE lmdz_thermcell_ini, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level

      USE lmdz_thermcell_ini, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals

      USE lmdz_thermcell_ini, ONLY: iflag_thermals_down,fact_thermals_down

      USE lmdz_thermcell_ini, ONLY: RD,RG

      USE lmdz_thermcell_down, ONLY: thermcell_updown_dq

      USE lmdz_thermcell_closure, ONLY: thermcell_closure

      USE lmdz_thermcell_dq, ONLY: thermcell_dq

      USE lmdz_thermcell_dry, ONLY: thermcell_dry

      USE lmdz_thermcell_dv2, ONLY: thermcell_dv2

      USE lmdz_thermcell_env, ONLY: thermcell_env

      USE lmdz_thermcell_flux2, ONLY: thermcell_flux2

      USE lmdz_thermcell_height, ONLY: thermcell_height

      USE lmdz_thermcell_plume, ONLY: thermcell_plume

      USE lmdz_thermcell_plume_6A, ONLY: thermcell_plume_6A,thermcell_plume_5B

#ifdef ISO

  USE infotrac_phy, ONLY : ntiso

#ifdef ISOVERIF

  USE isotopes_mod, ONLY : iso_eau,iso_HDO

  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &

        iso_verif_aberrant_encadre

#endif

#endif

      IMPLICIT NONE

!=======================================================================

!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu

!   Version du 09.02.07

!   Calcul du transport vertical dans la couche limite en presence

!   de "thermiques" explicitement representes avec processus nuageux

!

!   Reecriture a partir d'un listing papier a Habas, le 14/02/00

!

!   le thermique est suppose homogene et dissipe par melange avec

!   son environnement. la longueur l_mix controle l'efficacite du

!   melange

!

!   Le calcul du transport des differentes especes se fait en prenant

!   en compte:

!     1. un flux de masse montant

!     2. un flux de masse descendant

!     3. un entrainement

!     4. un detrainement

!

! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)

!    Introduction of an implicit computation of vertical advection in

!    the environment of thermal plumes in thermcell_dq

!    impl =     0 : explicit, 1 : implicit, -1 : old version

!    controled by iflag_thermals =

!       15, 16 run with impl=-1 : numerical convergence with NPv3

!       17, 18 run with impl=1  : more stable

!    15 and 17 correspond to the activation of the stratocumulus "bidouille"

!

! Using

!    abort_physic

!    iso_verif_aberrant_encadre

!    iso_verif_egalite

!    test_ltherm

!    thermcell_closure

!    thermcell_dq

!    thermcell_dry

!    thermcell_dv2

!    thermcell_env

!    thermcell_flux2

!    thermcell_height

!    thermcell_plume

!    thermcell_plume_5B

!    thermcell_plume_6A

!

!=======================================================================

!-----------------------------------------------------------------------

!   declarations:

!   -------------

!   arguments:

!   ----------

      integer, intent(in) :: itap,ngrid,nlay

      real, intent(in) ::  ptimestep

      real, intent(in), dimension(ngrid,nlay)    :: pt,pu,pv,pplay,pphi

! ATTENTION : po et zpspsk sont inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)

      real, intent(inout), dimension(ngrid,nlay)    :: po

      real, intent(out), dimension(ngrid,nlay)    :: zpspsk

      real, intent(in), dimension(ngrid,nlay+1)  :: pplev

      integer, intent(out), dimension(ngrid) :: lmax

      real, intent(out), dimension(ngrid,nlay)   :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0

      real, intent(out), dimension(ngrid,nlay)   :: ztla,zqla,zqta,zqsatth,zthl

      real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca

      real, intent(inout), dimension(ngrid) :: zmax0,f0

      real, intent(out), dimension(ngrid,nlay) :: ztva,ztv

      logical, intent(in) :: debut

      real,intent(out), dimension(ngrid,nlay) :: ratqscth,ratqsdiff

      real, intent(out), dimension(ngrid) :: pcon

      real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3

      real, intent(out), dimension(ngrid) :: wmax_sec

      integer,intent(out), dimension(ngrid) :: lalim

      real, intent(out), dimension(ngrid,nlay+1) :: fm

      real, intent(out), dimension(ngrid,nlay) :: alim_star

      real, intent(out), dimension(ngrid) :: zmax

!   local:

!   ------

      integer,save :: igout=1

!$OMP THREADPRIVATE(igout)

      integer,save :: lunout1=6

!$OMP THREADPRIVATE(lunout1)

      integer,save :: lev_out=10

!$OMP THREADPRIVATE(lev_out)

      real lambda, zf,zf2,var,vardiff,CHI

      integer ig,k,l,ierr,ll

      logical sorties

      character (len=20) :: modname='thermcell_main'

      character (len=80) :: abort_message

#ifdef ISO

      REAL xtpo(ntiso,ngrid,nlay),xtpdoadj(ntiso,ngrid,nlay)

      REAL xtzo(ntiso,ngrid,nlay)

      REAL xtpdoadj_tmp(ngrid,nlay)

      REAL xtpo_tmp(ngrid,nlay)

      REAL xtzo_tmp(ngrid,nlay)

      integer ixt

#endif

!

!-----------------------------------------------------------------------

!   initialisation:

!   ---------------

!

       sorties=.true.

      IF(ngrid.NE.ngrid) THEN

         PRINT*

         PRINT*,'STOP dans convadj'

         PRINT*,'ngrid    =',ngrid

         PRINT*,'ngrid  =',ngrid

      ENDIF

!

!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'

!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2

     enddo

       enddo

      endif

!-----------------------------------------------------------------------

! Calcul de T,q,ql a partir de Tl et qT dans l environnement

!   --------------------------------------------------------------------

!

      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &

!------------------------------------------------------------------------

!                       --------------------

!

!

!                       + + + + + + + + + + +

!

!

!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz

!  wh,wt,wo ...

!

!                       + + + + + + + + + + +  zh,zu,zv,zo,rho

!

!

!                       --------------------   zlev(1)

!                       \\\\\\\\\\\\\\\\\\\\

!

!

!-----------------------------------------------------------------------

!   Calcul des altitudes des couches

!-----------------------------------------------------------------------

      enddo

      enddo

      enddo

!-----------------------------------------------------------------------

!   Calcul des densites et masses

!-----------------------------------------------------------------------

      enddo

      enddo

!------------------------------------------------------------------

!

!             /|\

!    --------  |  F_k+1 -------

!                              ----> D_k

!             /|\              <---- E_k , A_k

!    --------  |  F_k ---------

!                              ----> D_k-1

!                              <---- E_k-1 , A_k-1

!

!

!

!

!

!    ---------------------------

!

!    ----- F_lmax+1=0 ----------         \

!            lmax     (zmax)              |

!    ---------------------------          |

!                                         |

!    ---------------------------          |

!                                         |

!    ---------------------------          |

!                                         |

!    ---------------------------          |

!                                         |

!    ---------------------------          |

!                                         |  E

!    ---------------------------          |  D

!                                         |

!    ---------------------------          |

!                                         |

!    ---------------------------  \       |

!            lalim                 |      |

!    ---------------------------   |      |

!                                  |      |

!    ---------------------------   |      |

!                                  | A    |

!    ---------------------------   |      |

!                                  |      |

!    ---------------------------   |      |

!    lmin  (=1 pour le moment)     |      |

!    ----- F_lmin=0 ------------  /      /

!

!    ---------------------------

!    //////////////////////////

!

!

!=============================================================================

!  Calculs initiaux ne faisant pas intervenir les changements de phase

!=============================================================================

!------------------------------------------------------------------

!  1. alim_star est le profil vertical de l'alimentation a la base du

!     panache thermique, calcule a partir de la flotabilite de l'air sec

!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star

!------------------------------------------------------------------

!

!-----------------------------------------------------------------------------

!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un

!     panache sec conservatif (e=d=0) alimente selon alim_star

!     Il s'agit d'un calcul de type CAPE

!     zmax_sec est utilise pour determiner la geometrie du thermique.

!------------------------------------------------------------------------------

!---------------------------------------------------------------------------------

!calcul du melange et des variables dans le thermique

!--------------------------------------------------------------------------------

!

!=====================================================================

! Old version of thermcell_plume in thermcell_plume_6A.F90

! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding

! to the 5B and 6A versions used for CMIP5 and CMIP6.

! The latest was previously named thermcellV1_plume.

! The new thermcell_plume is a clean version (removing obsolete

! options) of thermcell_plume_6A.

! The 3 versions are controled by

! flag_thermals_ed <= 9 thermcell_plume_6A

!                  <= 19 thermcell_plume_5B

!                  else thermcell_plume (default 20 for convergence with 6A)

! Fredho

!=====================================================================

!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'

         CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&

     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &

     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &

     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &

!        print*,'THERM RIO et al 2010, version d Arnaud'

         CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&

     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &

     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &

     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &

      else

         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&

     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &

     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &

     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &

      endif

      endif

!-------------------------------------------------------------------------------

! Calcul des caracteristiques du thermique:zmax,zmix,wmax

!-------------------------------------------------------------------------------

!

      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &

! Attention, w2 est transforme en sa racine carree dans cette routine

! Le probleme vient du fait que linter et lmix sont souvent egaux a 1.

      enddo

!     print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax

!-------------------------------------------------------------------------------

! Fermeture,determination de f

!-------------------------------------------------------------------------------

!

!

      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &

      endif

! Choix de la fonction d'alimentation utilisee pour la fermeture.

! Apparemment sans importance

!

!CR Appel de la fermeture seche

     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Appel avec les zmax et wmax tenant compte de la condensation

! Semble moins bien marcher

     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &

     endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      else

      endif

! Test valable seulement en 1D mais pas genant

      endif

!-------------------------------------------------------------------------------

!deduction des flux

      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &

     &       lalim,lmax,alim_star,  &

     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &

!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)

!------------------------------------------------------------------

!   On ne prend pas directement les profils issus des calculs precedents

!   mais on s'autorise genereusement une relaxation vers ceci avec

!   une constante de temps tau_thermals (typiquement 1800s).

!------------------------------------------------------------------

      else

      endif

!------------------------------------------------------------------

! Calcul de la fraction de l'ascendance

!------------------------------------------------------------------

      enddo

            else

            endif

         enddo

      enddo

!c------------------------------------------------------------------

!   calcul du transport vertical

!------------------------------------------------------------------

        ! we want to transport potential temperature, total water and momentum

      ELSE

      !--------------------------------------------------------------

        call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &

        call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &

#ifdef ISO

        ! C Risi: on utilise directement la meme routine

        do ixt=1,ntiso

          do ll=1,nlay

            DO ig=1,ngrid

                xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)

                xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)

            enddo

          enddo

          call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &

     &                   xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)

          do ll=1,nlay

            DO ig=1,ngrid

                xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)

            enddo

          enddo

        enddo

#endif

#ifdef ISO

#ifdef ISOVERIF

      DO  ll=1,nlay

        DO ig=1,ngrid

          if (iso_eau.gt.0) then

              call iso_verif_egalite(xtpo(iso_eau,ig,ll), &

     &          po(ig,ll),'thermcell_main 594')

              call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &

     &          pdoadj(ig,ll),'thermcell_main 596')

          endif

          if (iso_HDO.gt.0) then

              call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &

     &           /po(ig,ll),'thermcell_main 610')

          endif

        enddo

      enddo !DO  ll=1,nlay

      write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq'

#endif

#endif

!------------------------------------------------------------------

!  calcul du transport vertical du moment horizontal

!------------------------------------------------------------------

!IM 090508

! Calcul du transport de V tenant compte d'echange par gradient

! de pression horizontal avec l'environnement

         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &

!    &    ,fraca*dvdq,zmax &

     &    ,fraca,zmax &

      else

! calcul purement conservatif pour le transport de V

         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &

         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &

      endif

    ENDIF

!     print*,'13 OK convect8'

         enddo

      enddo

!------------------------------------------------------------------

!   Calculs de diagnostiques pour les sorties

!------------------------------------------------------------------

!calcul de fraca pour les sorties

      if (sorties) then

! calcul du niveau de condensation

! initialisation

      enddo

!nouveau calcul

      enddo

!IM   do k=1,nlay

         if ((pcon(ig).le.pplay(ig,k))  &

         endif

         enddo

      enddo

!IM

      ierr=0

           ierr=1

        endif

      enddo

      endif

            endif

         enddo

      enddo

!calcul des moments

!initialisation

         enddo

      enddo

!

            else

            endif

            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &

!test: on calcul q2/po=ratqsc

         enddo

      enddo

!calcul des flux: q, thetal et thetav

         enddo

      enddo

!calcul du ratqscdiff

      var=0.

      vardiff=0.

            endif

         enddo

      enddo

               zf2=zf/(1.-zf)

            endif

         enddo

      enddo

            enddo

         enddo

      endif

      end subroutine thermcell_main

!=============================================================================

!/////////////////////////////////////////////////////////////////////////////

!=============================================================================

    &            zqla,f_star,zw2,comment)                          ! in

!=============================================================================

      USE lmdz_thermcell_ini, ONLY: prt_level

      IMPLICIT NONE

      integer i, k, ngrid,nlay

      real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,po,ztva,zqla

      real, intent(in), dimension(ngrid,nlay) :: f_star,zw2

      integer, intent(in), dimension(ngrid) :: long

      real seuil

      character*21 comment

      seuil=0.25

      endif

      return

!  test sur la hauteur des thermiques ...

         do i=1,ngrid

!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then

           if (prt_level.ge.10) then

               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)

               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'

               do k=1,nlay

                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)

               enddo

           endif

         enddo

      return

      end

! nrlmd le 10/04/2012   Transport de la TKE par le thermique moyen pour la fermeture en ALP

!                       On transporte pbl_tke pour donner therm_tke

!                       Copie conforme de la subroutine DTKE dans physiq.F ecrite par Frederic Hourdin

!=======================================================================

!///////////////////////////////////////////////////////////////////////

!=======================================================================

     &     therm_tke_max)                                ! out

      USE lmdz_thermcell_ini, ONLY: prt_level

      implicit none

!=======================================================================

!

!   Calcul du transport verticale dans la couche limite en presence

!   de "thermiques" explicitement representes

!   calcul du dq/dt une fois qu'on connait les ascendances

!

!=======================================================================

      integer ngrid,nlay

      real, intent(in) :: ptimestep

      real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev

      real, intent(in), dimension(ngrid,nlay) :: entr0

      real, intent(in) :: rg

      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max

      integer lev_out                           ! niveau pour les print

      integer ig,k

      lev_out=0

!   calcul du detrainement

      enddo

! Decalage vertical des entrainements et detrainements.

      enddo

!!! nrlmd le 16/09/2010

      enddo

!!!

    if (1==1) then

     &         1.e-5*masse(ig,k)) then

         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &

            else

            endif

            if (qa(ig,k).lt.0.) then

!               print*,'qa<0!!!'

            endif

!               print*,'q<0!!!'

            endif

         enddo

      enddo

! Calcul du flux subsident

!               print*,'wqd<0!!!'

            endif

         enddo

      enddo

      enddo

! Calcul des tendances

            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &

     &               -wqd(ig,k)+wqd(ig,k+1))  &

         enddo

      enddo

 endif

!!! fin nrlmd le 10/04/2012

     end

END MODULE lmdz_thermcell_main