Changeset 4655


Ignore:
Timestamp:
Aug 30, 2023, 10:12:42 AM (9 months ago)
Author:
fhourdin
Message:

Retour en arriere sur une commission intempestive

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/physiqex_mod.F90

    r4654 r4655  
    99     &            debut,lafin,pdtphys, &
    1010     &            paprs,pplay,pphi,pphis,presnivs, &
    11      &            u,v,rot,temp,qx, &
     11     &            u,v,rot,t,qx, &
    1212     &            flxmass_w, &
    1313     &            d_u, d_v, d_t, d_qx, d_ps)
    1414
    15 
    16 USE dimphy, only : klon,klev
    17 USE infotrac_phy, only : nqtot
    18 USE geometry_mod, only : latitude
    19 USE ioipsl, only : ymds2ju
    20 USE phys_state_var_mod, only : phys_state_var_init
    21 USE phyetat0_mod, only: phyetat0
    22 USE output_physiqex_mod, ONLY: output_physiqex
    23 use vdif_ini, only : vdif_ini_
    24 USE lmdz_thermcell_ini, ONLY : thermcell_ini
    25 USE ioipsl_getin_p_mod, ONLY : getin_p
    26 USE wxios, ONLY: missing_val, using_xios
    27 USE lscp_mod, ONLY : lscp
    28 USE lscp_ini_mod, ONLY : lscp_ini
    29 USE add_phys_tend_mod, ONLY : add_phys_tend
    30 
     15      USE dimphy, only : klon,klev
     16      USE infotrac_phy, only : nqtot
     17      USE geometry_mod, only : latitude
     18!      USE comcstphy, only : rg
     19      USE ioipsl, only : ymds2ju
     20      USE phys_state_var_mod, only : phys_state_var_init
     21      USE phyetat0_mod, only: phyetat0
     22      USE output_physiqex_mod, ONLY: output_physiqex
    3123
    3224      IMPLICIT none
    33 
    34 include "YOETHF.h"
    35 
    36 
    37 
    38 
    3925!
    4026! Routine argument:
    … …  
    4632      logical,intent(in) :: lafin ! signals last call to physics
    4733      real,intent(in) :: pdtphys ! physics time step (s)
    48       real,dimension(klon,klev+1),intent(in) :: paprs ! interlayer pressure (Pa)
    49       real,dimension(klon,klev),intent(in) :: pplay ! mid-layer pressure (Pa)
    50       real,dimension(klon,klev),intent(in) :: pphi ! geopotential at mid-layer
    51       real,dimension(klon),intent(in) :: pphis ! surface geopotential
    52       real,dimension(klev),intent(in) :: presnivs ! pseudo-pressure (Pa) of mid-layers
    53       real,dimension(klon,klev),intent(in) :: u ! eastward zonal wind (m/s)
    54       real,dimension(klon,klev),intent(in) :: v ! northward meridional wind (m/s)
    55       real,dimension(klon,klev),intent(in) :: rot ! northward meridional wind (m/s)
    56       real,dimension(klon,klev),intent(in) :: temp ! temperature (K)
    57       real,dimension(klon,klev,nqtot),intent(in) :: qx  ! tracers (.../kg_air)
    58       real,dimension(klon,klev),intent(in) :: flxmass_w ! vertical mass flux
    59       real,dimension(klon,klev),intent(out) :: d_u ! physics tendency on u (m/s/s)
    60       real,dimension(klon,klev),intent(out) :: d_v ! physics tendency on v (m/s/s)
    61       real,dimension(klon,klev),intent(out) :: d_t ! physics tendency on t (K/s)
    62       real,dimension(klon,klev,nqtot),intent(out) :: d_qx ! physics tendency on tracers
    63       real,dimension(klon),intent(out) :: d_ps ! physics tendency on surface pressure
    64 
    65       real, dimension(klon,klev) :: u_loc
    66       real, dimension(klon,klev) :: v_loc
    67       real, dimension(klon,klev) :: t_loc
    68       real, dimension(klon,klev) :: h_loc
    69       real, dimension(klon,klev) :: d_u_loc,d_v_loc,d_t_loc,d_h_loc
    70 
    71       real, dimension(klon,klev) :: d_u_dyn,d_v_dyn,d_t_dyn
    72       real, dimension(klon,klev,nqtot) :: d_q_dyn
    73       real, allocatable, dimension(:,:), save :: u_prev,v_prev,t_prev
    74       real, allocatable, dimension(:,:,:), save :: q_prev
    75 !$OMP THREADPRIVATE(u_prev,v_prev,t_prev,q_prev)
    76 
    77 
    78 
    79       real, dimension(klon,klev) :: d_u_vdif,d_v_vdif,d_t_vdif,d_h_vdif
    80       real, dimension(klon,klev) :: d_u_the,d_v_the,d_t_the
    81       real, dimension(klon,klev,nqtot) :: q_loc,d_q_loc,d_q_vdif,d_q_the
    82 
    83 real, dimension(klon) :: capcal,z0m,z0h,dtsrf,emis,fluxsrf,cdh,cdv,tsrf_
    84 real, dimension(klon,klev) :: zzlay,masse,zpopsk
    85 real, dimension(klon,klev+1) :: zzlev,kz_v,kz_h,richardson
    86 
    87 real, save, allocatable, dimension(:) :: tsrf,f0,zmax0
    88 real, save, allocatable, dimension(:,:) :: q2
    89 !$OMP THREADPRIVATE(tsrf,q2,f0,zmax0)
    90 
    91     real,save :: ratqsbas=0.002,ratqshaut=0.3,ratqsp0=50000.,ratqsdp=20000.
    92     !$OMP THREADPRIVATE(ratqsbas,ratqshaut,ratqsp0,ratqsdp)
    93 
    94 
    95 real :: z1,z2,tau_thermals
    96 logical :: lwrite
    97 integer :: iflag_replay
    98 
    99 integer :: iflag_thermals=18
    100 
    101 !--------------------------------------------------------------
    102 ! Declaration lscp
    103 !--------------------------------------------------------------
    104   INTEGER                         :: iflag_cld_th    ! flag that determines the distribution of convective clouds    ! IN
    105   INTEGER                         :: iflag_ice_thermo! flag to activate the ice thermodynamics                       ! IN
    106   LOGICAL                         :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated           ! IN
    107   LOGICAL, DIMENSION(klon,klev)   :: ptconv          ! grid points where deep convection scheme is active            ! IN
    108   REAL, DIMENSION(klon,klev)      :: ztv             ! virtual potential temperature [K]                             ! IN
    109   REAL, DIMENSION(klon,klev)      :: zqta            ! specific humidity within thermals [kg/kg]                     ! IN
    110   REAL, DIMENSION(klon,klev+1)      :: frac_the,fm_the
    111   REAL, DIMENSION(klon,klev)      :: zpspsk          ! exner potential (p/100000)**(R/cp)                            ! IN
    112   REAL, DIMENSION(klon,klev)      :: ztla            ! liquid temperature within thermals [K]                        ! IN
    113   REAL, DIMENSION(klon,klev)         :: zthl         ! liquid potential temperature [K]                              ! INOUT
    114   REAL, DIMENSION(klon,klev)      :: ratqs            ! function of pressure that sets the large-scale               ! INOUT
    115   REAL, DIMENSION(klon,klev)      :: beta             ! conversion rate of condensed water                           ! INOUT
    116   REAL, DIMENSION(klon,klev)      :: rneb_seri        ! fraction nuageuse en memoire                                 ! INOUT
    117   REAL, DIMENSION(klon,klev)      :: d_t_lscp         ! temperature increment [K]                                    ! OUT
    118   REAL, DIMENSION(klon,klev)      :: d_q_lscp         ! specific humidity increment [kg/kg]                          ! OUT
    119   REAL, DIMENSION(klon,klev)      :: d_ql_lscp        ! liquid water increment [kg/kg]                               ! OUT
    120   REAL, DIMENSION(klon,klev)      :: d_qi_lscp        ! cloud ice mass increment [kg/kg]                             ! OUT
    121   REAL, DIMENSION(klon,klev)      :: rneb             ! cloud fraction [-]                                           ! OUT
    122   REAL, DIMENSION(klon,klev)      :: rneblsvol        ! cloud fraction per unit volume [-]                           ! OUT
    123   REAL, DIMENSION(klon,klev)      :: pfraclr          ! precip fraction clear-sky part [-]                           ! OUT
    124   REAL, DIMENSION(klon,klev)      :: pfracld          ! precip fraction cloudy part [-]                              ! OUT
    125   REAL, DIMENSION(klon,klev)      :: radocond         ! condensed water used in the radiation scheme [kg/kg]         ! OUT
    126   REAL, DIMENSION(klon,klev)      :: radicefrac       ! ice fraction of condensed water for radiation scheme         ! OUT
    127   REAL, DIMENSION(klon,klev)      :: rhcl             ! clear-sky relative humidity [-]                              ! OUT
    128   REAL, DIMENSION(klon)           :: rain             ! surface large-scale rainfall [kg/s/m2]                       ! OUT
    129   REAL, DIMENSION(klon)           :: snow             ! surface large-scale snowfall [kg/s/m2]                       ! OUT
    130   REAL, DIMENSION(klon,klev)      :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]              ! OUT
    131   REAL, DIMENSION(klon,klev)      :: qsats            ! saturation specific humidity wrt ice [kg/kg]                 ! OUT
    132   REAL, DIMENSION(klon,klev+1)    :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]            ! OUT
    133   REAL, DIMENSION(klon,klev+1)    :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]            ! OUT
    134   REAL, DIMENSION(klon,klev)      :: distcltop        ! distance to cloud top [m]                                    ! OUT
    135   REAL, DIMENSION(klon,klev)      :: temp_cltop       ! temperature of cloud top [K]                                 ! OUT
    136   REAL, DIMENSION(klon,klev)      :: frac_impa        ! scavenging fraction due tu impaction [-]                     ! OUT
    137   REAL, DIMENSION(klon,klev)      :: frac_nucl        ! scavenging fraction due tu nucleation [-]                    ! OUT
    138   REAL, DIMENSION(klon,klev)      :: qclr             ! specific total water content in clear sky region [kg/kg]     ! OUT
    139   REAL, DIMENSION(klon,klev)      :: qcld             ! specific total water content in cloudy region [kg/kg]        ! OUT
    140   REAL, DIMENSION(klon,klev)      :: qss              ! specific total water content in supersat region [kg/kg]      ! OUT
    141   REAL, DIMENSION(klon,klev)      :: qvc              ! specific vapor content in clouds [kg/kg]                     ! OUT
    142   REAL, DIMENSION(klon,klev)      :: rnebclr          ! mesh fraction of clear sky [-]                               ! OUT
    143   REAL, DIMENSION(klon,klev)      :: rnebss           ! mesh fraction of ISSR [-]                                    ! OUT
    144   REAL, DIMENSION(klon,klev)      :: gamma_ss         ! coefficient governing the ice nucleation RHi threshold [-]   ! OUT   
    145   REAL, DIMENSION(klon,klev)      :: Tcontr           ! threshold temperature for contrail formation [K]             ! OUT
    146   REAL, DIMENSION(klon,klev)      :: qcontr           ! threshold humidity for contrail formation [kg/kg]            ! OUT
    147   REAL, DIMENSION(klon,klev)      :: qcontr2          ! // (2nd expression more consistent with LMDZ expression of q)! OUT         
    148   REAL, DIMENSION(klon,klev)      :: fcontrN          ! fraction of grid favourable to non-persistent contrails      ! OUT
    149   REAL, DIMENSION(klon,klev)      :: fcontrP          ! fraction of grid favourable to persistent contrails          ! OUT
    150 !--------------------------------------------------------------
    151 
    152   REAL, DIMENSION(klon,klev)      :: d_t_eva,d_q_eva,d_ql_eva,d_qi_eva
    153 include "YOMCST.h"
     34      real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
     35      real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
     36      real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
     37      real,intent(in) :: pphis(klon) ! surface geopotential
     38      real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
     39      real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
     40      real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
     41      real,intent(in) :: rot(klon,klev) ! northward meridional wind (m/s)
     42      real,intent(in) :: t(klon,klev) ! temperature (K)
     43      real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
     44      real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
     45      real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
     46      real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
     47      real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
     48      real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
     49      real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
    15450
    15551!    include "clesphys.h"
    … …  
    16258
    16359
    164 real,dimension(klon,klev) :: temp_newton
    165 integer :: i,k,iq
    166 INTEGER, SAVE :: itap=0
    167 !$OMP THREADPRIVATE(itap)
    168 INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
    169 !$OMP THREADPRIVATE(abortphy)
    170 
    171 integer, save :: iflag_reevap=1,iflag_newton=0,iflag_vdif=1,iflag_lscp=1,iflag_cloudth_vert=3,iflag_ratqs=4
    172 !$OMP THREADPRIVATE(iflag_reevap,iflag_newton,iflag_vdif,iflag_lscp,iflag_cloudth_vert,iflag_ratqs)
    173 
     60real :: temp_newton(klon,klev)
     61integer :: k
    17462logical, save :: first=.true.
    17563!$OMP THREADPRIVATE(first)
     64
     65real,save :: rg=9.81
     66!$OMP THREADPRIVATE(rg)
    17667
    17768! For I/Os
    17869integer :: itau0
    17970real :: zjulian
    180 real,dimension(klon,klev) :: du0,dv0,dqbs0
    181 real,dimension(klon,klev) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
    18271
    18372
    … …  
    20493  ! Initialize IOIPSL output file
    20594#endif
    206 call suphel
    207 call vdif_ini_(klon,RCPD,RD,RG,RKAPPA)
    208 ! Pourquoi ce tau_thermals en argument ??? AFAIRE
    209 tau_thermals=0.
    210 call getin_p('iflag_thermals',iflag_thermals)
    211 
    212 call getin_p('iflag_newton',iflag_newton)
    213 call getin_p('iflag_reevap',iflag_reevap)
    214 call getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
    215 call getin_p('iflag_ratqs',iflag_ratqs)
    216 call getin_p('iflag_vdif',iflag_vdif)
    217 call getin_p('iflag_lscp',iflag_lscp)
    218 call getin_p('ratqsbas',ratqsbas)
    219 call getin_p('ratqshaut',ratqshaut)
    220 call getin_p('ratqsp0',ratqsp0)
    221 call getin_p('ratqsdp',ratqsdp)
    222 CALL thermcell_ini(iflag_thermals,0,tau_thermals,6, &
    223    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
    224 CALL lscp_ini(pdtphys,.false.,iflag_ratqs, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
    225 
    226 
    227 
    228 allocate(tsrf(klon),q2(klon,klev+1),f0(klon),zmax0(klon))
    229 allocate(u_prev(klon,klev),v_prev(klon,klev),t_prev(klon,klev),q_prev(klon,klev,nqtot))
    230 
    231 u_prev(:,:)=u(:,:)
    232 v_prev(:,:)=v(:,:)
    233 t_prev(:,:)=temp(:,:)
    234 q_prev(:,:,:)=qx(:,:,:)
    235 
    236 q2=1.e-10
    237 tsrf=temp(:,1)
    238 f0=0.
    239 zmax0=0.
    240 
    241 iflag_replay=0
    242 call getin_p('iflag_replay',iflag_replay)
    243 if ( iflag_replay >= 0 ) CALL iophys_ini(pdtphys)
    244 
    24595
    24696endif ! of if (debut)
    … …  
    250100!------------------------------------------------------------
    251101
    252 d_u_dyn(:,:)=(u(:,:)-u_prev(:,:))/pdtphys
    253 d_v_dyn(:,:)=(v(:,:)-v_prev(:,:))/pdtphys
    254 d_t_dyn(:,:)=(temp(:,:)-t_prev(:,:))/pdtphys
    255 d_q_dyn(:,:,:)=(qx(:,:,:)-q_prev(:,:,:))/pdtphys
    256102
    257103! set all tendencies to zero
    … …  
    262108d_ps(1:klon)=0.
    263109
    264 u_loc(1:klon,1:klev)=u(1:klon,1:klev)
    265 v_loc(1:klon,1:klev)=v(1:klon,1:klev)
    266 t_loc(1:klon,1:klev)=temp(1:klon,1:klev)
    267 d_u_loc(1:klon,1:klev)=0.
    268 d_v_loc(1:klon,1:klev)=0.
    269 d_t_loc(1:klon,1:klev)=0.
    270 do iq=1,nqtot
    271    do k=1,klev
    272       do i=1,klon
    273          q_loc(i,k,iq)=qx(i,k,iq)
    274       enddo
    275    enddo
    276 enddo
    277 
    278 du0(1:klon,1:klev)=0.
    279 dv0(1:klon,1:klev)=0.
    280 dqbs0(1:klon,1:klev)=0.
    281 
    282 
    283 
    284110!------------------------------------------------------------
    285111! Calculs
    286112!------------------------------------------------------------
    287113
    288 !------------------------------------------------------------
    289 ! Rappel en temperature et frottement dans la premiere chouche
    290 !------------------------------------------------------------
    291 
    292 if ( iflag_newton == 1 ) then
    293     ! compute tendencies to return to the dynamics:
    294     ! "friction" on the first layer
    295     d_u(1:klon,1)=-u(1:klon,1)/86400.
    296     d_v(1:klon,1)=-v(1:klon,1)/86400.
    297     ! newtonian relaxation towards temp_newton()
    298     do k=1,klev
    299        temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
    300        d_t(1:klon,k)=(temp_newton(1:klon,k)-temp(1:klon,k))/5.e5
    301     enddo
    302 else
    303    temp_newton(:,:)=0.
    304 endif
    305 
    306 
    307 !------------------------------------------------------------
    308 ! Reevaporation de la pluie
    309 !------------------------------------------------------------
    310 
    311 iflag_ice_thermo=1
    312 if ( iflag_reevap == 1 ) then
    313       CALL reevap (klon,klev,iflag_ice_thermo,t_loc,q_loc(:,:,1),q_loc(:,:,2),q_loc(:,:,3), &
    314 &                  d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
    315      do k=1,klev
    316         do i=1,klon
    317            t_loc(i,k)=t_loc(i,k)+d_t_eva(i,k)
    318            q_loc(i,k,1)=q_loc(i,k,1)+d_q_eva(i,k)
    319            q_loc(i,k,2)=q_loc(i,k,2)+d_ql_eva(i,k)
    320            q_loc(i,k,3)=q_loc(i,k,3)+d_qi_eva(i,k)
    321            q_loc(i,k,2)=0.
    322            q_loc(i,k,3)=0.
    323         enddo
    324      enddo
    325 else
    326      d_t_eva(:,:)=0.
    327      d_q_eva(:,:)=0.
    328      d_ql_eva(:,:)=0.
    329      d_qi_eva(:,:)=0.
    330 endif
    331 
    332 
    333 
    334 !-----------------------------------------------------------------------
    335 ! Variables intermédiaires (altitudes, temperature potentielle ...)
    336 !-----------------------------------------------------------------------
    337 
    338 DO k=1,klev
    339    DO i=1,klon
    340       zzlay(i,k)=pphi(i,k)/rg
    341    ENDDO
    342 ENDDO
    343 DO i=1,klon
    344    zzlev(i,1)=0.
    345 ENDDO
    346 DO k=2,klev
    347    DO i=1,klon
    348       z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
    349       z2=(paprs(i,k)+pplay(i,k))/(paprs(i,k)-pplay(i,k))
    350       zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
    351    ENDDO
    352 ENDDO
    353 
    354 !   Transformation de la temperature en temperature potentielle
    355 DO k=1,klev
    356    DO i=1,klon
    357       zpopsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
    358       masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
    359    ENDDO
    360 ENDDO
    361 DO k=1,klev
    362    DO i=1,klon
    363       h_loc(i,k)=t_loc(i,k)/zpopsk(i,k)
    364       d_h_loc(i,k)=d_t_loc(i,k)/zpopsk(i,k)
    365       d_q_loc(i,k,1)=0.
    366    ENDDO
    367 ENDDO
    368 
    369 !-----------------------------------------------------------------------
    370 ! Diffusion verticale
    371 !-----------------------------------------------------------------------
    372 
    373 if ( iflag_vdif == 1 ) then
    374       emis(:)=1.
    375       !tsrf=300.
    376       z0m=0.035
    377       z0h=0.035
    378       capcal=1e2
    379       lwrite=.false.
    380       print*,'lwrite ',lwrite
    381             call vdif(klon,klev,                                               &
    382      &                pdtphys,capcal,z0m,z0h,                                  &
    383      &                pplay,paprs,zzlay,zzlev,                                 &
    384      &                u_loc,v_loc,t_loc,h_loc,q_loc,tsrf,emis,                 &
    385      &                d_u_loc,d_v_loc,d_h_loc,d_q_loc,fluxsrf,                 &
    386      &                d_u_vdif,d_v_vdif,d_h_vdif,d_q_vdif,dtsrf,q2,kz_v,kz_h,  &
    387      &                richardson,cdv,cdh,                                      &
    388      &                lwrite)
    389       do k=1,klev
    390          do i=1,klon
    391             d_t_vdif(i,k)=d_h_vdif(i,k)*zpopsk(i,k)
    392             t_loc(i,k)=t_loc(i,k)+d_t_vdif(i,k)*pdtphys
    393             u_loc(i,k)=u_loc(i,k)+d_u_vdif(i,k)*pdtphys
    394             v_loc(i,k)=v_loc(i,k)+d_v_vdif(i,k)*pdtphys
    395             q_loc(i,k,1)=q_loc(i,k,1)+d_q_vdif(i,k,1)*pdtphys
    396          enddo
    397       enddo
    398       do i=1,klon
    399          tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys
    400       enddo
    401 else
    402       d_u_vdif(:,:)=0.
    403       d_v_vdif(:,:)=0.
    404       d_t_vdif(:,:)=0.
    405       d_h_vdif(:,:)=0.
    406       d_q_vdif(:,:,1)=0.
    407       kz_v(:,:)=0.
    408       kz_h(:,:)=0.
    409       richardson(:,:)=0.
    410 endif
    411 
    412 !-----------------------------------------------------------------------
    413 ! Thermiques
    414 !-----------------------------------------------------------------------
    415 
     114! compute tendencies to return to the dynamics:
     115! "friction" on the first layer
     116d_u(1:klon,1)=-u(1:klon,1)/86400.
     117d_v(1:klon,1)=-v(1:klon,1)/86400.
     118! newtonian relaxation towards temp_newton()
    416119do k=1,klev
    417    do i=1,klon
    418       d_u_the(i,k)=0.
    419       d_v_the(i,k)=0.
    420       d_t_the(i,k)=0.
    421       d_q_the(i,k,1)=0.
    422    enddo
     120  temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
     121  d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
    423122enddo
    424123
    425 if ( iflag_thermals > 0 ) then
    426 
    427 
    428           zqta(:,:)=q_loc(:,:,1)
    429           call caltherm(pdtphys                            &
    430          &      ,pplay,paprs,pphi                          &
    431          &      ,u_loc,v_loc,t_loc,q_loc,debut             &
    432          &      ,f0,zmax0,d_u_the,d_v_the,d_t_the,d_q_the  &
    433          &     ,frac_the,fm_the,zqta,ztv,zpspsk,ztla,zthl  &
    434          &   )
    435    
    436           do k=1,klev
    437              do i=1,klon
    438                 t_loc(i,k)=t_loc(i,k)+d_t_the(i,k)
    439                 u_loc(i,k)=u_loc(i,k)+d_u_the(i,k)
    440                 v_loc(i,k)=v_loc(i,k)+d_v_the(i,k)
    441                 q_loc(i,k,1)=q_loc(i,k,1)+d_q_the(i,k,1)
    442              enddo
    443           enddo
    444 
    445 else
    446           frac_the(:,:)=0.
    447           fm_the(:,:)=0.
    448           ztv(:,:)=t_loc(:,:)
    449           zqta(:,:)=q_loc(:,:,1)
    450           ztla(:,:)=0.
    451           zthl(:,:)=0.
    452           zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
    453  
    454 endif
    455 
    456 !-----------------------------------------------------------------------
    457 ! Condensation grande échelle
    458 !-----------------------------------------------------------------------
    459 
    460 iflag_cld_th=5
    461 ok_ice_sursat=.false.
    462 ptconv(:,:)=.false.
    463 distcltop=0.
    464 temp_cltop=0.
    465 beta(:,:)=1.
    466 rneb_seri(:,:)=0.
    467 do k=1,klev
    468   ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
    469   *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
    470 enddo
    471 
    472 
    473 if ( iflag_lscp == 1 ) then
    474 
    475     call lscp(klon,klev,pdtphys,missing_val,            &
    476      paprs,pplay,t_loc,q_loc,ptconv,ratqs,                      &
    477      d_t_lscp, d_q_lscp, d_ql_lscp, d_qi_lscp, rneb, rneblsvol, rneb_seri,  &
    478      pfraclr,pfracld,                                   &
    479      radocond, radicefrac, rain, snow,                  &
    480      frac_impa, frac_nucl, beta,                        &
    481      prfl, psfl, rhcl, zqta, frac_the,                     &
    482      ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
    483      iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
    484      distcltop,temp_cltop,                               &
    485      qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
    486      Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
    487      cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    488 
    489 
    490           do k=1,klev
    491              do i=1,klon
    492                 t_loc(i,k)=t_loc(i,k)+d_t_lscp(i,k)
    493                 q_loc(i,k,1)=q_loc(i,k,1)+d_q_lscp(i,k)
    494                 q_loc(i,k,2)=q_loc(i,k,2)+d_ql_lscp(i,k)
    495                 q_loc(i,k,3)=q_loc(i,k,3)+d_qi_lscp(i,k)
    496              enddo
    497           enddo
    498 
    499 else
    500      d_t_lscp(:,:)=0.
    501      d_q_lscp(:,:)=0.
    502      d_ql_lscp(:,:)=0.
    503      d_qi_lscp(:,:)=0.
    504      rneb(:,:)=0.
    505      rneblsvol(:,:)=0.
    506      pfraclr(:,:)=0.
    507      pfracld(:,:)=0.
    508      radocond(:,:)=0.
    509      rain(:)=0.
    510      snow(:)=0.
    511      radicefrac(:,:)=0.
    512      rhcl      (:,:)=0.
    513      qsatl     (:,:)=0.
    514      qsats     (:,:)=0.
    515      prfl      (:,:)=0.
    516      psfl      (:,:)=0.
    517      distcltop (:,:)=0.
    518      temp_cltop(:,:)=0.
    519      frac_impa (:,:)=0.
    520      frac_nucl (:,:)=0.
    521      qclr      (:,:)=0.
    522      qcld      (:,:)=0.
    523      qss       (:,:)=0.
    524      qvc       (:,:)=0.
    525      rnebclr   (:,:)=0.
    526      rnebss    (:,:)=0.
    527      gamma_ss  (:,:)=0.
    528      Tcontr    (:,:)=0.
    529      qcontr    (:,:)=0.
    530      qcontr2   (:,:)=0.
    531      fcontrN   (:,:)=0.
    532      fcontrP   (:,:)=0.
    533 endif
    534 
    535 
    536 d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys
    537 d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys
    538 d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys
    539 d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys
    540124
    541125!------------------------------------------------------------
    … …  
    544128
    545129
    546 tsrf_(:)=tsrf(:)
    547 if ( iflag_replay == -1 ) then
    548    call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,temp,qx,0.*u,0.*u,0.*u,0.*u,q2,0.*u)
    549 else if (iflag_replay == 0 ) then
    550      ! En mode replay, on sort aussi les variables de base
    551 ! Les lignes qui suivent ont été générées automatiquement avec :
    552 ! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon,klev' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",klev,"","",'$i')' ; done ) > physiqex_out.h
    553 ! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon)' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",1,"","",'$i')' ; done ) >> physiqex_out.h
    554 include "physiqex_out.h"
    555 
    556 endif
     130call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t)
    557131
    558132
    … …  
    562136endif
    563137
    564 print*,'Fin physiqex'
    565138
    566139end subroutine physiqex
Note: See TracChangeset for help on using the changeset viewer.