Changeset 4658 for LMDZ6/trunk/libf


Ignore:
Timestamp:
Aug 30, 2023, 6:33:15 PM (15 months ago)
Author:
fhourdin
Message:

Petite correction de la commission precedente.

Deso

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4657 r4658  
    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
    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"
    … …  
    16157    !$OMP THREADPRIVATE(clesphy0)
    16258
    163 real,dimension(klon,klev) :: temp_newton
    164 integer :: i,k,iq
    165 INTEGER, SAVE :: itap=0
    166 !$OMP THREADPRIVATE(itap)
    167 INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
    168 !$OMP THREADPRIVATE(abortphy)
    16959
    170 integer, save :: iflag_reevap=1,iflag_newton=0,iflag_vdif=1,iflag_lscp=1,iflag_cloudth_vert=3,iflag_ratqs=4
    171 !$OMP THREADPRIVATE(iflag_reevap,iflag_newton,iflag_vdif,iflag_lscp,iflag_cloudth_vert,iflag_ratqs)
    172 
     60real :: temp_newton(klon,klev)
     61integer :: k
    17362logical, save :: first=.true.
    17463!$OMP THREADPRIVATE(first)
     64
     65real,save :: rg=9.81
     66!$OMP THREADPRIVATE(rg)
    17567
    17668! For I/Os
    17769integer :: itau0
    17870real :: zjulian
    179 real,dimension(klon,klev) :: du0,dv0,dqbs0
    180 real,dimension(klon,klev) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
    18171
    18272
    … …  
    20393  ! Initialize IOIPSL output file
    20494#endif
    205 call suphel
    206 call vdif_ini_(klon,RCPD,RD,RG,RKAPPA)
    207 ! Pourquoi ce tau_thermals en argument ??? AFAIRE
    208 tau_thermals=0.
    209 call getin_p('iflag_thermals',iflag_thermals)
    210 
    211 call getin_p('iflag_newton',iflag_newton)
    212 call getin_p('iflag_reevap',iflag_reevap)
    213 call getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
    214 call getin_p('iflag_ratqs',iflag_ratqs)
    215 call getin_p('iflag_vdif',iflag_vdif)
    216 call getin_p('iflag_lscp',iflag_lscp)
    217 call getin_p('ratqsbas',ratqsbas)
    218 call getin_p('ratqshaut',ratqshaut)
    219 call getin_p('ratqsp0',ratqsp0)
    220 call getin_p('ratqsdp',ratqsdp)
    221 CALL thermcell_ini(iflag_thermals,0,tau_thermals,6, &
    222    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
    223 CALL lscp_ini(pdtphys,.false.,iflag_ratqs, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
    224 
    225 
    226 
    227 allocate(tsrf(klon),q2(klon,klev+1),f0(klon),zmax0(klon))
    228 allocate(u_prev(klon,klev),v_prev(klon,klev),t_prev(klon,klev),q_prev(klon,klev,nqtot))
    229 
    230 u_prev(:,:)=u(:,:)
    231 v_prev(:,:)=v(:,:)
    232 t_prev(:,:)=temp(:,:)
    233 q_prev(:,:,:)=qx(:,:,:)
    234 
    235 q2=1.e-10
    236 tsrf=temp(:,1)
    237 f0=0.
    238 zmax0=0.
    239 
    240 iflag_replay=0
    241 call getin_p('iflag_replay',iflag_replay)
    242 if ( iflag_replay >= 0 ) 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 ! Variables intermédiaires (altitudes, temperature potentielle ...)
    333 !-----------------------------------------------------------------------
    334 
    335 DO k=1,klev
    336    DO i=1,klon
    337       zzlay(i,k)=pphi(i,k)/rg
    338    ENDDO
    339 ENDDO
    340 DO i=1,klon
    341    zzlev(i,1)=0.
    342 ENDDO
    343 DO k=2,klev
    344    DO i=1,klon
    345       z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
    346       z2=(paprs(i,k)+pplay(i,k))/(paprs(i,k)-pplay(i,k))
    347       zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
    348    ENDDO
    349 ENDDO
    350 
    351 !   Transformation de la temperature en temperature potentielle
    352 DO k=1,klev
    353    DO i=1,klon
    354 !      zpspsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
    355 zpspsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
    356 masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
    357 ENDDO
    358 ENDDO
    359 DO k=1,klev
    360 DO i=1,klon
    361 h_loc(i,k)=t_loc(i,k)/zpspsk(i,k)
    362 d_h_loc(i,k)=d_t_loc(i,k)/zpspsk(i,k)
    363 d_q_loc(i,k,1)=0.
    364 ENDDO
    365 ENDDO
    366 
    367 !-----------------------------------------------------------------------
    368 ! Diffusion verticale
    369 !-----------------------------------------------------------------------
    370 
    371 if ( iflag_vdif == 1 ) then
    372 emis(:)=1.
    373 !tsrf=300.
    374 z0m=0.035
    375 z0h=0.035
    376 capcal=1e2
    377 lwrite=.false.
    378 print*,'lwrite ',lwrite
    379     call vdif(klon,klev,                                               &
    380 &                pdtphys,capcal,z0m,z0h,                                  &
    381 &                pplay,paprs,zzlay,zzlev,                                 &
    382 &                u_loc,v_loc,t_loc,h_loc,q_loc,tsrf,emis,                 &
    383 &                d_u_loc,d_v_loc,d_h_loc,d_q_loc,fluxsrf,                 &
    384 &                d_u_vdif,d_v_vdif,d_h_vdif,d_q_vdif,dtsrf,q2,kz_v,kz_h,  &
    385 &                richardson,cdv,cdh,                                      &
    386 &                lwrite)
     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()
    387119do k=1,klev
    388  do i=1,klon
    389     d_t_vdif(i,k)=d_h_vdif(i,k)*zpspsk(i,k)
    390     t_loc(i,k)=t_loc(i,k)+d_t_vdif(i,k)*pdtphys
    391     u_loc(i,k)=u_loc(i,k)+d_u_vdif(i,k)*pdtphys
    392     v_loc(i,k)=v_loc(i,k)+d_v_vdif(i,k)*pdtphys
    393     q_loc(i,k,1)=q_loc(i,k,1)+d_q_vdif(i,k,1)*pdtphys
    394  enddo
    395 enddo
    396 do i=1,klon
    397  tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys
    398 enddo
    399 else
    400 d_u_vdif(:,:)=0.
    401 d_v_vdif(:,:)=0.
    402 d_t_vdif(:,:)=0.
    403 d_h_vdif(:,:)=0.
    404 d_q_vdif(:,:,1)=0.
    405 kz_v(:,:)=0.
    406 kz_h(:,:)=0.
    407 richardson(:,:)=0.
    408 endif
    409 
    410 !-----------------------------------------------------------------------
    411 ! Thermiques
    412 !-----------------------------------------------------------------------
    413 
    414 do k=1,klev
    415 do i=1,klon
    416 d_u_the(i,k)=0.
    417 d_v_the(i,k)=0.
    418 d_t_the(i,k)=0.
    419 d_q_the(i,k,1)=0.
    420 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
    421122enddo
    422123
    423 if ( iflag_thermals > 0 ) then
    424 
    425 
    426   zqta(:,:)=q_loc(:,:,1)
    427   call caltherm(pdtphys                            &
    428  &      ,pplay,paprs,pphi                          &
    429  &      ,u_loc,v_loc,t_loc,q_loc,debut             &
    430  &      ,f0,zmax0,d_u_the,d_v_the,d_t_the,d_q_the  &
    431  &     ,frac_the,fm_the,zqta,ztv,zpspsk,ztla,zthl  &
    432  &   )
    433 
    434   do k=1,klev
    435      do i=1,klon
    436         t_loc(i,k)=t_loc(i,k)+d_t_the(i,k)
    437         u_loc(i,k)=u_loc(i,k)+d_u_the(i,k)
    438         v_loc(i,k)=v_loc(i,k)+d_v_the(i,k)
    439         q_loc(i,k,1)=q_loc(i,k,1)+d_q_the(i,k,1)
    440      enddo
    441   enddo
    442 
    443 else
    444   frac_the(:,:)=0.
    445   fm_the(:,:)=0.
    446   ztv(:,:)=t_loc(:,:)
    447   zqta(:,:)=q_loc(:,:,1)
    448   ztla(:,:)=0.
    449   zthl(:,:)=0.
    450 endif
    451 
    452 !-----------------------------------------------------------------------
    453 ! Condensation grande échelle
    454 !-----------------------------------------------------------------------
    455 
    456 iflag_cld_th=5
    457 ok_ice_sursat=.false.
    458 ptconv(:,:)=.false.
    459 distcltop=0.
    460 temp_cltop=0.
    461 beta(:,:)=1.
    462 rneb_seri(:,:)=0.
    463 do k=1,klev
    464   ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
    465   *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
    466 enddo
    467 
    468 
    469 if ( iflag_lscp == 1 ) then
    470 
    471     call lscp(klon,klev,pdtphys,missing_val,            &
    472      paprs,pplay,t_loc,q_loc,ptconv,ratqs,                      &
    473      d_t_lscp, d_q_lscp, d_ql_lscp, d_qi_lscp, rneb, rneblsvol, rneb_seri,  &
    474      pfraclr,pfracld,                                   &
    475      radocond, radicefrac, rain, snow,                  &
    476      frac_impa, frac_nucl, beta,                        &
    477      prfl, psfl, rhcl, zqta, frac_the,                     &
    478      ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
    479      iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
    480      distcltop,temp_cltop,                               &
    481      qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
    482      Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
    483      cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    484 
    485 
    486           do k=1,klev
    487              do i=1,klon
    488                 t_loc(i,k)=t_loc(i,k)+d_t_lscp(i,k)
    489                 q_loc(i,k,1)=q_loc(i,k,1)+d_q_lscp(i,k)
    490                 q_loc(i,k,2)=q_loc(i,k,2)+d_ql_lscp(i,k)
    491                 q_loc(i,k,3)=q_loc(i,k,3)+d_qi_lscp(i,k)
    492              enddo
    493           enddo
    494 
    495 else
    496      d_t_lscp(:,:)=0.
    497      d_q_lscp(:,:)=0.
    498      d_ql_lscp(:,:)=0.
    499      d_qi_lscp(:,:)=0.
    500      rneb(:,:)=0.
    501      rneblsvol(:,:)=0.
    502      pfraclr(:,:)=0.
    503      pfracld(:,:)=0.
    504      radocond(:,:)=0.
    505      rain(:)=0.
    506      snow(:)=0.
    507      radicefrac(:,:)=0.
    508      rhcl      (:,:)=0.
    509      qsatl     (:,:)=0.
    510      qsats     (:,:)=0.
    511      prfl      (:,:)=0.
    512      psfl      (:,:)=0.
    513      distcltop (:,:)=0.
    514      temp_cltop(:,:)=0.
    515      frac_impa (:,:)=0.
    516      frac_nucl (:,:)=0.
    517      qclr      (:,:)=0.
    518      qcld      (:,:)=0.
    519      qss       (:,:)=0.
    520      qvc       (:,:)=0.
    521      rnebclr   (:,:)=0.
    522      rnebss    (:,:)=0.
    523      gamma_ss  (:,:)=0.
    524      Tcontr    (:,:)=0.
    525      qcontr    (:,:)=0.
    526      qcontr2   (:,:)=0.
    527      fcontrN   (:,:)=0.
    528      fcontrP   (:,:)=0.
    529 endif
    530 
    531 
    532 d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys
    533 d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys
    534 d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys
    535 d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys
    536124
    537125!------------------------------------------------------------
    … …  
    540128
    541129
    542 tsrf_(:)=tsrf(:)
    543 if ( iflag_replay == -1 ) then
    544    call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,temp,qx,0.*u,0.*u,0.*u,0.*u,q2,0.*u)
    545 else if (iflag_replay == 0 ) then
    546      ! En mode replay, on sort aussi les variables de base
    547 ! Les lignes qui suivent ont été générées automatiquement avec :
    548 ! ( 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
    549 ! ( 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
    550 include "physiqex_out.h"
    551 
    552 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)
    553131
    554132
    … …  
    558136endif
    559137
    560 print*,'Fin physiqex'
    561138
    562139end subroutine physiqex
Note: See TracChangeset for help on using the changeset viewer.