Changeset 4657


Ignore:
Timestamp:
Aug 30, 2023, 6:24:49 PM (16 months ago)
Author:
fhourdin
Message:

Poursuite nettoyage replauisation

Location:
LMDZ6/trunk/libf/phylmd
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r4590 r4657  
    7878      real fmc_therm(klon,klev+1),zqasc(klon,klev)
    7979      real zqla(klon,klev)
    80       real zqta(klon,klev)
    8180      real ztv(klon,klev),ztva(klon,klev)
    8281      real zpspsk(klon,klev)

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

    r4654 r4657  
    159159!!$!          ycdragm(1:knon) = ycoefh(1:knon,1)
    160160!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    161        CALL ustarhb(knon,yu,yv,ycdragm, yustar)
     161
     162       ! Normalement, on peut passer dans les codes avec knon=0
     163       ! Mais ca fait planter le replay.
     164       ! En attendant une réécriture, on a joute des if (Fredho)
     165       if ( klon>1 .or. (klon==1 .and. knon==1) ) then
     166          CALL ustarhb(knon,klev,knon,yu,yv,ycdragm, yustar)
     167       endif
    162168     
    163169       IF (prt_level > 9) THEN
    … …  
    167173!   iflag_pbl peut etre utilise comme longuer de melange
    168174       IF (iflag_pbl.GE.31) THEN
    169           if (knon>0) then
    170           CALL vdif_kcay(knon,klev,dtime,RG,RD,ypaprs,yt, &
     175          if ( klon>1 .or. (klon==1 .and. knon==1) ) then
     176          CALL vdif_kcay(knon,klev,knon,dtime,RG,RD,ypaprs,yt, &
    171177               yzlev,yzlay,yu,yv,yteta, &
    172178               ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &

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

    r4653 r4657  
     1
    12!
    23! $Id$
    … …  
    35133514       !  ====
    35143515       IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
     3516       WRITE(lunout,*) 'WARNING : running without dry convection. Somme intermediate variables are not properly defined in physiq_mod.F90'
     3517       ! Reprendre proprement les initialisation ci dessouds si on veut vraiment utiliser l'option (FH)
     3518          fraca(:,:)=0.
     3519          fm_therm(:,:)=0.
     3520          ztv(:,:)=t_seri(:,:)
     3521          zqasc(:,:)=q_seri(:,:)
     3522          ztla(:,:)=0.
     3523          zthl(:,:)=0.
     3524          zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
     3525
    35153526
    35163527

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

    r4655 r4657  
    99     &            debut,lafin,pdtphys, &
    1010     &            paprs,pplay,pphi,pphis,presnivs, &
    11      &            u,v,rot,t,qx, &
     11     &            u,v,rot,temp,qx, &
    1212     &            flxmass_w, &
    1313     &            d_u, d_v, d_t, d_qx, d_ps)
    1414
    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
     15
     16USE dimphy, only : klon,klev
     17USE infotrac_phy, only : nqtot
     18USE geometry_mod, only : latitude
     19USE ioipsl, only : ymds2ju
     20USE phys_state_var_mod, only : phys_state_var_init
     21USE phyetat0_mod, only: phyetat0
     22USE output_physiqex_mod, ONLY: output_physiqex
     23use vdif_ini, only : vdif_ini_
     24USE lmdz_thermcell_ini, ONLY : thermcell_ini
     25USE ioipsl_getin_p_mod, ONLY : getin_p
     26USE wxios, ONLY: missing_val, using_xios
     27USE lscp_mod, ONLY : lscp
     28USE lscp_ini_mod, ONLY : lscp_ini
     29USE add_phys_tend_mod, ONLY : add_phys_tend
     30
    2331
    2432      IMPLICIT none
     33
     34include "YOETHF.h"
     35
     36
     37
     38
    2539!
    2640! Routine argument:
    … …  
    3246      logical,intent(in) :: lafin ! signals last call to physics
    3347      real,intent(in) :: pdtphys ! physics time step (s)
    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
     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
     83real, dimension(klon) :: capcal,z0m,z0h,dtsrf,emis,fluxsrf,cdh,cdv,tsrf_
     84real, dimension(klon,klev) :: zzlay,masse
     85real, dimension(klon,klev+1) :: zzlev,kz_v,kz_h,richardson
     86
     87real, save, allocatable, dimension(:) :: tsrf,f0,zmax0
     88real, 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
     95real :: z1,z2,tau_thermals
     96logical :: lwrite
     97integer :: iflag_replay
     98
     99integer :: 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
     153include "YOMCST.h"
    50154
    51155!    include "clesphys.h"
    … …  
    57161    !$OMP THREADPRIVATE(clesphy0)
    58162
    59 
    60 real :: temp_newton(klon,klev)
    61 integer :: k
     163real,dimension(klon,klev) :: temp_newton
     164integer :: i,k,iq
     165INTEGER, SAVE :: itap=0
     166!$OMP THREADPRIVATE(itap)
     167INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
     168!$OMP THREADPRIVATE(abortphy)
     169
     170integer, 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
    62173logical, save :: first=.true.
    63174!$OMP THREADPRIVATE(first)
    64 
    65 real,save :: rg=9.81
    66 !$OMP THREADPRIVATE(rg)
    67175
    68176! For I/Os
    69177integer :: itau0
    70178real :: zjulian
     179real,dimension(klon,klev) :: du0,dv0,dqbs0
     180real,dimension(klon,klev) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
    71181
    72182
    … …  
    93203  ! Initialize IOIPSL output file
    94204#endif
     205call suphel
     206call vdif_ini_(klon,RCPD,RD,RG,RKAPPA)
     207! Pourquoi ce tau_thermals en argument ??? AFAIRE
     208tau_thermals=0.
     209call getin_p('iflag_thermals',iflag_thermals)
     210
     211call getin_p('iflag_newton',iflag_newton)
     212call getin_p('iflag_reevap',iflag_reevap)
     213call getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
     214call getin_p('iflag_ratqs',iflag_ratqs)
     215call getin_p('iflag_vdif',iflag_vdif)
     216call getin_p('iflag_lscp',iflag_lscp)
     217call getin_p('ratqsbas',ratqsbas)
     218call getin_p('ratqshaut',ratqshaut)
     219call getin_p('ratqsp0',ratqsp0)
     220call getin_p('ratqsdp',ratqsdp)
     221CALL thermcell_ini(iflag_thermals,0,tau_thermals,6, &
     222   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
     223CALL lscp_ini(pdtphys,.false.,iflag_ratqs, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
     224
     225
     226
     227allocate(tsrf(klon),q2(klon,klev+1),f0(klon),zmax0(klon))
     228allocate(u_prev(klon,klev),v_prev(klon,klev),t_prev(klon,klev),q_prev(klon,klev,nqtot))
     229
     230u_prev(:,:)=u(:,:)
     231v_prev(:,:)=v(:,:)
     232t_prev(:,:)=temp(:,:)
     233q_prev(:,:,:)=qx(:,:,:)
     234
     235q2=1.e-10
     236tsrf=temp(:,1)
     237f0=0.
     238zmax0=0.
     239
     240iflag_replay=0
     241call getin_p('iflag_replay',iflag_replay)
     242if ( iflag_replay >= 0 ) CALL iophys_ini(pdtphys)
     243
    95244
    96245endif ! of if (debut)
    … …  
    100249!------------------------------------------------------------
    101250
     251d_u_dyn(:,:)=(u(:,:)-u_prev(:,:))/pdtphys
     252d_v_dyn(:,:)=(v(:,:)-v_prev(:,:))/pdtphys
     253d_t_dyn(:,:)=(temp(:,:)-t_prev(:,:))/pdtphys
     254d_q_dyn(:,:,:)=(qx(:,:,:)-q_prev(:,:,:))/pdtphys
    102255
    103256! set all tendencies to zero
    … …  
    108261d_ps(1:klon)=0.
    109262
     263u_loc(1:klon,1:klev)=u(1:klon,1:klev)
     264v_loc(1:klon,1:klev)=v(1:klon,1:klev)
     265t_loc(1:klon,1:klev)=temp(1:klon,1:klev)
     266d_u_loc(1:klon,1:klev)=0.
     267d_v_loc(1:klon,1:klev)=0.
     268d_t_loc(1:klon,1:klev)=0.
     269do 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
     275enddo
     276
     277du0(1:klon,1:klev)=0.
     278dv0(1:klon,1:klev)=0.
     279dqbs0(1:klon,1:klev)=0.
     280
     281
     282
    110283!------------------------------------------------------------
    111284! Calculs
    112285!------------------------------------------------------------
    113286
    114 ! compute tendencies to return to the dynamics:
    115 ! "friction" on the first layer
    116 d_u(1:klon,1)=-u(1:klon,1)/86400.
    117 d_v(1:klon,1)=-v(1:klon,1)/86400.
    118 ! newtonian relaxation towards temp_newton()
     287!------------------------------------------------------------
     288! Rappel en temperature et frottement dans la premiere chouche
     289!------------------------------------------------------------
     290
     291if ( 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
     301else
     302   temp_newton(:,:)=0.
     303endif
     304
     305
     306!------------------------------------------------------------
     307! Reevaporation de la pluie
     308!------------------------------------------------------------
     309
     310iflag_ice_thermo=1
     311if ( 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
     324else
     325     d_t_eva(:,:)=0.
     326     d_q_eva(:,:)=0.
     327     d_ql_eva(:,:)=0.
     328     d_qi_eva(:,:)=0.
     329endif
     330
     331!-----------------------------------------------------------------------
     332! Variables intermédiaires (altitudes, temperature potentielle ...)
     333!-----------------------------------------------------------------------
     334
     335DO k=1,klev
     336   DO i=1,klon
     337      zzlay(i,k)=pphi(i,k)/rg
     338   ENDDO
     339ENDDO
     340DO i=1,klon
     341   zzlev(i,1)=0.
     342ENDDO
     343DO 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
     349ENDDO
     350
     351!   Transformation de la temperature en temperature potentielle
     352DO k=1,klev
     353   DO i=1,klon
     354!      zpspsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
     355zpspsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
     356masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
     357ENDDO
     358ENDDO
     359DO k=1,klev
     360DO i=1,klon
     361h_loc(i,k)=t_loc(i,k)/zpspsk(i,k)
     362d_h_loc(i,k)=d_t_loc(i,k)/zpspsk(i,k)
     363d_q_loc(i,k,1)=0.
     364ENDDO
     365ENDDO
     366
     367!-----------------------------------------------------------------------
     368! Diffusion verticale
     369!-----------------------------------------------------------------------
     370
     371if ( iflag_vdif == 1 ) then
     372emis(:)=1.
     373!tsrf=300.
     374z0m=0.035
     375z0h=0.035
     376capcal=1e2
     377lwrite=.false.
     378print*,'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)
    119387do k=1,klev
    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
     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
    122395enddo
    123 
     396do i=1,klon
     397 tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys
     398enddo
     399else
     400d_u_vdif(:,:)=0.
     401d_v_vdif(:,:)=0.
     402d_t_vdif(:,:)=0.
     403d_h_vdif(:,:)=0.
     404d_q_vdif(:,:,1)=0.
     405kz_v(:,:)=0.
     406kz_h(:,:)=0.
     407richardson(:,:)=0.
     408endif
     409
     410!-----------------------------------------------------------------------
     411! Thermiques
     412!-----------------------------------------------------------------------
     413
     414do k=1,klev
     415do i=1,klon
     416d_u_the(i,k)=0.
     417d_v_the(i,k)=0.
     418d_t_the(i,k)=0.
     419d_q_the(i,k,1)=0.
     420enddo
     421enddo
     422
     423if ( 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
     443else
     444  frac_the(:,:)=0.
     445  fm_the(:,:)=0.
     446  ztv(:,:)=t_loc(:,:)
     447  zqta(:,:)=q_loc(:,:,1)
     448  ztla(:,:)=0.
     449  zthl(:,:)=0.
     450endif
     451
     452!-----------------------------------------------------------------------
     453! Condensation grande échelle
     454!-----------------------------------------------------------------------
     455
     456iflag_cld_th=5
     457ok_ice_sursat=.false.
     458ptconv(:,:)=.false.
     459distcltop=0.
     460temp_cltop=0.
     461beta(:,:)=1.
     462rneb_seri(:,:)=0.
     463do k=1,klev
     464  ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
     465  *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
     466enddo
     467
     468
     469if ( 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
     495else
     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.
     529endif
     530
     531
     532d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys
     533d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys
     534d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys
     535d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys
    124536
    125537!------------------------------------------------------------
    … …  
    128540
    129541
    130 call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t)
     542tsrf_(:)=tsrf(:)
     543if ( 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)
     545else 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
     550include "physiqex_out.h"
     551
     552endif
    131553
    132554
    … …  
    136558endif
    137559
     560print*,'Fin physiqex'
    138561
    139562end subroutine physiqex

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

    r2346 r4657  
    22! $Header$
    33
    4 SUBROUTINE ustarhb(knon, u, v, cd_m, ustar)
    5   USE dimphy
     4SUBROUTINE ustarhb(klon, klev, knon, u, v, cd_m, ustar)
    65  IMPLICIT NONE
    76  ! ======================================================================
    … …  
    1817  ! model. J. of Climate, vol. 6, 1825-1842.
    1918  ! ======================================================================
    20   include "YOMCST.h"
    2119
    2220  ! Arguments:
    2321
    24   INTEGER knon ! nombre de points a calculer
    25   REAL u(klon, klev) ! vitesse U (m/s)
    26   REAL v(klon, klev) ! vitesse V (m/s)
    27   REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
    28   REAL ustar(klon)
     22  INTEGER, INTENT(IN) :: klon, klev, knon ! nombre de points a calculer
     23  REAL, DIMENSION(klon, klev), INTENT(IN) :: u,v ! vent horizontal (m/s)
     24  REAL, DIMENSION(klon), INTENT(IN) :: cd_m ! coefficient de friction au sol pour vitesse
     25  REAL, DIMENSION(klon), INTENT(OUT) :: ustar
    2926
    30   INTEGER i, k
    31   REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
    32   REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
    33   LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
    34   LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
    35   LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
     27  INTEGER :: i, k
     28  REAL :: zxt, zxq, zxu, zxv, zxmod, taux, tauy
     29  REAL :: zx_alf1, zx_alf2 ! parametres pour extrapolation
    3630
    37   include "YOETHF.h"
    38   include "FCTTRE.h"
    3931  DO i = 1, knon
    4032    zx_alf1 = 1.0
    … …  
    4638    tauy = zxv*zxmod*cd_m(i)
    4739    ustar(i) = sqrt(taux**2+tauy**2)
    48     ! print*,'Ust ',zxu,zxmod,taux,ustar(i)
    4940  END DO
    5041

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

    r4654 r4657  
    22! $Header$
    33
    4 SUBROUTINE vdif_kcay(klon,klev, dt, g, rconst, plev, temp, zlev, zlay, u, v, &
     4SUBROUTINE vdif_kcay(klon,klev,ngrid,dt, g, rconst, plev, temp, zlev, zlay, u, v, &
    55    teta, cd, q2, q2diag, km, kn, ustar, l_mix)
    66
    … …  
    2828
    2929  ! .......................................................................
    30   INTEGER, INTENT(IN) :: klon,klev
     30  INTEGER, INTENT(IN) :: klon,klev,ngrid
    3131  REAL,INTENT(IN) :: dt, g, rconst
    3232  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: plev
    3333  REAL,DIMENSION(klon,klev),INTENT(IN) :: temp
    34   REAL,DIMENSION(klon),INTENT(IN) :: ustar(klon)
     34  REAL,DIMENSION(klon),INTENT(IN) :: ustar
    3535  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: zlev
    3636  REAL,DIMENSION(klon,klev),INTENT(IN) :: zlay
    … …  
    3838  REAL,DIMENSION(klon,klev),INTENT(IN) :: v
    3939  REAL,DIMENSION(klon,klev),INTENT(IN) :: teta
    40   REAL,DIMENSION(klon),INTENT(IN) :: cd(klon)
     40  REAL,DIMENSION(klon),INTENT(IN) :: cd
    4141  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2
    42   REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2diag
     42  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: q2diag
    4343  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: km
    4444  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: kn
    … …  
    6262
    6363  ! .......................................................................
    64   INTEGER :: nlay, nlev, ngrid
     64  INTEGER :: nlay, nlev
    6565  REAL, DIMENSION(klon,klev) :: unsdz
    6666  REAL, DIMENSION(klon, klev+1) :: unsdzdec,q
    … …  
    200200
    201201  LOGICAL :: first
    202   SAVE first
    203   DATA first/.TRUE./
    204   !$OMP THREADPRIVATE(first)
     202!  SAVE first
     203!  DATA first/.TRUE./
     204!  !$OMP THREADPRIVATE(first)
    205205  ! .......................................................................
    206206  ! traitment des valeur de q2 en entree
    … …  
    208208
    209209  ! Initialisation de q2
    210 ngrid=klon
    211210  nlay = klev
    212211  nlev = klev + 1
    … …  
    360359      gninf = .FALSE.
    361360      gnsup = .FALSE.
    362       long(igrid, ilev) = long(igrid, ilev)
    363       long(igrid, ilev) = long(igrid, ilev)
    364361
    365362      IF (gn<gnmin) THEN
Note: See TracChangeset for help on using the changeset viewer.