Ignore:
Timestamp:
Dec 14, 2015, 11:43:09 AM (9 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2298:2396 into testing branch

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/physiq.F90

    r2298 r2408  
    33
    44SUBROUTINE physiq (nlon,nlev, &
    5      debut,lafin,jD_cur, jH_cur,pdtphys, &
     5     debut,lafin,jD_cur_,jH_cur_,pdtphys_, &
    66     paprs,pplay,pphi,pphis,presnivs, &
    7      u,v,t,qx, &
     7     u,v,rot,t,qx, &
    88     flxmass_w, &
    99     d_u, d_v, d_t, d_qx, d_ps &
    … …  
    1212  USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
    1313       histwrite, ju2ymds, ymds2ju, getin
    14   USE comgeomphy
     14  USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
    1515  USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, year_cur, &
    16        mth_cur, phys_cal_update
     16       mth_cur,jD_cur, jH_cur, jD_ref, phys_cal_update
    1717  USE write_field_phy
    1818  USE dimphy
    19   USE infotrac
    20   USE mod_grid_phy_lmdz
     19  USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
     20  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo
    2121  USE mod_phys_lmdz_para
    2222  USE iophy
    23   USE misc_mod, mydebug=>debug
     23  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
     24  USE phystokenc_mod, ONLY: offline, phystokenc
     25  USE time_phylmdz_mod, only: raz_date, day_step_phy
    2426  USE vampir
    2527  USE pbl_surface_mod, ONLY : pbl_surface
    … …  
    4648  use radlwsw_m, only: radlwsw
    4749  use phyaqua_mod, only: zenang_an
    48   USE control_mod
     50  USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
     51                              start_time, pdtphys, day_ini
     52  USE tracinca_mod, ONLY: config_inca
    4953#ifdef CPP_XIOS
    5054  USE wxios, ONLY: missing_val, missing_val_omp
    … …  
    6064  USE YOERAD   , ONLY : NRADLP
    6165#endif
     66  USE ioipsl_getin_p_mod, ONLY : getin_p
     67
    6268
    6369  !IM stations CFMIP
    6470  USE CFMIP_point_locations
    6571  use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
     72  use ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
    6673
    6774  IMPLICIT none
    … …  
    109116  !! d_ps----output-R-tendance physique de la pression au sol
    110117  !!======================================================================
    111   include "dimensions.h"
    112118  integer jjmp1
    113   parameter (jjmp1=jjm+1-1/jjm)
    114   integer iip1
    115   parameter (iip1=iim+1)
     119!  parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
     120!  integer iip1
     121!  parameter (iip1=iim+1)
    116122
    117123  include "regdim.h"
    118124  include "dimsoil.h"
    119125  include "clesphys.h"
    120   include "temps.h"
    121   include "iniprint.h"
    122126  include "thermcell.h"
    123127  !======================================================================
    … …  
    139143  !======================================================================
    140144  ! Clef controlant l'activation du cycle diurne:
    141   !cc      LOGICAL cycle_diurne
    142   !cc      PARAMETER (cycle_diurne=.FALSE.)
     145  ! en attente du codage des cles par Fred
     146        INTEGER iflag_cycle_diurne
     147        PARAMETER (iflag_cycle_diurne=1)
    143148  !======================================================================
    144149  ! Modele thermique du sol, a activer pour le cycle diurne:
    … …  
    215220  INTEGER nlon
    216221  INTEGER nlev
    217   REAL, intent(in):: jD_cur, jH_cur
    218 
    219   REAL pdtphys
     222  REAL, intent(in):: jD_cur_, jH_cur_
     223! JD_cur and JH_cur to be used in physics are in phys_cal_mod
     224  REAL,INTENT(IN) :: pdtphys_
     225! NB: pdtphys to be used in physics is in time_phylmdz_mod
    220226  LOGICAL debut, lafin
    221227  REAL paprs(klon,klev+1)
    … …  
    229235  REAL u(klon,klev)
    230236  REAL v(klon,klev)
     237
     238  REAL, intent(in):: rot(klon, klev)
     239  ! relative vorticity, in s-1, needed for frontal waves
     240
    231241  REAL t(klon,klev),thetal(klon,klev)
    232242  ! thetal: ligne suivante a decommenter si vous avez les fichiers     MPL 20130625
    … …  
    404414  REAL dt_a(klon,klev)
    405415  REAL dq_a(klon,klev)
     416  REAL d_t_adjwk(klon,klev)                !jyg
     417  REAL d_q_adjwk(klon,klev)                !jyg
     418  LOGICAL,SAVE :: ok_adjwk=.FALSE.
     419  !$OMP THREADPRIVATE(ok_adjwk)
    406420  REAL, dimension(klon) :: www
    407421  REAL, SAVE :: alp_offset
    … …  
    440454  integer :: tau_trig(klon)
    441455
     456  REAL,SAVE :: random_notrig_max=1.
     457  !$OMP THREADPRIVATE(random_notrig_max)
     458
    442459  !--------Statistical Boundary Layer Closure: ALP_BL--------
    443460  !---Profils de TKE dans et hors du thermique
    … …  
    494511  SAVE lmt_pas                ! frequence de mise a jour
    495512  !$OMP THREADPRIVATE(lmt_pas)
    496   real zmasse(klon, llm),exner(klon, llm)
     513  real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
    497514  !     (column-density of mass of air in a cell, in kg m-2)
    498515  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
    … …  
    528545  EXTERNAL suphel    ! initialiser certaines constantes
    529546  EXTERNAL transp    ! transport total de l'eau et de l'energie
    530   EXTERNAL ecribina  ! ecrire le fichier binaire global
    531   EXTERNAL ecribins  ! ecrire le fichier binaire global
    532   EXTERNAL ecrirega  ! ecrire le fichier binaire regional
    533   EXTERNAL ecriregs  ! ecrire le fichier binaire regional
    534547  !IM
    535548  EXTERNAL haut2bas  !variables de haut en bas
    … …  
    574587  !  REAL zxsnow(klon)
    575588  REAL zxsnow_dummy(klon)
     589  REAL zsav_tsol(klon)
    576590  !
    577591  REAL dist, rmu0(klon), fract(klon)
    578   REAL zdtime, zlongi
     592  REAL zrmu0(klon), zfract(klon)
     593  REAL zdtime, zdtime1, zdtime2, zlongi
    579594  !
    580595  REAL qcheck
    … …  
    682697  REAL tabcntr0( length       )
    683698  !
    684   INTEGER ndex2d(iim*jjmp1)
     699  INTEGER ndex2d(nbp_lon*nbp_lat)
    685700  !IM
    686701  !
    … …  
    691706  REAL zustrli(klon), zvstrli(klon)
    692707  REAL zustrph(klon), zvstrph(klon)
    693   REAL zustrhi(klon), zvstrhi(klon)
    694708  REAL aam, torsfc
    695709  !IM 141004 END
    696710  !IM 190504 BEG
    697   INTEGER ij, imp1jmp1
    698   PARAMETER(imp1jmp1=(iim+1)*jjmp1)
     711  INTEGER ij
     712!  INTEGER imp1jmp1
     713!  PARAMETER(imp1jmp1=(iim+1)*jjmp1)
    699714  !ym A voir plus tard
    700   REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
    701   REAL padyn(iim+1,jjmp1,klev+1)
    702   REAL dudyn(iim+1,jjmp1,klev)
    703   REAL rlatdyn(iim+1,jjmp1)
     715  REAL zx_tmp((nbp_lon+1)*nbp_lat)
     716  REAL airedyn(nbp_lon+1,nbp_lat)
     717  REAL padyn(nbp_lon+1,nbp_lat,klev+1)
     718  REAL dudyn(nbp_lon+1,nbp_lat,klev)
     719  REAL rlatdyn(nbp_lon+1,nbp_lat)
    704720  !IM 190504 END
    705721  LOGICAL ok_msk
    … …  
    713729  REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
    714730  REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
    715   REAL zx_tmp_2d(iim,jjmp1)
    716   REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
     731  REAL zx_tmp_2d(nbp_lon,nbp_lat)
     732  REAL zx_lon(nbp_lon,nbp_lat)
     733  REAL zx_lat(nbp_lon,nbp_lat)
    717734  !
    718735  INTEGER nid_day_seri, nid_ctesGCM
    … …  
    889906!albedo SB <<<
    890907
     908  ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
     909  jjmp1=nbp_lat
     910
    891911  !======================================================================
    892912  ! Gestion calendrier : mise a jour du module phys_cal_mod
    893913  !
     914  JD_cur=JD_cur_
     915  JH_cur=JH_cur_
     916  pdtphys=pdtphys_
    894917  CALL phys_cal_update(jD_cur,jH_cur)
    895918
    … …  
    902925     igout=klon/2+1/klon
    903926     write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
    904      write(lunout,*) 'igout, rlat, rlon ',igout, rlatd(igout)*180./3.141593, rlond(igout)*180./3.141593
     927     write(lunout,*) 'igout, rlat, rlon ',igout, latitude_deg(igout), longitude_deg(igout)
    905928     write(lunout,*) &
    906929          'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
    … …  
    975998  pde_u(:,:) = 0.
    976999  aam=0.
     1000  d_t_adjwk(:,:)=0
     1001  d_q_adjwk(:,:)=0
    9771002
    9781003  alp_bl_conv(:)=0.
    9791004
    9801005  torsfc=0.
    981   forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
     1006  forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
    9821007
    9831008
    … …  
    9931018  IF (debut) THEN
    9941019     CALL suphel ! initialiser constantes et parametres phys.
     1020     CALL getin_p('random_notrig_max',random_notrig_max)
     1021     CALL getin_p('ok_adjwk',ok_adjwk)
    9951022  ENDIF
    9961023
    … …  
    10331060        piz_aero(:,:,:,:) = 0.
    10341061        cg_aero(:,:,:,:) = 0.
     1062
     1063        config_inca='none' ! default
     1064        CALL getin_p('config_inca',config_inca)
     1065
    10351066     END IF
    10361067
    … …  
    10451076     print*,'iflag_coupl,iflag_clos,iflag_wake', &
    10461077          iflag_coupl,iflag_clos,iflag_wake
    1047      print*,'CYCLE_DIURNE', cycle_diurne
     1078     print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
    10481079     !
    10491080     IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
    10501081        abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
    1051         CALL abort_gcm (modname,abort_message,1)
     1082        CALL abort_physic (modname,abort_message,1)
    10521083     ENDIF
    10531084     !
    … …  
    10731104     ! pour obtenir le meme resultat.
    10741105     dtime=pdtphys
    1075      radpas = NINT( 86400./dtime/nbapp_rad)
     1106     IF (MOD(INT(86400./dtime),nbapp_rad).EQ.0) THEN
     1107       radpas = NINT( 86400./dtime/nbapp_rad)
     1108     ELSE
     1109       WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un multiple de nbapp_rad'
     1110       WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test mais 1+1<>2'
     1111       abort_message='nbre de pas de temps physique n est pas multiple de nbapp_rad'
     1112       call abort_physic(modname,abort_message,1)
     1113     ENDIF
    10761114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    10771115
    … …  
    10971135     ENDIF
    10981136
    1099      !IM cf. AM 081204 BEG
    1100      PRINT*,'cycle_diurne3 =',cycle_diurne
    1101      !IM cf. AM 081204 END
    1102      !
    11031137     CALL printflag( tabcntr0,radpas,ok_journe, &
    11041138          ok_instan, ok_region )
    … …  
    11081142             pdtphys
    11091143        abort_message='Pas physique n est pas correct '
    1110         !           call abort_gcm(modname,abort_message,1)
     1144        !           call abort_physic(modname,abort_message,1)
    11111145        dtime=pdtphys
    11121146     ENDIF
    … …  
    11151149             klon
    11161150        abort_message='nlon et klon ne sont pas coherents'
    1117         call abort_gcm(modname,abort_message,1)
     1151        call abort_physic(modname,abort_message,1)
    11181152     ENDIF
    11191153     IF (nlev .NE. klev) THEN
    … …  
    11211155             klev
    11221156        abort_message='nlev et klev ne sont pas coherents'
    1123         call abort_gcm(modname,abort_message,1)
     1157        call abort_physic(modname,abort_message,1)
    11241158     ENDIF
    11251159     !
    1126      IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
     1160     IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
    11271161        WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
    11281162        WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
    11291163        abort_message='Nbre d appels au rayonnement insuffisant'
    1130         call abort_gcm(modname,abort_message,1)
     1164        call abort_physic(modname,abort_message,1)
    11311165     ENDIF
    11321166     WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
    … …  
    11831217           IF(nCFMIP.GT.npCFMIP) THEN
    11841218              print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
    1185               call abort_gcm("physiq", "", 1)
     1219              call abort_physic("physiq", "", 1)
    11861220           else
    11871221              print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
    … …  
    13311365             rg, &
    13321366             ra, &
    1333              airephy, &
     1367             cell_area, &
    13341368             rlat, &
    13351369             rlon, &
    … …  
    15381572  IF (ip_ebil_phy.ge.1) THEN
    15391573     ztit='after dynamic'
    1540      CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
     1574     CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
    15411575          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    15421576          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    … …  
    15451579     !     est egale a la variation de la physique au pas de temps precedent.
    15461580     !     Donc la somme de ces 2 variations devrait etre nulle.
    1547      call diagphy(airephy,ztit,ip_ebil_phy &
     1581     call diagphy(cell_area,ztit,ip_ebil_phy &
    15481582          , zero_v, zero_v, zero_v, zero_v, zero_v &
    15491583          , zero_v, zero_v, zero_v, ztsol &
    … …  
    16881722        IF (read_climoz/=-1) THEN
    16891723           abort_message ='read_climoz=-1 is recommended when solarlong0=1000.'
    1690            CALL abort_gcm (modname,abort_message,1)
     1724           CALL abort_physic (modname,abort_message,1)
    16911725        ENDIF
    16921726     ELSE
    … …  
    17491783  IF (ip_ebil_phy.ge.2) THEN
    17501784     ztit='after reevap'
    1751      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
     1785     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,1,dtime &
    17521786          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    17531787          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    1754      call diagphy(airephy,ztit,ip_ebil_phy &
     1788     call diagphy(cell_area,ztit,ip_ebil_phy &
    17551789          , zero_v, zero_v, zero_v, zero_v, zero_v &
    17561790          , zero_v, zero_v, zero_v, ztsol &
    … …  
    17961830  ! non nul aux poles.
    17971831  IF (abs(solarlong0-1000.)<1.e-4) then
    1798      call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
     1832     call zenang_an(iflag_cycle_diurne.GE.1,jH_cur,rlat,rlon,rmu0,fract)
     1833     JrNt = 1.0
    17991834  ELSE
    1800      !  Avec ou sans cycle diurne
    1801      IF (cycle_diurne) THEN
     1835  ! recode par Olivier Boucher en sept 2015
     1836     SELECT CASE (iflag_cycle_diurne)
     1837     CASE(0) 
     1838     !  Sans cycle diurne
     1839        CALL angle(zlongi, rlat, fract, rmu0)
     1840        swradcorr = 1.0
     1841        JrNt = 1.0
     1842        zrmu0 = rmu0
     1843     CASE(1) 
     1844     !  Avec cycle diurne sans application des poids
     1845     !  bit comparable a l ancienne formulation cycle_diurne=true
     1846     !  on integre entre gmtime et gmtime+radpas
    18021847        zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
    1803         CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
    1804      ELSE
    1805         CALL angle(zlongi, rlat, fract, rmu0)
    1806      ENDIF
     1848        CALL zenang(zlongi,jH_cur,0.0,zdtime,rlat,rlon,rmu0,fract)
     1849        zrmu0 = rmu0
     1850        swradcorr = 1.0
     1851     ! Calcul du flag jour-nuit
     1852        JrNt = 0.0
     1853        WHERE (fract.GT.0.0) JrNt = 1.0
     1854     CASE(2) 
     1855     !  Avec cycle diurne sans application des poids
     1856     !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
     1857     !  Comme cette routine est appele a tous les pas de temps de la physique
     1858     !  meme si le rayonnement n'est pas appele je remonte en arriere les
     1859     !  radpas-1 pas de temps suivant. Petite ruse avec MOD pour prendre en
     1860     !  compte le premier pas de temps de la physique pendant lequel itaprad=0
     1861        zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)     
     1862        zdtime2=dtime*REAL(radpas-MOD(itaprad,radpas)-1)
     1863        CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,rmu0,fract)
     1864     !
     1865     ! Calcul des poids
     1866     !
     1867        zdtime1=-dtime !--on corrige le rayonnement pour representer le
     1868        zdtime2=0.0    !--pas de temps de la physique qui se termine
     1869        CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,zrmu0,zfract)
     1870        swradcorr = 0.0
     1871        WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) swradcorr=zfract/fract*zrmu0/rmu0
     1872     ! Calcul du flag jour-nuit
     1873        JrNt = 0.0
     1874        WHERE (zfract.GT.0.0) JrNt = 1.0
     1875     END SELECT
    18071876  ENDIF
    18081877
    1809   ! AI Janv 2014
    1810   do i = 1, klon
    1811      if (fract(i).le.0.) then
    1812         JrNt(i)=0.
    1813      else
    1814         JrNt(i)=1.
    1815      endif
    1816   enddo
    1817 
    18181878  if (mydebug) then
    1819      call writefield_phy('u_seri',u_seri,llm)
    1820      call writefield_phy('v_seri',v_seri,llm)
    1821      call writefield_phy('t_seri',t_seri,llm)
    1822      call writefield_phy('q_seri',q_seri,llm)
     1879     call writefield_phy('u_seri',u_seri,nbp_lev)
     1880     call writefield_phy('v_seri',v_seri,nbp_lev)
     1881     call writefield_phy('t_seri',t_seri,nbp_lev)
     1882     call writefield_phy('q_seri',q_seri,nbp_lev)
    18231883  endif
    18241884
    … …  
    18871947          dtime,     date0,     itap,    days_elapsed+1, &
    18881948          debut,     lafin, &
    1889           rlon,      rlat,      rugoro,  rmu0,      &
     1949          rlon,      rlat,      rugoro,  zrmu0,      &
    18901950          zsig,      sollwdown, pphi,    cldt,      &
    18911951          rain_fall, snow_fall, solsw,   sollw,     &
    … …  
    19632023
    19642024     if (mydebug) then
    1965         call writefield_phy('u_seri',u_seri,llm)
    1966         call writefield_phy('v_seri',v_seri,llm)
    1967         call writefield_phy('t_seri',t_seri,llm)
    1968         call writefield_phy('q_seri',q_seri,llm)
     2025        call writefield_phy('u_seri',u_seri,nbp_lev)
     2026        call writefield_phy('v_seri',v_seri,nbp_lev)
     2027        call writefield_phy('t_seri',t_seri,nbp_lev)
     2028        call writefield_phy('q_seri',q_seri,nbp_lev)
    19692029     endif
    19702030
    … …  
    20072067     IF (ip_ebil_phy.ge.2) THEN
    20082068        ztit='after surface_main'
    2009         CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     2069        CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    20102070             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    20112071             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2012         call diagphy(airephy,ztit,ip_ebil_phy &
     2072        call diagphy(cell_area,ztit,ip_ebil_phy &
    20132073             , zero_v, zero_v, zero_v, zero_v, sens &
    20142074             , evap  , zero_v, zero_v, ztsol &
    … …  
    20572117  ENDDO
    20582118  IF (check) THEN
    2059      za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
     2119     za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
    20602120     WRITE(lunout,*) "avantcon=", za
    20612121  ENDIF
    … …  
    20772137  DO k = 1, klev
    20782138     DO i = 1, klon
    2079         omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
     2139        omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
    20802140     END DO
    20812141  END DO
    … …  
    20852145  IF (iflag_con.EQ.1) THEN
    20862146     abort_message ='reactiver le call conlmd dans physiq.F'
    2087      CALL abort_gcm (modname,abort_message,1)
     2147     CALL abort_physic (modname,abort_message,1)
    20882148     !     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
    20892149     !    .             d_t_con, d_q_con,
    … …  
    21332193        enddo
    21342194     enddo
     2195!
     2196!jyg<
     2197     ! Perform dry adiabatic adjustment on wake profile
     2198     ! The corresponding tendencies are added to the convective tendencies
     2199     ! after the call to the convective scheme.
     2200     IF (iflag_wake>=1) then
     2201      IF (ok_adjwk) THEN
     2202        limbas(:) = 1
     2203        CALL ajsec(paprs, pplay, t_wake, q_wake, limbas, &
     2204                  d_t_adjwk, d_q_adjwk)
     2205      ENDIF
     2206!
     2207      DO k=1,klev
     2208        DO i=1,klon
     2209          IF (wake_s(i) .GT. 1.e-3) THEN
     2210            t_wake(i,k) = t_wake(i,k) + d_t_adjwk(i,k)
     2211            q_wake(i,k) = q_wake(i,k) + d_q_adjwk(i,k)
     2212            wake_deltat(i,k) = wake_deltat(i,k) + d_t_adjwk(i,k)
     2213            wake_deltaq(i,k) = wake_deltaq(i,k) + d_q_adjwk(i,k)
     2214          ENDIF
     2215        ENDDO
     2216      ENDDO
     2217     ENDIF ! (iflag_wake>=1)
     2218!>jyg
     2219!
    21352220
    21362221     ! Calcul de l'energie disponible ALE (J/kg) et de la puissance
    … …  
    22232308           else
    22242309       abort_message ='Ne pas passer la car www non calcule'
    2225        CALL abort_gcm (modname,abort_message,1)
     2310       CALL abort_physic (modname,abort_message,1)
    22262311
    22272312!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    23222407           if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
    23232408        enddo
    2324 
     2409!
     2410!jyg<
     2411!    Add the tendency due to the dry adjustment of the wake profile
     2412      IF (iflag_wake>=1) THEN
     2413        DO k=1,klev
     2414          DO i=1,klon
     2415            ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
     2416            fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
     2417            d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
     2418            d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
     2419          ENDDO
     2420        ENDDO
     2421      ENDIF
     2422!>jyg
     2423!
    23252424     ELSE ! ok_cvl
    23262425
    … …  
    24032502  ELSE
    24042503     WRITE(lunout,*) "iflag_con non-prevu", iflag_con
    2405      call abort_gcm("physiq", "", 1)
     2504     call abort_physic("physiq", "", 1)
    24062505  ENDIF
    24072506
    … …  
    24152514
    24162515  if (mydebug) then
    2417      call writefield_phy('u_seri',u_seri,llm)
    2418      call writefield_phy('v_seri',v_seri,llm)
    2419      call writefield_phy('t_seri',t_seri,llm)
    2420      call writefield_phy('q_seri',q_seri,llm)
     2516     call writefield_phy('u_seri',u_seri,nbp_lev)
     2517     call writefield_phy('v_seri',v_seri,nbp_lev)
     2518     call writefield_phy('t_seri',t_seri,nbp_lev)
     2519     call writefield_phy('q_seri',q_seri,nbp_lev)
    24212520  endif
    24222521
    … …  
    24242523  IF (ip_ebil_phy.ge.2) THEN
    24252524     ztit='after convect'
    2426      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     2525     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    24272526          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    24282527          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2429      call diagphy(airephy,ztit,ip_ebil_phy &
     2528     call diagphy(cell_area,ztit,ip_ebil_phy &
    24302529          , zero_v, zero_v, zero_v, zero_v, zero_v &
    24312530          , zero_v, rain_con, snow_con, ztsol &
    … …  
    24352534  !
    24362535  IF (check) THEN
    2437      za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
     2536     za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
    24382537     WRITE(lunout,*)"aprescon=", za
    24392538     zx_t = 0.0
    24402539     za = 0.0
    24412540     DO i = 1, klon
    2442         za = za + airephy(i)/REAL(klon)
     2541        za = za + cell_area(i)/REAL(klon)
    24432542        zx_t = zx_t + (rain_con(i)+ &
    2444              snow_con(i))*airephy(i)/REAL(klon)
     2543             snow_con(i))*cell_area(i)/REAL(klon)
    24452544     ENDDO
    24462545     zx_t = zx_t/za*dtime
    … …  
    25892688  IF (ip_ebil_phy.ge.2) THEN
    25902689     ztit='after wake'
    2591      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     2690     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    25922691          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    25932692          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2594      call diagphy(airephy,ztit,ip_ebil_phy &
     2693     call diagphy(cell_area,ztit,ip_ebil_phy &
    25952694          , zero_v, zero_v, zero_v, zero_v, zero_v &
    25962695          , zero_v, zero_v, zero_v, ztsol &
    … …  
    26822781             ,ztv,zpspsk,ztla,zthl &
    26832782             !cc nrlmd le 10/04/2012
    2684              ,pbl_tke_input,pctsrf,omega,airephy &
     2783             ,pbl_tke_input,pctsrf,omega,cell_area &
    26852784             ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
    26862785             ,n2,s2,ale_bl_stat &
    … …  
    27342833              proba_notrig(i)=1.
    27352834              random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
     2835              if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
    27362836              if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
    27372837                 tau_trig(i)=tau_trig_shallow
    … …  
    29053005  IF (ip_ebil_phy.ge.2) THEN
    29063006     ztit='after dry_adjust'
    2907      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     3007     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    29083008          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    29093009          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2910      call diagphy(airephy,ztit,ip_ebil_phy &
     3010     call diagphy(cell_area,ztit,ip_ebil_phy &
    29113011          , zero_v, zero_v, zero_v, zero_v, zero_v &
    29123012          , zero_v, zero_v, zero_v, ztsol &
    … …  
    29583058  ENDDO
    29593059  IF (check) THEN
    2960      za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
     3060     za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
    29613061     WRITE(lunout,*)"apresilp=", za
    29623062     zx_t = 0.0
    29633063     za = 0.0
    29643064     DO i = 1, klon
    2965         za = za + airephy(i)/REAL(klon)
     3065        za = za + cell_area(i)/REAL(klon)
    29663066        zx_t = zx_t + (rain_lsc(i) &
    2967              + snow_lsc(i))*airephy(i)/REAL(klon)
     3067             + snow_lsc(i))*cell_area(i)/REAL(klon)
    29683068     ENDDO
    29693069     zx_t = zx_t/za*dtime
    … …  
    29733073  IF (ip_ebil_phy.ge.2) THEN
    29743074     ztit='after fisrt'
    2975      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     3075     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    29763076          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    29773077          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2978      call diagphy(airephy,ztit,ip_ebil_phy &
     3078     call diagphy(cell_area,ztit,ip_ebil_phy &
    29793079          , zero_v, zero_v, zero_v, zero_v, zero_v &
    29803080          , zero_v, rain_lsc, snow_lsc, ztsol &
    … …  
    29843084
    29853085  if (mydebug) then
    2986      call writefield_phy('u_seri',u_seri,llm)
    2987      call writefield_phy('v_seri',v_seri,llm)
    2988      call writefield_phy('t_seri',t_seri,llm)
    2989      call writefield_phy('q_seri',q_seri,llm)
     3086     call writefield_phy('u_seri',u_seri,nbp_lev)
     3087     call writefield_phy('v_seri',v_seri,nbp_lev)
     3088     call writefield_phy('t_seri',t_seri,nbp_lev)
     3089     call writefield_phy('q_seri',q_seri,nbp_lev)
    29903090  endif
    29913091
    … …  
    30563156     !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
    30573157     IF (flag_aerosol .gt. 0) THEN
    3058         IF (iflag_rrtm .EQ. 0) THEN !--old radiation
     3158         IF (iflag_rrtm .EQ. 0) THEN !--old radiation
    30593159           IF (.NOT. aerosol_couple) THEN
    30603160              !
    … …  
    30693169           IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
    30703170            abort_message='config_inca=aero et rrtm=1 impossible'
    3071             call abort_gcm(modname,abort_message,1)
     3171            call abort_physic(modname,abort_message,1)
    30723172           ELSE
    30733173!
    30743174#ifdef CPP_RRTM
     3175           IF (NSW.EQ.6) THEN
     3176!--new aerosol properties
     3177!
    30753178             CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
    30763179             new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
    … …  
    30803183             tausum_aero, tau3d_aero)
    30813184
    3082              CALL aeropt_lw_rrtm
     3185           ELSE IF (NSW.EQ.2) THEN
     3186!--for now we use the old aerosol properties
     3187!
     3188              CALL readaerosol_optic( &
     3189                   debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
     3190                   pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
     3191                   mass_solu_aero, mass_solu_aero_pi,  &
     3192                   tau_aero, piz_aero, cg_aero,  &
     3193                   tausum_aero, tau3d_aero)
     3194!
     3195                   !--natural aerosols
     3196                   tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
     3197                   piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
     3198                   cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
     3199                   !--all aerosols
     3200                   tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
     3201                   piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
     3202                   cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
     3203           ELSE
     3204              abort_message='Only NSW=2 or 6 are possible with aerosols and iflag_rrtm=1'
     3205              call abort_physic(modname,abort_message,1)
     3206           ENDIF
     3207
     3208           CALL aeropt_lw_rrtm
     3209!
    30833210#else
    3084 
    3085               abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
    3086               call abort_gcm(modname,abort_message,1)
     3211           abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
     3212           call abort_physic(modname,abort_message,1)
    30873213#endif
    30883214              !
    … …  
    31173243
    31183244           abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
    3119            call abort_gcm(modname,abort_message,1)
     3245           call abort_physic(modname,abort_message,1)
    31203246#endif
    31213247        ENDIF
    … …  
    32173343  IF (ip_ebil_phy.ge.2) THEN
    32183344     ztit="after diagcld"
    3219      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     3345     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    32203346          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    32213347          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    3222      call diagphy(airephy,ztit,ip_ebil_phy &
     3348     call diagphy(cell_area,ztit,ip_ebil_phy &
    32233349          , zero_v, zero_v, zero_v, zero_v, zero_v &
    32243350          , zero_v, zero_v, zero_v, ztsol &
    … …  
    32863412        CALL AEROSOL_METEO_CALC( &
    32873413             calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
    3288              prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
     3414             prfl,psfl,pctsrf,cell_area,rlat,rlon,u10m,v10m)
    32893415     END IF
    32903416
    … …  
    32973423          rlat, &
    32983424          rlon, &
    3299           airephy, &
     3425          cell_area, &
    33003426          paprs, &
    33013427          pplay, &
    … …  
    33163442          ibas_con, &
    33173443          cldfra, &
    3318           iim, &
    3319           jjm, &
     3444          nbp_lon, &
     3445          nbp_lat-1, &
    33203446          tr_seri, &
    33213447          ftsol, &
    … …  
    33463472        IF (ok_cdnc.AND.NRADLP.NE.3) THEN
    33473473           abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
    3348            call abort_gcm(modname,abort_message,1)
     3474           call abort_physic(modname,abort_message,1)
    33493475        endif
    33503476#else
    33513477
    33523478        abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
    3353         call abort_gcm(modname,abort_message,1)
     3479        call abort_physic(modname,abort_message,1)
    33543480#endif
    33553481     ENDIF
    … …  
    33763502  cldtaupirad = cldtaupi
    33773503  cldemirad   = cldemi
     3504  cldfrarad   = cldfra
    33783505
    33793506  !
    … …  
    34443571
    34453572     if (mydebug) then
    3446         call writefield_phy('u_seri',u_seri,llm)
    3447         call writefield_phy('v_seri',v_seri,llm)
    3448         call writefield_phy('t_seri',t_seri,llm)
    3449         call writefield_phy('q_seri',q_seri,llm)
     3573        call writefield_phy('u_seri',u_seri,nbp_lev)
     3574        call writefield_phy('v_seri',v_seri,nbp_lev)
     3575        call writefield_phy('t_seri',t_seri,nbp_lev)
     3576        call writefield_phy('q_seri',q_seri,nbp_lev)
    34503577     endif
     3578
     3579!
     3580!sonia :   If Iflag_radia >=2, pertubation of some variables input to radiation
     3581!(DICE)
     3582!
     3583      IF (iflag_radia .ge. 2) THEN
     3584        zsav_tsol (:) = zxtsol(:)
     3585        call perturb_radlwsw(zxtsol,iflag_radia)
     3586      ENDIF
    34513587
    34523588     IF (aerosol_couple) THEN
    … …  
    34573593             wo(:, :, 1), &
    34583594             cldfrarad, cldemirad, cldtaurad, &
    3459              heat,heat0,cool,cool0,radsol,albpla, &
     3595             heat,heat0,cool,cool0,albpla, &
    34603596             topsw,toplw,solsw,sollw, &
    34613597             sollwdown, &
    … …  
    35033639             zqsat, flwc, fiwc, &
    35043640             ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
    3505              heat,heat0,cool,cool0,radsol,albpla, &
     3641             heat,heat0,cool,cool0,albpla, &
    35063642             topsw,toplw,solsw,sollw, &
    35073643             sollwdown, &
    … …  
    35593695                   zqsat, flwc, fiwc, &
    35603696                   ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
    3561                    heatp,heat0p,coolp,cool0p,radsolp,albplap, &
     3697                   heatp,heat0p,coolp,cool0p,albplap, &
    35623698                   topswp,toplwp,solswp,sollwp, &
    35633699                   sollwdownp, &
    … …  
    35833719     ENDIF ! aerosol_couple
    35843720     itaprad = 0
     3721!
     3722!  If Iflag_radia >=2, reset pertubed variables
     3723!
     3724      IF (iflag_radia .ge. 2) THEN
     3725        zxtsol(:) = zsav_tsol (:)
     3726      ENDIF
    35853727  ENDIF ! MOD(itaprad,radpas)
    35863728  itaprad = itaprad + 1
    … …  
    36003742     swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
    36013743     swup0=0.
    3602      swdn=0.
    3603      swdn0=0.
    36043744     lwup=0.
    36053745     lwup0=0.
    … …  
    36093749
    36103750  !
     3751  ! Calculer radsol a l'exterieur de radlwsw
     3752  ! pour prendre en compte le cycle diurne
     3753  ! recode par Olivier Boucher en sept 2015
     3754  !
     3755  radsol=solsw*swradcorr+sollw
     3756  if (ok_4xCO2atm) then
     3757    radsolp=solswp*swradcorr+sollwp
     3758  endif
     3759
     3760  !
    36113761  ! Ajouter la tendance des rayonnements (tous les pas)
    3612   !
    3613   d_t_swr(:,:)=heat(:,:)*dtime/RDAY
    3614   d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
    3615   d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
    3616   d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
     3762  ! avec une correction pour le cycle diurne dans le SW
     3763  !
     3764 
     3765  DO k=1, klev
     3766    d_t_swr(:,k)=swradcorr(:)*heat(:,k)*dtime/RDAY
     3767    d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*dtime/RDAY
     3768    d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
     3769    d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
     3770  ENDDO
     3771
    36173772  CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy)
    36183773  CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy)
    … …  
    36203775  !
    36213776  if (mydebug) then
    3622      call writefield_phy('u_seri',u_seri,llm)
    3623      call writefield_phy('v_seri',v_seri,llm)
    3624      call writefield_phy('t_seri',t_seri,llm)
    3625      call writefield_phy('q_seri',q_seri,llm)
     3777     call writefield_phy('u_seri',u_seri,nbp_lev)
     3778     call writefield_phy('v_seri',v_seri,nbp_lev)
     3779     call writefield_phy('t_seri',t_seri,nbp_lev)
     3780     call writefield_phy('q_seri',q_seri,nbp_lev)
    36263781  endif
    36273782
    … …  
    36293784  IF (ip_ebil_phy.ge.2) THEN
    36303785     ztit='after rad'
    3631      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     3786     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    36323787          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    36333788          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    3634      call diagphy(airephy,ztit,ip_ebil_phy &
     3789     call diagphy(cell_area,ztit,ip_ebil_phy &
    36353790          , topsw, toplw, solsw, sollw, zero_v &
    36363791          , zero_v, zero_v, zero_v, ztsol &
    … …  
    37053860  !
    37063861  if (mydebug) then
    3707      call writefield_phy('u_seri',u_seri,llm)
    3708      call writefield_phy('v_seri',v_seri,llm)
    3709      call writefield_phy('t_seri',t_seri,llm)
    3710      call writefield_phy('q_seri',q_seri,llm)
     3862     call writefield_phy('u_seri',u_seri,nbp_lev)
     3863     call writefield_phy('v_seri',v_seri,nbp_lev)
     3864     call writefield_phy('t_seri',t_seri,nbp_lev)
     3865     call writefield_phy('q_seri',q_seri,nbp_lev)
    37113866  endif
    37123867
    … …  
    37423897             d_t_lif, d_u_lif, d_v_lif)
    37433898     ENDIF
    3744      !   
    3745      !-----------------------------------------------------------------------------------------
     3899
    37463900     ! ajout des tendances de la portance de l'orographie
    3747      CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif',abortphy)
    3748      !-----------------------------------------------------------------------------------------
    3749      !
     3901     CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
     3902          'lif', abortphy)
    37503903  ENDIF ! fin de test sur ok_orolf
    3751   !  HINES GWD PARAMETRIZATION
    37523904
    37533905  IF (ok_hines) then
    3754 
    3755      CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
    3756           rlat,t_seri,u_seri,v_seri, &
    3757           zustrhi,zvstrhi, &
    3758           d_t_hin, d_u_hin, d_v_hin)
    3759      !
    3760      !  ajout des tendances
    3761      CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin',abortphy)
    3762 
     3906     !  HINES GWD PARAMETRIZATION
     3907     east_gwstress=0.
     3908     west_gwstress=0.
     3909     du_gwd_hines=0.
     3910     dv_gwd_hines=0.
     3911     CALL hines_gwd(klon, klev, dtime, paprs, pplay, rlat, t_seri, u_seri, &
     3912          v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, du_gwd_hines, &
     3913          dv_gwd_hines)
     3914     zustr_gwd_hines=0.
     3915     zvstr_gwd_hines=0.
     3916     DO k = 1, klev
     3917        zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
     3918             * (paprs(:, k)-paprs(:, k+1))/rg
     3919        zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
     3920             * (paprs(:, k)-paprs(:, k+1))/rg
     3921     ENDDO
     3922
     3923     d_t_hin(:, :)=0.
     3924     CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, dqi0, &
     3925          paprs, 'hin', abortphy)
     3926  ENDIF
     3927
     3928  IF (.not. ok_hines .and. ok_gwd_rando) then
     3929     CALL acama_GWD_rando(DTIME, pplay, rlat, t_seri, u_seri, v_seri, rot, &
     3930          zustr_gwd_front, zvstr_gwd_front, du_gwd_front, dv_gwd_front, &
     3931          east_gwstress, west_gwstress)
     3932     zustr_gwd_front=0.
     3933     zvstr_gwd_front=0.
     3934     DO k = 1, klev
     3935        zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
     3936             * (paprs(:, k)-paprs(:, k+1))/rg
     3937        zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
     3938             * (paprs(:, k)-paprs(:, k+1))/rg
     3939     ENDDO
     3940
     3941     CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
     3942          paprs, 'front_gwd_rando', abortphy)
    37633943  ENDIF
    37643944
    … …  
    37663946     call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
    37673947          rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
    3768           du_gwd_rando, dv_gwd_rando)
    3769      CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
    3770           'flott_gwd_rando',abortphy)
     3948          du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
     3949     CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
     3950          paprs, 'flott_gwd_rando', abortphy)
     3951     zustr_gwd_rando=0.
     3952     zvstr_gwd_rando=0.
     3953     DO k = 1, klev
     3954        zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
     3955             * (paprs(:, k)-paprs(:, k+1))/rg
     3956        zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
     3957             * (paprs(:, k)-paprs(:, k+1))/rg
     3958     ENDDO
    37713959  end if
    37723960
    … …  
    37743962
    37753963  if (mydebug) then
    3776      call writefield_phy('u_seri',u_seri,llm)
    3777      call writefield_phy('v_seri',v_seri,llm)
    3778      call writefield_phy('t_seri',t_seri,llm)
    3779      call writefield_phy('q_seri',q_seri,llm)
     3964     call writefield_phy('u_seri',u_seri,nbp_lev)
     3965     call writefield_phy('v_seri',v_seri,nbp_lev)
     3966     call writefield_phy('t_seri',t_seri,nbp_lev)
     3967     call writefield_phy('q_seri',q_seri,nbp_lev)
    37803968  endif
    37813969
    … …  
    38083996  IF (ip_ebil_phy.ge.2) THEN
    38093997     ztit='after orography'
    3810      CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
     3998     CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
    38113999          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    38124000          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    3813      call diagphy(airephy,ztit,ip_ebil_phy &
     4001     call diagphy(cell_area,ztit,ip_ebil_phy &
    38144002          , zero_v, zero_v, zero_v, zero_v, zero_v &
    38154003          , zero_v, zero_v, zero_v, ztsol &
    … …  
    38224010     CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
    38234011  ! ajout de la tendance d'humidite due au methane
    3824      CALL add_phys_tend(du0,dv0,dt0,d_q_ch4*dtime,dql0,'q_ch4',abortphy)
     4012     CALL add_phys_tend(du0, dv0, dt0, d_q_ch4*dtime, dql0, dqi0, paprs, &
     4013          'q_ch4', abortphy)
    38254014  END IF
    38264015  !
    … …  
    39154104          cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
    39164105          frac_impa, frac_nucl, &
    3917           pphis,airephy,dtime,itap, &
     4106          pphis,cell_area,dtime,itap, &
    39184107          qx(:,:,ivap),da,phi,mp,upwd,dnwd)
    39194108
    … …  
    39464135
    39474136  d_t_ec(:,:)=0.
    3948   forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
     4137  forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
    39494138  CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
    39504139       u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
    … …  
    39554144  IF (ip_ebil_phy.ge.1) THEN
    39564145     ztit='after physic'
    3957      CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
     4146     CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
    39584147          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
    39594148          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    … …  
    39634152     !     Donc la somme de ces 2 variations devrait etre nulle.
    39644153
    3965      call diagphy(airephy,ztit,ip_ebil_phy &
     4154     call diagphy(cell_area,ztit,ip_ebil_phy &
    39664155          , topsw, toplw, solsw, sollw, sens &
    39674156          , evap, rain_fall, snow_fall, ztsol &
    … …  
    39984187  include "calcul_STDlev.h"
    39994188  !
    4000   ! slp sea level pressure
    4001   slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
     4189  ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
     4190  CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
    40024191  !
    40034192  !cc prw = eau precipitable
    … …  
    40234212          paprs, &
    40244213          q_seri, &
    4025           airephy, &
     4214          cell_area, &
    40264215          pphi, &
    40274216          pphis, &
    … …  
    40424231  !
    40434232  if (mydebug) then
    4044      call writefield_phy('u_seri',u_seri,llm)
    4045      call writefield_phy('v_seri',v_seri,llm)
    4046      call writefield_phy('t_seri',t_seri,llm)
    4047      call writefield_phy('q_seri',q_seri,llm)
     4233     call writefield_phy('u_seri',u_seri,nbp_lev)
     4234     call writefield_phy('v_seri',v_seri,nbp_lev)
     4235     call writefield_phy('t_seri',t_seri,nbp_lev)
     4236     call writefield_phy('q_seri',q_seri,nbp_lev)
    40484237  endif
    40494238
    … …  
    41334322     enddo
    41344323
    4135      write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
     4324!jyg<     (En attendant de statuer sur le sort de d_t_oli)
     4325!jyg!     write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
     4326!jyg!     do k=1,klev
     4327!jyg!        write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
     4328!jyg!             d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
     4329!jyg!     enddo
     4330     write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
    41364331     do k=1,klev
    4137         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
     4332        write(lunout,*) d_t_vdf(igout,k), &
    41384333             d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
    41394334     enddo
     4335!>jyg
    41404336
    41414337     write(lunout,*) 'd_ps ',d_ps(igout)
    … …  
    42384434    IF (abortphy==1) THEN
    42394435       abort_message ='Plantage hgardfou'
    4240        CALL abort_gcm (modname,abort_message,1)
     4436       CALL abort_physic (modname,abort_message,1)
    42414437    ENDIF
    42424438
    … …  
    42484444  !====================================================================
    42494445  !
    4250 
    4251   !        -----------------------------------------------------------------
    4252   !        WSTATS: Saving statistics
    4253   !        -----------------------------------------------------------------
    4254   !        ("stats" stores and accumulates 8 key variables in file "stats.nc"
    4255   !        which can later be used to make the statistic files of the run:
    4256   !        "stats")          only possible in 3D runs !
    4257 
    4258 
    4259   IF (callstats) THEN
    4260 
    4261      call wstats(klon,o_psol%name,"Surface pressure","Pa" &
    4262           ,2,paprs(:,1))
    4263      call wstats(klon,o_tsol%name,"Surface temperature","K", &
    4264           2,zxtsol)
    4265      zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
    4266      call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
    4267           "kg/(s*m2)",2,zx_tmp_fi2d)
    4268      zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
    4269      call wstats(klon,o_plul%name,"Large-scale Precip", &
    4270           "kg/(s*m2)",2,zx_tmp_fi2d)
    4271      zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
    4272      call wstats(klon,o_pluc%name,"Convective Precip", &
    4273           "kg/(s*m2)",2,zx_tmp_fi2d)
    4274      call wstats(klon,o_sols%name,"Solar rad. at surf.", &
    4275           "W/m2",2,solsw)
    4276      call wstats(klon,o_soll%name,"IR rad. at surf.", &
    4277           "W/m2",2,sollw)
    4278      zx_tmp_fi2d(:) = topsw(:)-toplw(:)
    4279      call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
    4280           "W/m2",2,zx_tmp_fi2d)
    4281 
    4282 
    4283 
    4284      call wstats(klon,o_temp%name,"Air temperature","K", &
    4285           3,t_seri)
    4286      call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
    4287           3,u_seri)
    4288      call wstats(klon,o_vitv%name,"Meridional wind", &
    4289           "m.s-1",3,v_seri)
    4290      call wstats(klon,o_vitw%name,"Vertical wind", &
    4291           "m.s-1",3,omega)
    4292      call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
    4293           3,q_seri)
    4294 
    4295 
    4296 
    4297      IF(lafin) THEN
    4298         write (*,*) "Writing stats..."
    4299         call mkstats(ierr)
    4300      ENDIF
    4301 
    4302   ENDIF !if callstats
    43034446
    43044447  IF (lafin) THEN
Note: See TracChangeset for help on using the changeset viewer.