Changeset 4652


Ignore:
Timestamp:
Aug 28, 2023, 4:41:57 PM (15 months ago)
Author:
fhourdin
Message:

Petite correction de la commission precedente

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4651 r4652  
    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,iflag_ratqs
    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=1
    242 if ( iflag_replay == 1 ) CALL iophys_ini(pdtphys)
    243 
    24495
    24596endif ! of if (debut)
    … …  
    249100!------------------------------------------------------------
    250101
    251 d_u_dyn(:,:)=(u(:,:)-u_prev(:,:))/pdtphys
    252 d_v_dyn(:,:)=(v(:,:)-v_prev(:,:))/pdtphys
    253 d_t_dyn(:,:)=(temp(:,:)-t_prev(:,:))/pdtphys
    254 d_q_dyn(:,:,:)=(qx(:,:,:)-q_prev(:,:,:))/pdtphys
    255102
    256103! set all tendencies to zero
    … …  
    261108d_ps(1:klon)=0.
    262109
    263 u_loc(1:klon,1:klev)=u(1:klon,1:klev)
    264 v_loc(1:klon,1:klev)=v(1:klon,1:klev)
    265 t_loc(1:klon,1:klev)=temp(1:klon,1:klev)
    266 d_u_loc(1:klon,1:klev)=0.
    267 d_v_loc(1:klon,1:klev)=0.
    268 d_t_loc(1:klon,1:klev)=0.
    269 do iq=1,nqtot
    270    do k=1,klev
    271       do i=1,klon
    272          q_loc(i,k,iq)=qx(i,k,iq)
    273       enddo
    274    enddo
    275 enddo
    276 
    277 du0(1:klon,1:klev)=0.
    278 dv0(1:klon,1:klev)=0.
    279 dqbs0(1:klon,1:klev)=0.
    280 
    281 
    282 
    283110!------------------------------------------------------------
    284111! Calculs
    285112!------------------------------------------------------------
    286113
    287 !------------------------------------------------------------
    288 ! Rappel en temperature et frottement dans la premiere chouche
    289 !------------------------------------------------------------
    290 
    291 if ( iflag_newton == 1 ) then
    292     ! compute tendencies to return to the dynamics:
    293     ! "friction" on the first layer
    294     d_u(1:klon,1)=-u(1:klon,1)/86400.
    295     d_v(1:klon,1)=-v(1:klon,1)/86400.
    296     ! newtonian relaxation towards temp_newton()
    297     do k=1,klev
    298        temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
    299        d_t(1:klon,k)=(temp_newton(1:klon,k)-temp(1:klon,k))/5.e5
    300     enddo
    301 else
    302    temp_newton(:,:)=0.
    303 endif
    304 
    305 
    306 !------------------------------------------------------------
    307 ! Reevaporation de la pluie
    308 !------------------------------------------------------------
    309 
    310 iflag_ice_thermo=1
    311 if ( iflag_reevap == 1 ) then
    312       CALL reevap (klon,klev,iflag_ice_thermo,t_loc,q_loc(:,:,1),q_loc(:,:,2),q_loc(:,:,3), &
    313 &                  d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
    314      do k=1,klev
    315         do i=1,klon
    316            t_loc(i,k)=t_loc(i,k)+d_t_eva(i,k)
    317            q_loc(i,k,1)=q_loc(i,k,1)+d_q_eva(i,k)
    318            q_loc(i,k,2)=q_loc(i,k,2)+d_ql_eva(i,k)
    319            q_loc(i,k,3)=q_loc(i,k,3)+d_qi_eva(i,k)
    320            q_loc(i,k,2)=0.
    321            q_loc(i,k,3)=0.
    322         enddo
    323      enddo
    324 else
    325      d_t_eva(:,:)=0.
    326      d_q_eva(:,:)=0.
    327      d_ql_eva(:,:)=0.
    328      d_qi_eva(:,:)=0.
    329 endif
    330 
    331 
    332 
    333 !-----------------------------------------------------------------------
    334 ! Variables intermédiaires (altitudes, temperature potentielle ...)
    335 !-----------------------------------------------------------------------
    336 
    337 DO k=1,klev
    338    DO i=1,klon
    339       zzlay(i,k)=pphi(i,k)/rg
    340    ENDDO
    341 ENDDO
    342 DO i=1,klon
    343    zzlev(i,1)=0.
    344 ENDDO
    345 DO k=2,klev
    346    DO i=1,klon
    347       z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
    348       z2=(paprs(i,k)+pplay(i,k))/(paprs(i,k)-pplay(i,k))
    349       zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
    350    ENDDO
    351 ENDDO
    352 
    353 !   Transformation de la temperature en temperature potentielle
    354 DO k=1,klev
    355    DO i=1,klon
    356       zpopsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
    357       masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
    358    ENDDO
    359 ENDDO
    360 DO k=1,klev
    361    DO i=1,klon
    362       h_loc(i,k)=t_loc(i,k)/zpopsk(i,k)
    363       d_h_loc(i,k)=d_t_loc(i,k)/zpopsk(i,k)
    364       d_q_loc(i,k,1)=0.
    365    ENDDO
    366 ENDDO
    367 
    368 !-----------------------------------------------------------------------
    369 ! Diffusion verticale
    370 !-----------------------------------------------------------------------
    371 
    372 if ( iflag_vdif == 1 ) then
    373       emis(:)=1.
    374       !tsrf=300.
    375       z0m=0.01
    376       z0h=0.01
    377       capcal=1e5
    378       lwrite=.false.
    379       print*,'lwrite ',lwrite
    380             call vdif(klon,klev,                                   &
    381      &                pdtphys,capcal,z0m,z0h,                  &
    382      &                pplay,paprs,zzlay,zzlev,                      &
    383      &                u_loc,v_loc,t_loc,h_loc,q_loc,tsrf,emis,                         &
    384      &                d_u_loc,d_v_loc,d_h_loc,d_q_loc,fluxsrf,                   &
    385      &                d_u_vdif,d_v_vdif,d_h_vdif,d_q_vdif,dtsrf,q2,kz_v,kz_h,     &
    386      &                richardson,cdv,cdh,                         &
    387      &                lwrite)
    388       do k=1,klev
    389          do i=1,klon
    390             d_t_vdif(i,k)=d_h_vdif(i,k)*zpopsk(i,k)
    391             t_loc(i,k)=t_loc(i,k)+d_t_vdif(i,k)*pdtphys
    392             u_loc(i,k)=u_loc(i,k)+d_u_vdif(i,k)*pdtphys
    393             v_loc(i,k)=v_loc(i,k)+d_v_vdif(i,k)*pdtphys
    394             q_loc(i,k,1)=q_loc(i,k,1)+d_q_vdif(i,k,1)*pdtphys
    395          enddo
    396       enddo
    397       do i=1,klon
    398          tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys
    399       enddo
    400 else
    401       d_u_vdif(:,:)=0.
    402       d_v_vdif(:,:)=0.
    403       d_t_vdif(:,:)=0.
    404       d_h_vdif(:,:)=0.
    405       d_q_vdif(:,:,1)=0.
    406       kz_v(:,:)=0.
    407       kz_h(:,:)=0.
    408       richardson(:,:)=0.
    409 endif
    410 
    411 !-----------------------------------------------------------------------
    412 ! Thermiques
    413 !-----------------------------------------------------------------------
    414 
     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()
    415119do k=1,klev
    416    do i=1,klon
    417       d_u_the(i,k)=0.
    418       d_v_the(i,k)=0.
    419       d_t_the(i,k)=0.
    420       d_q_the(i,k,1)=0.
    421    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
    422122enddo
    423123
    424 if ( iflag_thermals > 0 ) then
    425 
    426 
    427           zqta(:,:)=q_loc(:,:,1)
    428           call caltherm(pdtphys  &
    429          &      ,pplay,paprs,pphi  &
    430          &      ,u_loc,v_loc,t_loc,q_loc,debut  &
    431          &      ,f0,zmax0,d_u_the,d_v_the,d_t_the,d_q_the  &
    432          &     ,frac_the,fm_the,zqta,ztv,zpspsk,ztla,zthl &
    433          &   )
    434    
    435           do k=1,klev
    436              do i=1,klon
    437                 t_loc(i,k)=t_loc(i,k)+d_t_the(i,k)
    438                 u_loc(i,k)=u_loc(i,k)+d_u_the(i,k)
    439                 v_loc(i,k)=v_loc(i,k)+d_v_the(i,k)
    440                 q_loc(i,k,1)=q_loc(i,k,1)+d_q_the(i,k,1)
    441              enddo
    442           enddo
    443 
    444 else
    445           frac_the(:,:)=0.
    446           fm_the(:,:)=0.
    447           ztv(:,:)=t_loc(:,:)
    448           zqta(:,:)=q_loc(:,:,1)
    449           ztla(:,:)=0.
    450           zthl(:,:)=0.
    451           zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
    452  
    453 endif
    454 
    455 !-----------------------------------------------------------------------
    456 ! Condensation grande échelle
    457 !-----------------------------------------------------------------------
    458 
    459 iflag_cld_th=5
    460 ok_ice_sursat=.false.
    461 ptconv(:,:)=.false.
    462 distcltop=0.
    463 temp_cltop=0.
    464 beta(:,:)=1.
    465 rneb_seri(:,:)=0.
    466 do k=1,klev
    467   ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
    468   *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
    469 enddo
    470 
    471 
    472 if ( iflag_lscp == 1 ) then
    473 
    474     call lscp(klon,klev,pdtphys,missing_val,            &
    475      paprs,pplay,t_loc,q_loc,ptconv,ratqs,                      &
    476      d_t_lscp, d_q_lscp, d_ql_lscp, d_qi_lscp, rneb, rneblsvol, rneb_seri,  &
    477      pfraclr,pfracld,                                   &
    478      radocond, radicefrac, rain, snow,                  &
    479      frac_impa, frac_nucl, beta,                        &
    480      prfl, psfl, rhcl, zqta, frac_the,                     &
    481      ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
    482      iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
    483      distcltop,temp_cltop,                               &
    484      qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
    485      Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
    486      cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    487 
    488 
    489           do k=1,klev
    490              do i=1,klon
    491                 t_loc(i,k)=t_loc(i,k)+d_t_lscp(i,k)
    492                 q_loc(i,k,1)=q_loc(i,k,1)+d_q_lscp(i,k)
    493                 q_loc(i,k,2)=q_loc(i,k,2)+d_ql_lscp(i,k)
    494                 q_loc(i,k,3)=q_loc(i,k,3)+d_qi_lscp(i,k)
    495              enddo
    496           enddo
    497 
    498 else
    499      d_t_lscp(:,:)=0.
    500      d_q_lscp(:,:)=0.
    501      d_ql_lscp(:,:)=0.
    502      d_qi_lscp(:,:)=0.
    503      rneb(:,:)=0.
    504      rneblsvol(:,:)=0.
    505      pfraclr(:,:)=0.
    506      pfracld(:,:)=0.
    507      radocond(:,:)=0.
    508      rain(:)=0.
    509      snow(:)=0.
    510      radicefrac(:,:)=0.
    511      rhcl      (:,:)=0.
    512      qsatl     (:,:)=0.
    513      qsats     (:,:)=0.
    514      prfl      (:,:)=0.
    515      psfl      (:,:)=0.
    516      distcltop (:,:)=0.
    517      temp_cltop(:,:)=0.
    518      frac_impa (:,:)=0.
    519      frac_nucl (:,:)=0.
    520      qclr      (:,:)=0.
    521      qcld      (:,:)=0.
    522      qss       (:,:)=0.
    523      qvc       (:,:)=0.
    524      rnebclr   (:,:)=0.
    525      rnebss    (:,:)=0.
    526      gamma_ss  (:,:)=0.
    527      Tcontr    (:,:)=0.
    528      qcontr    (:,:)=0.
    529      qcontr2   (:,:)=0.
    530      fcontrN   (:,:)=0.
    531      fcontrP   (:,:)=0.
    532 endif
    533 
    534 
    535 d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys
    536 d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys
    537 d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys
    538 d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys
    539124
    540125!------------------------------------------------------------
    … …  
    543128
    544129
    545 tsrf_(:)=tsrf(:)
    546 if ( iflag_replay == 0 ) then
    547    call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,temp,qx,0.*u,0.*u,0.*u,0.*u,q2,0.*u)
    548 else
    549      ! En mode replay, on sort aussi les variables de base
    550 ! Les lignes qui suivent ont été générées automatiquement avec :
    551 ! ( 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
    552 ! ( 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
    553 include "physiqex_out.h"
    554 
    555 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)
    556131
    557132
    … …  
    561136endif
    562137
    563 print*,'Fin physiqex'
    564138
    565139end subroutine physiqex
Note: See TracChangeset for help on using the changeset viewer.