/[lmdze]/trunk/phylmd/physiq.f
ViewVC logotype

Diff of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/physiq.f90 revision 73 by guez, Fri Nov 15 17:48:30 2013 UTC trunk/Sources/phylmd/physiq.f revision 174 by guez, Wed Nov 25 20:14:19 2015 UTC
# Line 4  module physiq_m Line 4  module physiq_m
4    
5  contains  contains
6    
7    SUBROUTINE physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &    SUBROUTINE physiq(lafin, dayvrai, time, paprs, play, pphi, pphis, u, v, t, &
8         u, v, t, qx, omega, d_u, d_v, d_t, d_qx, d_ps, dudyn, PVteta)         qx, omega, d_u, d_v, d_t, d_qx)
9    
10      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28      ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11      ! (subversion revision 678)      ! (subversion revision 678)
12    
13      ! Author: Z.X. Li (LMD/CNRS) 1993      ! Author: Z. X. Li (LMD/CNRS) 1993
14    
15      ! This is the main procedure for the "physics" part of the program.      ! This is the main procedure for the "physics" part of the program.
16    
# Line 18  contains Line 18  contains
18      USE abort_gcm_m, ONLY: abort_gcm      USE abort_gcm_m, ONLY: abort_gcm
19      use aeropt_m, only: aeropt      use aeropt_m, only: aeropt
20      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
     USE calendar, ONLY: ymds2ju  
21      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
22      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &      USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &
23           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin           ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin
24      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &      USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &
25           ok_orodr, ok_orolf, soil_model           ok_orodr, ok_orolf
26      USE clmain_m, ONLY: clmain      USE clmain_m, ONLY: clmain
27      use clouds_gno_m, only: clouds_gno      use clouds_gno_m, only: clouds_gno
28      USE comgeomphy, ONLY: airephy, cuphy, cvphy      use comconst, only: dtphys
29        USE comgeomphy, ONLY: airephy
30      USE concvl_m, ONLY: concvl      USE concvl_m, ONLY: concvl
31      USE conf_gcm_m, ONLY: offline, raz_date      USE conf_gcm_m, ONLY: offline, raz_date, day_step, iphysiq
32      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
33      use conflx_m, only: conflx      use conflx_m, only: conflx
34      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
35      use diagcld2_m, only: diagcld2      use diagcld2_m, only: diagcld2
36      use diagetpq_m, only: diagetpq      use diagetpq_m, only: diagetpq
37      use diagphy_m, only: diagphy      use diagphy_m, only: diagphy
38      USE dimens_m, ONLY: iim, jjm, llm, nqmx      USE dimens_m, ONLY: llm, nqmx
39      USE dimphy, ONLY: klon, nbtr      USE dimphy, ONLY: klon
40      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
41      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
42        use dynetat0_m, only: day_ref, annee_ref
43      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
44      use fisrtilp_m, only: fisrtilp      use fisrtilp_m, only: fisrtilp
45      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
     USE histsync_m, ONLY: histsync  
     USE histwrite_m, ONLY: histwrite  
46      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
47           nbsrf           nbsrf
     USE ini_histhf_m, ONLY: ini_histhf  
     USE ini_histday_m, ONLY: ini_histday  
48      USE ini_histins_m, ONLY: ini_histins      USE ini_histins_m, ONLY: ini_histins
49        use netcdf95, only: NF95_CLOSE
50      use newmicro_m, only: newmicro      use newmicro_m, only: newmicro
51      USE oasis_m, ONLY: ok_oasis      USE orbite_m, ONLY: orbite
     USE orbite_m, ONLY: orbite, zenang  
52      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
53      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0, rlat, rlon
54      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
55        USE phyredem0_m, ONLY: phyredem0
56      USE phystokenc_m, ONLY: phystokenc      USE phystokenc_m, ONLY: phystokenc
57      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
58      USE qcheck_m, ONLY: qcheck      USE qcheck_m, ONLY: qcheck
59      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
60      use readsulfate_m, only: readsulfate      use readsulfate_m, only: readsulfate
61      use sugwd_m, only: sugwd      use readsulfate_preind_m, only: readsulfate_preind
62      USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt      use yoegwd, only: sugwd
63      USE temps, ONLY: annee_ref, day_ref, itau_phy      USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt
64        USE temps, ONLY: itau_phy
65        use transp_m, only: transp
66      use unit_nml_m, only: unit_nml      use unit_nml_m, only: unit_nml
67        USE ymds2ju_m, ONLY: ymds2ju
68      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
69        use zenang_m, only: zenang
70    
71      ! Arguments:      logical, intent(in):: lafin ! dernier passage
72    
73      REAL, intent(in):: rdayvrai      integer, intent(in):: dayvrai
74      ! (elapsed time since January 1st 0h of the starting year, in days)      ! current day number, based at value 1 on January 1st of annee_ref
75    
76      REAL, intent(in):: time ! heure de la journée en fraction de jour      REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
     REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)  
     logical, intent(in):: lafin ! dernier passage  
77    
78      REAL, intent(in):: paprs(klon, llm + 1)      REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
79      ! (pression pour chaque inter-couche, en Pa)      ! pression pour chaque inter-couche, en Pa
80    
81      REAL, intent(in):: play(klon, llm)      REAL, intent(in):: play(:, :) ! (klon, llm)
82      ! (input pression pour le mileu de chaque couche (en Pa))      ! pression pour le mileu de chaque couche (en Pa)
83    
84      REAL, intent(in):: pphi(klon, llm)      REAL, intent(in):: pphi(:, :) ! (klon, llm)
85      ! (input geopotentiel de chaque couche (g z) (reference sol))      ! géopotentiel de chaque couche (référence sol)
86    
87      REAL, intent(in):: pphis(klon) ! input geopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
88    
89      REAL, intent(in):: u(klon, llm)      REAL, intent(in):: u(:, :) ! (klon, llm)
90      ! vitesse dans la direction X (de O a E) en m/s      ! vitesse dans la direction X (de O a E) en m/s
91    
92      REAL, intent(in):: v(klon, llm) ! vitesse Y (de S a N) en m/s      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m/s
93      REAL, intent(in):: t(klon, llm) ! input temperature (K)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
94    
95      REAL, intent(in):: qx(klon, llm, nqmx)      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
96      ! (humidité spécifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
97    
98      REAL omega(klon, llm) ! input vitesse verticale en Pa/s      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
99      REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
100      REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m/s/s)      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
101      REAL, intent(out):: d_t(klon, llm) ! tendance physique de "t" (K/s)      REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K/s)
     REAL d_qx(klon, llm, nqmx) ! output tendance physique de "qx" (kg/kg/s)  
     REAL d_ps(klon) ! output tendance physique de la pression au sol  
102    
103      LOGICAL:: firstcal = .true.      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
104        ! tendance physique de "qx" (s-1)
105    
106      INTEGER nbteta      ! Local:
     PARAMETER(nbteta = 3)  
107    
108      REAL PVteta(klon, nbteta)      LOGICAL:: firstcal = .true.
     ! (output vorticite potentielle a des thetas constantes)  
109    
110      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
111      PARAMETER (ok_gust = .FALSE.)      PARAMETER (ok_gust = .FALSE.)
112    
113      LOGICAL check ! Verifier la conservation du modele en eau      LOGICAL, PARAMETER:: check = .FALSE.
114      PARAMETER (check = .FALSE.)      ! Verifier la conservation du modele en eau
115    
116      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
117      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
118    
     ! Parametres lies au coupleur OASIS:  
     INTEGER, SAVE:: npas, nexca  
     logical rnpb  
     parameter(rnpb = .true.)  
   
     character(len = 6):: ocean = 'force '  
     ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")  
   
119      ! "slab" ocean      ! "slab" ocean
120      REAL, save:: tslab(klon) ! temperature of ocean slab      REAL, save:: tslab(klon) ! temperature of ocean slab
121      REAL, save:: seaice(klon) ! glace de mer (kg/m2)      REAL, save:: seaice(klon) ! glace de mer (kg/m2)
122      REAL fluxo(klon) ! flux turbulents ocean-glace de mer      REAL fluxo(klon) ! flux turbulents ocean-glace de mer
123      REAL fluxg(klon) ! flux turbulents ocean-atmosphere      REAL fluxg(klon) ! flux turbulents ocean-atmosphere
124    
     ! Modele thermique du sol, a activer pour le cycle diurne:  
     logical:: ok_veget = .false. ! type de modele de vegetation utilise  
   
125      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.      logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
126      ! sorties journalieres, mensuelles et instantanees dans les      ! sorties journalieres, mensuelles et instantanees dans les
127      ! fichiers histday, histmth et histins      ! fichiers histday, histmth et histins
# Line 148  contains Line 134  contains
134      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
135      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
136    
137      INTEGER ivap ! indice de traceurs pour vapeur d'eau      INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
138      PARAMETER (ivap = 1)      INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
     INTEGER iliq ! indice de traceurs pour eau liquide  
     PARAMETER (iliq = 2)  
139    
140      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
141      LOGICAL, save:: ancien_ok      LOGICAL, save:: ancien_ok
# Line 161  contains Line 145  contains
145    
146      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
147    
     !IM Amip2 PV a theta constante  
   
     CHARACTER(LEN = 3) ctetaSTD(nbteta)  
     DATA ctetaSTD/'350', '380', '405'/  
     REAL rtetaSTD(nbteta)  
     DATA rtetaSTD/350., 380., 405./  
   
     !MI Amip2 PV a theta constante  
   
148      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
149      REAL swup0(klon, llm + 1), swup(klon, llm + 1)      REAL swup0(klon, llm + 1), swup(klon, llm + 1)
150      SAVE swdn0, swdn, swup0, swup      SAVE swdn0, swdn, swup0, swup
# Line 178  contains Line 153  contains
153      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)      REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
154      SAVE lwdn0, lwdn, lwup0, lwup      SAVE lwdn0, lwdn, lwup0, lwup
155    
156      !IM Amip2      ! Amip2
157      ! variables a une pression donnee      ! variables a une pression donnee
158    
159      integer nlevSTD      integer nlevSTD
160      PARAMETER(nlevSTD = 17)      PARAMETER(nlevSTD = 17)
     real rlevSTD(nlevSTD)  
     DATA rlevSTD/100000., 92500., 85000., 70000., &  
          60000., 50000., 40000., 30000., 25000., 20000., &  
          15000., 10000., 7000., 5000., 3000., 2000., 1000./  
     CHARACTER(LEN = 4) clevSTD(nlevSTD)  
     DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &  
          '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &  
          '70 ', '50 ', '30 ', '20 ', '10 '/  
161    
162      ! prw: precipitable water      ! prw: precipitable water
163      real prw(klon)      real prw(klon)
# Line 203  contains Line 170  contains
170      INTEGER kmax, lmax      INTEGER kmax, lmax
171      PARAMETER(kmax = 8, lmax = 8)      PARAMETER(kmax = 8, lmax = 8)
172      INTEGER kmaxm1, lmaxm1      INTEGER kmaxm1, lmaxm1
173      PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)      PARAMETER(kmaxm1 = kmax - 1, lmaxm1 = lmax - 1)
   
     REAL zx_tau(kmaxm1), zx_pc(lmaxm1)  
     DATA zx_tau/0., 0.3, 1.3, 3.6, 9.4, 23., 60./  
     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./  
   
     ! cldtopres pression au sommet des nuages  
     REAL cldtopres(lmaxm1)  
     DATA cldtopres/50., 180., 310., 440., 560., 680., 800./  
   
     ! taulev: numero du niveau de tau dans les sorties ISCCP  
     CHARACTER(LEN = 4) taulev(kmaxm1)  
   
     DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/  
     CHARACTER(LEN = 3) pclev(lmaxm1)  
     DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/  
   
     CHARACTER(LEN = 28) cnameisccp(lmaxm1, kmaxm1)  
     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &  
          'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &  
          'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &  
          'pc= 680-800hPa, tau< 0.3', 'pc< 50hPa, tau= 0.3-1.3', &  
          'pc= 50-180hPa, tau= 0.3-1.3', 'pc= 180-310hPa, tau= 0.3-1.3', &  
          'pc= 310-440hPa, tau= 0.3-1.3', 'pc= 440-560hPa, tau= 0.3-1.3', &  
          'pc= 560-680hPa, tau= 0.3-1.3', 'pc= 680-800hPa, tau= 0.3-1.3', &  
          'pc< 50hPa, tau= 1.3-3.6', 'pc= 50-180hPa, tau= 1.3-3.6', &  
          'pc= 180-310hPa, tau= 1.3-3.6', 'pc= 310-440hPa, tau= 1.3-3.6', &  
          'pc= 440-560hPa, tau= 1.3-3.6', 'pc= 560-680hPa, tau= 1.3-3.6', &  
          'pc= 680-800hPa, tau= 1.3-3.6', 'pc< 50hPa, tau= 3.6-9.4', &  
          'pc= 50-180hPa, tau= 3.6-9.4', 'pc= 180-310hPa, tau= 3.6-9.4', &  
          'pc= 310-440hPa, tau= 3.6-9.4', 'pc= 440-560hPa, tau= 3.6-9.4', &  
          'pc= 560-680hPa, tau= 3.6-9.4', 'pc= 680-800hPa, tau= 3.6-9.4', &  
          'pc< 50hPa, tau= 9.4-23', 'pc= 50-180hPa, tau= 9.4-23', &  
          'pc= 180-310hPa, tau= 9.4-23', 'pc= 310-440hPa, tau= 9.4-23', &  
          'pc= 440-560hPa, tau= 9.4-23', 'pc= 560-680hPa, tau= 9.4-23', &  
          'pc= 680-800hPa, tau= 9.4-23', 'pc< 50hPa, tau= 23-60', &  
          'pc= 50-180hPa, tau= 23-60', 'pc= 180-310hPa, tau= 23-60', &  
          'pc= 310-440hPa, tau= 23-60', 'pc= 440-560hPa, tau= 23-60', &  
          'pc= 560-680hPa, tau= 23-60', 'pc= 680-800hPa, tau= 23-60', &  
          'pc< 50hPa, tau> 60.', 'pc= 50-180hPa, tau> 60.', &  
          'pc= 180-310hPa, tau> 60.', 'pc= 310-440hPa, tau> 60.', &  
          'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &  
          'pc= 680-800hPa, tau> 60.'/  
   
     !IM ISCCP simulator v3.4  
   
     integer nid_hf, nid_hf3d  
     save nid_hf, nid_hf3d  
174    
175      ! Variables propres a la physique      ! Variables propres a la physique
176    
177      INTEGER, save:: radpas      INTEGER, save:: radpas
178      ! (Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
179      ! "physiq".)      ! "physiq".
180    
181      REAL radsol(klon)      REAL radsol(klon)
182      SAVE radsol ! bilan radiatif au sol calcule par code radiatif      SAVE radsol ! bilan radiatif au sol calcule par code radiatif
183    
184      INTEGER, SAVE:: itap ! number of calls to "physiq"      INTEGER:: itap = 0 ! number of calls to "physiq"
185    
186      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
187    
# Line 272  contains Line 192  contains
192      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
193      SAVE fluxlat      SAVE fluxlat
194    
195      REAL fqsurf(klon, nbsrf)      REAL, save:: fqsurf(klon, nbsrf)
196      SAVE fqsurf ! humidite de l'air au contact de la surface      ! humidite de l'air au contact de la surface
197    
198      REAL, save:: qsol(klon) ! hauteur d'eau dans le sol      REAL, save:: qsol(klon)
199        ! column-density of water in soil, in kg m-2
200    
201      REAL fsnow(klon, nbsrf)      REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
202      SAVE fsnow ! epaisseur neigeuse      REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
203    
204      REAL falbe(klon, nbsrf)      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
     SAVE falbe ! albedo par type de surface  
     REAL falblw(klon, nbsrf)  
     SAVE falblw ! albedo par type de surface  
   
     ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :  
205      REAL, save:: zmea(klon) ! orographie moyenne      REAL, save:: zmea(klon) ! orographie moyenne
206      REAL, save:: zstd(klon) ! deviation standard de l'OESM      REAL, save:: zstd(klon) ! deviation standard de l'OESM
207      REAL, save:: zsig(klon) ! pente de l'OESM      REAL, save:: zsig(klon) ! pente de l'OESM
# Line 307  contains Line 223  contains
223      !KE43      !KE43
224      ! Variables liees a la convection de K. Emanuel (sb):      ! Variables liees a la convection de K. Emanuel (sb):
225    
     REAL bas, top ! cloud base and top levels  
     SAVE bas  
     SAVE top  
   
226      REAL Ma(klon, llm) ! undilute upward mass flux      REAL Ma(klon, llm) ! undilute upward mass flux
227      SAVE Ma      SAVE Ma
228      REAL qcondc(klon, llm) ! in-cld water content from convect      REAL qcondc(klon, llm) ! in-cld water content from convect
# Line 344  contains Line 256  contains
256      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)      REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
257      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
258    
259      REAL, save:: rain_fall(klon) ! pluie      REAL, save:: rain_fall(klon)
260      REAL, save:: snow_fall(klon) ! neige      ! liquid water mass flux (kg/m2/s), positive down
261    
262        REAL, save:: snow_fall(klon)
263        ! solid water mass flux (kg/m2/s), positive down
264    
265      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
266    
# Line 354  contains Line 269  contains
269      REAL dlw(klon) ! derivee infra rouge      REAL dlw(klon) ! derivee infra rouge
270      SAVE dlw      SAVE dlw
271      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
272      REAL fder(klon) ! Derive de flux (sensible et latente)      REAL, save:: fder(klon) ! Derive de flux (sensible et latente)
     save fder  
273      REAL ve(klon) ! integr. verticale du transport meri. de l'energie      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
274      REAL vq(klon) ! integr. verticale du transport meri. de l'eau      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
275      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
276      REAL uq(klon) ! integr. verticale du transport zonal de l'eau      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
277    
278      REAL frugs(klon, nbsrf) ! longueur de rugosite      REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
     save frugs  
279      REAL zxrugs(klon) ! longueur de rugosite      REAL zxrugs(klon) ! longueur de rugosite
280    
281      ! Conditions aux limites      ! Conditions aux limites
282    
283      INTEGER julien      INTEGER julien
   
284      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day      INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
285      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
286      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
287        REAL, save:: albsol(klon) ! albedo du sol total visible
     REAL albsol(klon)  
     SAVE albsol ! albedo du sol total  
     REAL albsollw(klon)  
     SAVE albsollw ! albedo du sol total  
   
288      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
289    
290      ! Declaration des procedures appelees      ! Declaration des procedures appelees
291    
     EXTERNAL alboc ! calculer l'albedo sur ocean  
     !KE43  
     EXTERNAL conema3 ! convect4.3  
292      EXTERNAL nuage ! calculer les proprietes radiatives      EXTERNAL nuage ! calculer les proprietes radiatives
     EXTERNAL transp ! transport total de l'eau et de l'energie  
293    
294      ! Variables locales      ! Variables locales
295    
# Line 411  contains Line 314  contains
314      REAL zxfluxu(klon, llm)      REAL zxfluxu(klon, llm)
315      REAL zxfluxv(klon, llm)      REAL zxfluxv(klon, llm)
316    
317      ! Le rayonnement n'est pas calculé tous les pas, il faut donc que      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
318      ! les variables soient rémanentes.      ! les variables soient r\'emanentes.
319      REAL, save:: heat(klon, llm) ! chauffage solaire      REAL, save:: heat(klon, llm) ! chauffage solaire
320      REAL heat0(klon, llm) ! chauffage solaire ciel clair      REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
321      REAL, save:: cool(klon, llm) ! refroidissement infrarouge      REAL, save:: cool(klon, llm) ! refroidissement infrarouge
322      REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair      REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
323      REAL, save:: topsw(klon), toplw(klon), solsw(klon)      REAL, save:: topsw(klon), toplw(klon), solsw(klon)
324      REAL, save:: sollw(klon) ! rayonnement infrarouge montant à la surface      REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
325      real, save:: sollwdown(klon) ! downward LW flux at surface      real, save:: sollwdown(klon) ! downward LW flux at surface
326      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
327      REAL albpla(klon)      REAL, save:: albpla(klon)
328      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
329      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
     SAVE albpla  
     SAVE heat0, cool0  
   
     INTEGER itaprad  
     SAVE itaprad  
330    
331      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
332      REAL conv_t(klon, llm) ! convergence of temperature (K/s)      REAL conv_t(klon, llm) ! convergence of temperature (K/s)
# Line 438  contains Line 336  contains
336    
337      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
338    
339      REAL dist, rmu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
340      REAL zdtime ! pas de temps du rayonnement (s)      real longi
     real zlongi  
341      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
342      REAL za, zb      REAL za, zb
343      REAL zx_t, zx_qs, zdelta, zcor      REAL zx_t, zx_qs, zcor
344      real zqsat(klon, llm)      real zqsat(klon, llm)
345      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
346      REAL, PARAMETER:: t_coup = 234.      REAL, PARAMETER:: t_coup = 234.
347      REAL zphi(klon, llm)      REAL zphi(klon, llm)
348    
349      !IM cf. AM Variables locales pour la CLA (hbtm2)      ! cf. AM Variables locales pour la CLA (hbtm2)
350    
351      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
352      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
# Line 472  contains Line 369  contains
369      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
370      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
371      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
     REAL tvp(klon, llm) ! virtual temp of lifted parcel  
372      REAL cape(klon) ! CAPE      REAL cape(klon) ! CAPE
373      SAVE cape      SAVE cape
374    
     REAL pbase(klon) ! cloud base pressure  
     SAVE pbase  
     REAL bbase(klon) ! cloud base buoyancy  
     SAVE bbase  
     REAL rflag(klon) ! flag fonctionnement de convect  
375      INTEGER iflagctrl(klon) ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
     ! -- convect43:  
     REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)  
     REAL dplcldt(klon), dplcldr(klon)  
376    
377      ! Variables du changement      ! Variables du changement
378    
379      ! con: convection      ! con: convection
380      ! lsc: large scale condensation      ! lsc: large scale condensation
381      ! ajs: ajustement sec      ! ajs: ajustement sec
382      ! eva: évaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
383      ! vdf: vertical diffusion in boundary layer      ! vdf: vertical diffusion in boundary layer
384      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
385      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL d_u_con(klon, llm), d_v_con(klon, llm)
# Line 534  contains Line 422  contains
422      integer:: iflag_cldcon = 1      integer:: iflag_cldcon = 1
423      logical ptconv(klon, llm)      logical ptconv(klon, llm)
424    
425      ! Variables locales pour effectuer les appels en série :      ! Variables locales pour effectuer les appels en s\'erie :
426    
427      REAL t_seri(klon, llm), q_seri(klon, llm)      REAL t_seri(klon, llm), q_seri(klon, llm)
428      REAL ql_seri(klon, llm), qs_seri(klon, llm)      REAL ql_seri(klon, llm)
429      REAL u_seri(klon, llm), v_seri(klon, llm)      REAL u_seri(klon, llm), v_seri(klon, llm)
430        REAL tr_seri(klon, llm, nqmx - 2)
     REAL tr_seri(klon, llm, nbtr)  
     REAL d_tr(klon, llm, nbtr)  
431    
432      REAL zx_rh(klon, llm)      REAL zx_rh(klon, llm)
433    
# Line 550  contains Line 436  contains
436      REAL zustrph(klon), zvstrph(klon)      REAL zustrph(klon), zvstrph(klon)
437      REAL aam, torsfc      REAL aam, torsfc
438    
     REAL dudyn(iim + 1, jjm + 1, llm)  
   
439      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique      REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
     REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)  
440    
441      INTEGER, SAVE:: nid_day, nid_ins      INTEGER, SAVE:: nid_ins
442    
443      REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.      REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
444      REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.      REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
445      REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.      REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
446      REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.      REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
447    
     REAL zsto  
   
     logical ok_sync  
448      real date0      real date0
449    
450      ! Variables liées au bilan d'énergie et d'enthalpie :      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
451      REAL ztsol(klon)      REAL ztsol(klon)
452      REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec      REAL d_h_vcol, d_qt, d_ec
453      REAL, SAVE:: d_h_vcol_phy      REAL, SAVE:: d_h_vcol_phy
     REAL fs_bound, fq_bound  
454      REAL zero_v(klon)      REAL zero_v(klon)
455      CHARACTER(LEN = 15) tit      CHARACTER(LEN = 20) tit
456      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
457      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
458    
459      REAL d_t_ec(klon, llm) ! tendance due à la conversion Ec -> E thermique      REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
460      REAL ZRCPD      REAL ZRCPD
461    
462      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m      REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
# Line 634  contains Line 513  contains
513      ! (column-density of mass of air in a cell, in kg m-2)      ! (column-density of mass of air in a cell, in kg m-2)
514    
515      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
516        integer, save:: ncid_startphy
517    
518      namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &      namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &
519           fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &           facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
520           ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, &           ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, nsplit_thermals
          nsplit_thermals  
521    
522      !----------------------------------------------------------------      !----------------------------------------------------------------
523    
524      IF (if_ebil >= 1) zero_v = 0.      IF (if_ebil >= 1) zero_v = 0.
     ok_sync = .TRUE.  
525      IF (nqmx < 2) CALL abort_gcm('physiq', &      IF (nqmx < 2) CALL abort_gcm('physiq', &
526           'eaux vapeur et liquide sont indispensables', 1)           'eaux vapeur et liquide sont indispensables')
527    
528      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
529         ! initialiser         ! initialiser
# Line 658  contains Line 536  contains
536         piz_ae = 0.         piz_ae = 0.
537         tau_ae = 0.         tau_ae = 0.
538         cg_ae = 0.         cg_ae = 0.
539         rain_con(:) = 0.         rain_con = 0.
540         snow_con(:) = 0.         snow_con = 0.
541         topswai(:) = 0.         topswai = 0.
542         topswad(:) = 0.         topswad = 0.
543         solswai(:) = 0.         solswai = 0.
544         solswad(:) = 0.         solswad = 0.
545    
546         d_u_con = 0.         d_u_con = 0.
547         d_v_con = 0.         d_v_con = 0.
# Line 696  contains Line 574  contains
574         ! Initialiser les compteurs:         ! Initialiser les compteurs:
575    
576         frugs = 0.         frugs = 0.
577         itap = 0         CALL phyetat0(pctsrf, ftsol, ftsoil, tslab, seaice, fqsurf, qsol, &
578         itaprad = 0              fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
579         CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &              radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
580              seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &              t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
581              snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, zmea, &              run_off_lic_0, sig1, w01, ncid_startphy)
             zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &  
             ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)  
582    
583         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
584         q2 = 1e-8         q2 = 1e-8
585    
586         radpas = NINT(86400. / dtphys / nbapp_rad)         lmt_pas = day_step / iphysiq
587           print *, 'Number of time steps of "physics" per day: ', lmt_pas
588    
589           radpas = lmt_pas / nbapp_rad
590    
591         ! on remet le calendrier a zero         ! On remet le calendrier a zero
592         IF (raz_date) itau_phy = 0         IF (raz_date) itau_phy = 0
593    
594         PRINT *, 'cycle_diurne = ', cycle_diurne         CALL printflag(radpas, ok_journe, ok_instan, ok_region)
        CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &  
             ok_instan, ok_region)  
   
        IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN  
           print *, "Au minimum 4 appels par jour si cycle diurne"  
           call abort_gcm('physiq', &  
                "Nombre d'appels au rayonnement insuffisant", 1)  
        ENDIF  
595    
596         ! Initialisation pour le schéma de convection d'Emanuel :         ! Initialisation pour le sch\'ema de convection d'Emanuel :
597         IF (iflag_con >= 3) THEN         IF (iflag_con >= 3) THEN
598            ibas_con = 1            ibas_con = 1
599            itop_con = 1            itop_con = 1
# Line 735  contains Line 606  contains
606            rugoro = 0.            rugoro = 0.
607         ENDIF         ENDIF
608    
        lmt_pas = NINT(86400. / dtphys) ! tous les jours  
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
609         ecrit_ins = NINT(ecrit_ins/dtphys)         ecrit_ins = NINT(ecrit_ins/dtphys)
610         ecrit_hf = NINT(ecrit_hf/dtphys)         ecrit_hf = NINT(ecrit_hf/dtphys)
611         ecrit_mth = NINT(ecrit_mth/dtphys)         ecrit_mth = NINT(ecrit_mth/dtphys)
612         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)         ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
613         ecrit_reg = NINT(ecrit_reg/dtphys)         ecrit_reg = NINT(ecrit_reg/dtphys)
614    
        ! Initialiser le couplage si necessaire  
   
        npas = 0  
        nexca = 0  
   
615         ! Initialisation des sorties         ! Initialisation des sorties
616    
        call ini_histhf(dtphys, nid_hf, nid_hf3d)  
        call ini_histday(dtphys, ok_journe, nid_day, nqmx)  
617         call ini_histins(dtphys, ok_instan, nid_ins)         call ini_histins(dtphys, ok_instan, nid_ins)
618         CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
619         ! Positionner date0 pour initialisation de ORCHIDEE         ! Positionner date0 pour initialisation de ORCHIDEE
620         print *, 'physiq date0: ', date0         print *, 'physiq date0: ', date0
621           CALL phyredem0(lmt_pas)
622      ENDIF test_firstcal      ENDIF test_firstcal
623    
624      ! Mettre a zero des variables de sortie (pour securite)      ! We will modify variables *_seri and we will not touch variables
625        ! u, v, t, qx:
626      DO i = 1, klon      t_seri = t
627         d_ps(i) = 0.      u_seri = u
628      ENDDO      v_seri = v
629      DO iq = 1, nqmx      q_seri = qx(:, :, ivap)
630         DO k = 1, llm      ql_seri = qx(:, :, iliq)
631            DO i = 1, klon      tr_seri = qx(:, :, 3:nqmx)
              d_qx(i, k, iq) = 0.  
           ENDDO  
        ENDDO  
     ENDDO  
     da = 0.  
     mp = 0.  
     phi = 0.  
   
     ! Ne pas affecter les valeurs entrées de u, v, h, et q :  
632    
633      DO k = 1, llm      ztsol = sum(ftsol * pctsrf, dim = 2)
        DO i = 1, klon  
           t_seri(i, k) = t(i, k)  
           u_seri(i, k) = u(i, k)  
           v_seri(i, k) = v(i, k)  
           q_seri(i, k) = qx(i, k, ivap)  
           ql_seri(i, k) = qx(i, k, iliq)  
           qs_seri(i, k) = 0.  
        ENDDO  
     ENDDO  
     IF (nqmx >= 3) THEN  
        tr_seri(:, :, :nqmx-2) = qx(:, :, 3:nqmx)  
     ELSE  
        tr_seri(:, :, 1) = 0.  
     ENDIF  
   
     DO i = 1, klon  
        ztsol(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
634    
635      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
636         tit = 'after dynamics'         tit = 'after dynamics'
637         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
638              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
639              d_ql, d_qs, d_ec)         ! Comme les tendances de la physique sont ajout\'es dans la
        ! Comme les tendances de la physique sont ajoutés dans la  
640         !  dynamique, la variation d'enthalpie par la dynamique devrait         !  dynamique, la variation d'enthalpie par la dynamique devrait
641         !  être égale à la variation de la physique au pas de temps         !  \^etre \'egale \`a la variation de la physique au pas de temps
642         !  précédent.  Donc la somme de ces 2 variations devrait être         !  pr\'ec\'edent.  Donc la somme de ces 2 variations devrait \^etre
643         !  nulle.         !  nulle.
644         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
645              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
646              d_qt, 0., fs_bound, fq_bound)              d_qt, 0.)
647      END IF      END IF
648    
649      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
# Line 845  contains Line 674  contains
674      ! Check temperatures:      ! Check temperatures:
675      CALL hgardfou(t_seri, ftsol)      CALL hgardfou(t_seri, ftsol)
676    
677      ! Incrementer le compteur de la physique      ! Incrémenter le compteur de la physique
678      itap = itap + 1      itap = itap + 1
679      julien = MOD(NINT(rdayvrai), 360)      julien = MOD(dayvrai, 360)
680      if (julien == 0) julien = 360      if (julien == 0) julien = 360
681    
682      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
   
     ! Mettre en action les conditions aux limites (albedo, sst etc.).  
683    
684      ! Prescrire l'ozone et calculer l'albedo sur l'ocean.      ! Prescrire l'ozone :
685      wo = ozonecm(REAL(julien), paprs)      wo = ozonecm(REAL(julien), paprs)
686    
687      ! Évaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
688      DO k = 1, llm      DO k = 1, llm
689         DO i = 1, klon         DO i = 1, klon
690            zb = MAX(0., ql_seri(i, k))            zb = MAX(0., ql_seri(i, k))
# Line 871  contains Line 698  contains
698      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
699         tit = 'after reevap'         tit = 'after reevap'
700         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
701              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
702         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
703              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
   
704      END IF      END IF
705    
706      ! Appeler la diffusion verticale (programme de couche limite)      frugs = MAX(frugs, 0.000015)
707        zxrugs = sum(frugs * pctsrf, dim = 2)
708    
709      DO i = 1, klon      ! Calculs nécessaires au calcul de l'albedo dans l'interface avec
710         zxrugs(i) = 0.      ! la surface.
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)  
        ENDDO  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
   
     ! calculs necessaires au calcul de l'albedo dans l'interface  
711    
712      CALL orbite(REAL(julien), zlongi, dist)      CALL orbite(REAL(julien), longi, dist)
713      IF (cycle_diurne) THEN      IF (cycle_diurne) THEN
714         zdtime = dtphys * REAL(radpas)         CALL zenang(longi, time, dtphys * radpas, mu0, fract)
        CALL zenang(zlongi, time, zdtime, rmu0, fract)  
715      ELSE      ELSE
716         rmu0 = -999.999         mu0 = - 999.999
717      ENDIF      ENDIF
718    
719      ! Calcul de l'abedo moyen par maille      ! Calcul de l'abedo moyen par maille
720      albsol(:) = 0.      albsol = sum(falbe * pctsrf, dim = 2)
     albsollw(:) = 0.  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)  
           albsollw(i) = albsollw(i) + falblw(i, nsrf) * pctsrf(i, nsrf)  
        ENDDO  
     ENDDO  
721    
722      ! Répartition sous maille des flux longwave et shortwave      ! R\'epartition sous maille des flux longwave et shortwave
723      ! Répartition du longwave par sous-surface linéarisée      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
724    
725      DO nsrf = 1, nbsrf      forall (nsrf = 1: nbsrf)
726         DO i = 1, klon         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
727            fsollw(i, nsrf) = sollw(i) &              * (ztsol - ftsol(:, nsrf))
728                 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
729            fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))      END forall
        ENDDO  
     ENDDO  
730    
731      fder = dlw      fder = dlw
732    
733      ! Couche limite:      ! Couche limite:
734    
735      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, &      CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
736           u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, &           v_seri, julien, mu0, co2_ppm, ftsol, cdmmax, cdhmax, ksta, ksta_ter, &
737           ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &           ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, &
738           qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &           fluxlat, rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, &
739           rain_fall, snow_fall, fsolsw, fsollw, fder, rlon, rlat, &           firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, &
740           frugs, firstcal, agesno, rugoro, d_t_vdf, &           fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, &
741           d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &           ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, &
742           cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &           pblT, therm, trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, &
743           pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &           run_off_lic_0, fluxo, fluxg, tslab)
          fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)  
744    
745      ! Incrémentation des flux      ! Incr\'ementation des flux
746    
747      zxfluxt = 0.      zxfluxt = 0.
748      zxfluxq = 0.      zxfluxq = 0.
# Line 959  contains Line 760  contains
760      END DO      END DO
761      DO i = 1, klon      DO i = 1, klon
762         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol         sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
763         evap(i) = - zxfluxq(i, 1) ! flux d'évaporation au sol         evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
764         fder(i) = dlw(i) + dsens(i) + devap(i)         fder(i) = dlw(i) + dsens(i) + devap(i)
765      ENDDO      ENDDO
766    
# Line 975  contains Line 776  contains
776      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
777         tit = 'after clmain'         tit = 'after clmain'
778         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
779              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
780         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
781              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
782      END IF      END IF
783    
784      ! Update surface temperature:      ! Update surface temperature:
# Line 1008  contains Line 807  contains
807    
808         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &         IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
809              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &              + pctsrf(i, is_sic) - 1.)  >  EPSFRA) print *, &
810              'physiq : problème sous surface au point ', i, pctsrf(i, 1 : nbsrf)              'physiq : probl\`eme sous surface au point ', i, &
811                pctsrf(i, 1 : nbsrf)
812      ENDDO      ENDDO
813      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
814         DO i = 1, klon         DO i = 1, klon
# Line 1036  contains Line 836  contains
836         ENDDO         ENDDO
837      ENDDO      ENDDO
838    
839      ! Si une sous-fraction n'existe pas, elle prend la temp. moyenne      ! Si une sous-fraction n'existe pas, elle prend la température moyenne :
   
840      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
841         DO i = 1, klon         DO i = 1, klon
842            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)            IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
# Line 1062  contains Line 861  contains
861         ENDDO         ENDDO
862      ENDDO      ENDDO
863    
864      ! Calculer la derive du flux infrarouge      ! Calculer la dérive du flux infrarouge
865    
866      DO i = 1, klon      DO i = 1, klon
867         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3         dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
868      ENDDO      ENDDO
869    
870      ! Appeler la convection (au choix)      IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
871    
872      DO k = 1, llm      ! Appeler la convection (au choix)
        DO i = 1, klon  
           conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys  
           conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys  
        ENDDO  
     ENDDO  
   
     IF (check) THEN  
        za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)  
        print *, "avantcon = ", za  
     ENDIF  
873    
874      if (iflag_con == 2) then      if (iflag_con == 2) then
875           conv_q = d_q_dyn + d_q_vdf / dtphys
876           conv_t = d_t_dyn + d_t_vdf / dtphys
877         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
878         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
879              q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &              q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
880              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
881              mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &              mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
882              kdtop, pmflxr, pmflxs)              kdtop, pmflxr, pmflxs)
883         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
884         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
# Line 1096  contains Line 887  contains
887      else      else
888         ! iflag_con >= 3         ! iflag_con >= 3
889    
890         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &         da = 0.
891              v_seri, tr_seri, sig1, w01, d_t_con, d_q_con, &         mp = 0.
892              d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &         phi = 0.
893              itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &         CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
894              pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &              w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, snow_con, &
895              wd, pmflxr, pmflxs, da, phi, mp, ntra=1)              ibas_con, itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, &
896         ! (number of tracers for the convection scheme of Kerry Emanuel:              qcondc, wd, pmflxr, pmflxs, da, phi, mp)
        ! la partie traceurs est faite dans phytrac  
        ! on met ntra = 1 pour limiter les appels mais on peut  
        ! supprimer les calculs / ftra.)  
   
897         clwcon0 = qcondc         clwcon0 = qcondc
898         mfu = upwd + dnwd         mfu = upwd + dnwd
899         IF (.NOT. ok_gust) wd = 0.         IF (.NOT. ok_gust) wd = 0.
900    
901         ! Calcul des propriétés des nuages convectifs         IF (thermcep) THEN
902              zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
903         DO k = 1, llm            zqsat = zqsat / (1. - retv * zqsat)
904            DO i = 1, klon         ELSE
905               zx_t = t_seri(i, k)            zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
906               IF (thermcep) THEN         ENDIF
                 zdelta = MAX(0., SIGN(1., rtt-zx_t))  
                 zx_qs = r2es * FOEEW(zx_t, zdelta) / play(i, k)  
                 zx_qs = MIN(0.5, zx_qs)  
                 zcor = 1./(1.-retv*zx_qs)  
                 zx_qs = zx_qs*zcor  
              ELSE  
                 IF (zx_t < t_coup) THEN  
                    zx_qs = qsats(zx_t)/play(i, k)  
                 ELSE  
                    zx_qs = qsatl(zx_t)/play(i, k)  
                 ENDIF  
              ENDIF  
              zqsat(i, k) = zx_qs  
           ENDDO  
        ENDDO  
907    
908         ! calcul des proprietes des nuages convectifs         ! Properties of convective clouds
909         clwcon0 = fact_cldcon * clwcon0         clwcon0 = fact_cldcon * clwcon0
910         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &         call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
911              rnebcon0)              rnebcon0)
# Line 1157  contains Line 929  contains
929      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
930         tit = 'after convect'         tit = 'after convect'
931         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
932              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
933         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
934              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
935      END IF      END IF
936    
937      IF (check) THEN      IF (check) THEN
938         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
939         print *, "aprescon = ", za         print *, "aprescon = ", za
940         zx_t = 0.         zx_t = 0.
941         za = 0.         za = 0.
# Line 1190  contains Line 960  contains
960         ENDDO         ENDDO
961      ENDIF      ENDIF
962    
963      ! Convection sèche (thermiques ou ajustement)      ! Convection s\`eche (thermiques ou ajustement)
964    
965      d_t_ajs = 0.      d_t_ajs = 0.
966      d_u_ajs = 0.      d_u_ajs = 0.
# Line 1213  contains Line 983  contains
983      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
984         tit = 'after dry_adjust'         tit = 'after dry_adjust'
985         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
986              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
987      END IF      END IF
988    
989      ! Caclul des ratqs      ! Caclul des ratqs
990    
991      ! ratqs convectifs à l'ancienne en fonction de (q(z = 0) - q) / q      ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
992      ! on écrase le tableau ratqsc calculé par clouds_gno      ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
993      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
994         do k = 1, llm         do k = 1, llm
995            do i = 1, klon            do i = 1, klon
# Line 1272  contains Line 1041  contains
1041         ENDDO         ENDDO
1042      ENDDO      ENDDO
1043      IF (check) THEN      IF (check) THEN
1044         za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)         za = qcheck(paprs, q_seri, ql_seri)
1045         print *, "apresilp = ", za         print *, "apresilp = ", za
1046         zx_t = 0.         zx_t = 0.
1047         za = 0.         za = 0.
# Line 1288  contains Line 1057  contains
1057      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1058         tit = 'after fisrt'         tit = 'after fisrt'
1059         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1060              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1061         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &         call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1062              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
1063      END IF      END IF
1064    
1065      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1066    
1067      ! 1. NUAGES CONVECTIFS      ! 1. NUAGES CONVECTIFS
1068    
1069      IF (iflag_cldcon <= -1) THEN      IF (iflag_cldcon <= - 1) THEN
1070         ! seulement pour Tiedtke         ! seulement pour Tiedtke
1071         snow_tiedtke = 0.         snow_tiedtke = 0.
1072         if (iflag_cldcon == -1) then         if (iflag_cldcon == - 1) then
1073            rain_tiedtke = rain_con            rain_tiedtke = rain_con
1074         else         else
1075            rain_tiedtke = 0.            rain_tiedtke = 0.
1076            do k = 1, llm            do k = 1, llm
1077               do i = 1, klon               do i = 1, klon
1078                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
1079                     rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1080                          *zmasse(i, k)                          *zmasse(i, k)
1081                  endif                  endif
1082               enddo               enddo
# Line 1329  contains Line 1096  contains
1096         ENDDO         ENDDO
1097      ELSE IF (iflag_cldcon == 3) THEN      ELSE IF (iflag_cldcon == 3) THEN
1098         ! On prend pour les nuages convectifs le maximum du calcul de         ! On prend pour les nuages convectifs le maximum du calcul de
1099         ! la convection et du calcul du pas de temps précédent diminué         ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1100         ! d'un facteur facttemps.         ! d'un facteur facttemps.
1101         facteur = dtphys * facttemps         facteur = dtphys * facttemps
1102         do k = 1, llm         do k = 1, llm
# Line 1369  contains Line 1136  contains
1136      ENDDO      ENDDO
1137    
1138      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1139           dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &           dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1140           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1141    
1142      ! Humidité relative pour diagnostic :      ! Humidit\'e relative pour diagnostic :
1143      DO k = 1, llm      DO k = 1, llm
1144         DO i = 1, klon         DO i = 1, klon
1145            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
1146            IF (thermcep) THEN            IF (thermcep) THEN
1147               zdelta = MAX(0., SIGN(1., rtt-zx_t))               zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
              zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)  
1148               zx_qs = MIN(0.5, zx_qs)               zx_qs = MIN(0.5, zx_qs)
1149               zcor = 1./(1.-retv*zx_qs)               zcor = 1./(1. - retv*zx_qs)
1150               zx_qs = zx_qs*zcor               zx_qs = zx_qs*zcor
1151            ELSE            ELSE
1152               IF (zx_t < t_coup) THEN               IF (zx_t < t_coup) THEN
# Line 1397  contains Line 1163  contains
1163      ! Introduce the aerosol direct and first indirect radiative forcings:      ! Introduce the aerosol direct and first indirect radiative forcings:
1164      IF (ok_ade .OR. ok_aie) THEN      IF (ok_ade .OR. ok_aie) THEN
1165         ! Get sulfate aerosol distribution :         ! Get sulfate aerosol distribution :
1166         CALL readsulfate(rdayvrai, firstcal, sulfate)         CALL readsulfate(dayvrai, time, firstcal, sulfate)
1167         CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)         CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)
1168    
1169         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &         CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1170              aerindex)              aerindex)
# Line 1408  contains Line 1174  contains
1174         cg_ae = 0.         cg_ae = 0.
1175      ENDIF      ENDIF
1176    
1177      ! Paramètres optiques des nuages et quelques paramètres pour diagnostics :      ! Param\`etres optiques des nuages et quelques param\`etres pour
1178        ! diagnostics :
1179      if (ok_newmicro) then      if (ok_newmicro) then
1180         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1181              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
# Line 1419  contains Line 1186  contains
1186              bl95_b1, cldtaupi, re, fl)              bl95_b1, cldtaupi, re, fl)
1187      endif      endif
1188    
1189      ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.      IF (MOD(itap - 1, radpas) == 0) THEN
1190      IF (MOD(itaprad, radpas) == 0) THEN         ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1191         DO i = 1, klon         ! Calcul de l'abedo moyen par maille
1192            albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &         albsol = sum(falbe * pctsrf, dim = 2)
1193                 + falbe(i, is_lic) * pctsrf(i, is_lic) &  
                + falbe(i, is_ter) * pctsrf(i, is_ter) &  
                + falbe(i, is_sic) * pctsrf(i, is_sic)  
           albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &  
                + falblw(i, is_lic) * pctsrf(i, is_lic) &  
                + falblw(i, is_ter) * pctsrf(i, is_ter) &  
                + falblw(i, is_sic) * pctsrf(i, is_sic)  
        ENDDO  
1194         ! Rayonnement (compatible Arpege-IFS) :         ! Rayonnement (compatible Arpege-IFS) :
1195         CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &         CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1196              albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1197              heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1198              sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1199              lwup, swdn0, swdn, swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, &              swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1200              cg_ae, topswad, solswad, cldtaupi, topswai, solswai)              solswad, cldtaupi, topswai, solswai)
        itaprad = 0  
1201      ENDIF      ENDIF
     itaprad = itaprad + 1  
1202    
1203      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
1204    
1205      DO k = 1, llm      DO k = 1, llm
1206         DO i = 1, klon         DO i = 1, klon
1207            t_seri(i, k) = t_seri(i, k) + (heat(i, k)-cool(i, k)) * dtphys/86400.            t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys/86400.
1208         ENDDO         ENDDO
1209      ENDDO      ENDDO
1210    
1211      IF (if_ebil >= 2) THEN      IF (if_ebil >= 2) THEN
1212         tit = 'after rad'         tit = 'after rad'
1213         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1214              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1215         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1216              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &              zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
1217      END IF      END IF
1218    
1219      ! Calculer l'hydrologie de la surface      ! Calculer l'hydrologie de la surface
# Line 1472  contains Line 1228  contains
1228         ENDDO         ENDDO
1229      ENDDO      ENDDO
1230    
1231      ! Calculer le bilan du sol et la dérive de température (couplage)      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1232    
1233      DO i = 1, klon      DO i = 1, klon
1234         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1235      ENDDO      ENDDO
1236    
1237      ! Paramétrisation de l'orographie à l'échelle sous-maille :      ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1238    
1239      IF (ok_orodr) THEN      IF (ok_orodr) THEN
1240         ! selection des points pour lesquels le shema est actif:         ! S\'election des points pour lesquels le sch\'ema est actif :
1241         igwd = 0         igwd = 0
1242         DO i = 1, klon         DO i = 1, klon
1243            itest(i) = 0            itest(i) = 0
1244            IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1245               itest(i) = 1               itest(i) = 1
1246               igwd = igwd + 1               igwd = igwd + 1
1247               idx(igwd) = i               idx(igwd) = i
# Line 1493  contains Line 1249  contains
1249         ENDDO         ENDDO
1250    
1251         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1252              zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1253              zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1254    
1255         ! ajout des tendances         ! ajout des tendances
1256         DO k = 1, llm         DO k = 1, llm
# Line 1507  contains Line 1263  contains
1263      ENDIF      ENDIF
1264    
1265      IF (ok_orolf) THEN      IF (ok_orolf) THEN
1266         ! Sélection des points pour lesquels le schéma est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
1267         igwd = 0         igwd = 0
1268         DO i = 1, klon         DO i = 1, klon
1269            itest(i) = 0            itest(i) = 0
1270            IF ((zpic(i) - zmea(i)) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
1271               itest(i) = 1               itest(i) = 1
1272               igwd = igwd + 1               igwd = igwd + 1
1273               idx(igwd) = i               idx(igwd) = i
# Line 1532  contains Line 1288  contains
1288         ENDDO         ENDDO
1289      ENDIF      ENDIF
1290    
1291      ! Stress nécessaires : toute la physique      ! Stress n\'ecessaires : toute la physique
1292    
1293      DO i = 1, klon      DO i = 1, klon
1294         zustrph(i) = 0.         zustrph(i) = 0.
# Line 1547  contains Line 1303  contains
1303         ENDDO         ENDDO
1304      ENDDO      ENDDO
1305    
1306      CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &      CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1307           zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)           zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1308    
1309      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &      IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1310           2, dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &           2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1311           d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)           d_qt, d_ec)
1312    
1313      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
1314      call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &      call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, t, &
1315           dtphys, u, t, paprs, play, mfu, mfd, pen_u, pde_u, pen_d, pde_d, &           paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &
1316           ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &           yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, &
1317           frac_nucl, pphis, albsol, rhcl, cldfra, rneb, diafra, cldliq, &           dnwd, tr_seri, zmasse, ncid_startphy, nid_ins)
1318           pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)  
1319        IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1320      IF (offline) THEN           pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1321         call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, pde_u, &           pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
             pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &  
             pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)  
     ENDIF  
1322    
1323      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
1324      CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &      CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
          ue, uq)  
1325    
1326      ! diag. bilKP      ! diag. bilKP
1327    
# Line 1592  contains Line 1344  contains
1344      IF (if_ebil >= 1) THEN      IF (if_ebil >= 1) THEN
1345         tit = 'after physic'         tit = 'after physic'
1346         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &         CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1347              ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &              ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
             d_ql, d_qs, d_ec)  
1348         ! Comme les tendances de la physique sont ajoute dans la dynamique,         ! Comme les tendances de la physique sont ajoute dans la dynamique,
1349         ! on devrait avoir que la variation d'entalpie par la dynamique         ! on devrait avoir que la variation d'entalpie par la dynamique
1350         ! est egale a la variation de la physique au pas de temps precedent.         ! est egale a la variation de la physique au pas de temps precedent.
1351         ! Donc la somme de ces 2 variations devrait etre nulle.         ! Donc la somme de ces 2 variations devrait etre nulle.
1352         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &         call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1353              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &              evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
             fs_bound, fq_bound)  
   
1354         d_h_vcol_phy = d_h_vcol         d_h_vcol_phy = d_h_vcol
   
1355      END IF      END IF
1356    
1357      ! SORTIES      ! SORTIES
# Line 1628  contains Line 1376  contains
1376         ENDDO         ENDDO
1377      ENDDO      ENDDO
1378    
1379      IF (nqmx >= 3) THEN      DO iq = 3, nqmx
1380         DO iq = 3, nqmx         DO k = 1, llm
1381            DO k = 1, llm            DO i = 1, klon
1382               DO i = 1, klon               d_qx(i, k, iq) = (tr_seri(i, k, iq - 2) - qx(i, k, iq)) / dtphys
                 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys  
              ENDDO  
1383            ENDDO            ENDDO
1384         ENDDO         ENDDO
1385      ENDIF      ENDDO
1386    
1387      ! Sauvegarder les valeurs de t et q a la fin de la physique:      ! Sauvegarder les valeurs de t et q a la fin de la physique:
1388      DO k = 1, llm      DO k = 1, llm
# Line 1646  contains Line 1392  contains
1392         ENDDO         ENDDO
1393      ENDDO      ENDDO
1394    
     ! Ecriture des sorties  
     call write_histhf  
     call write_histday  
1395      call write_histins      call write_histins
1396    
1397      ! Si c'est la fin, il faut conserver l'etat de redemarrage      IF (lafin) then
1398      IF (lafin) THEN         call NF95_CLOSE(ncid_startphy)
1399         itau_phy = itau_phy + itap         CALL phyredem(pctsrf, ftsol, ftsoil, tslab, seaice, fqsurf, qsol, &
1400         CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &              fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1401              tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &              radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1402              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &              t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1403              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              w01)
1404              q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)      end IF
     ENDIF  
1405    
1406      firstcal = .FALSE.      firstcal = .FALSE.
1407    
1408    contains    contains
1409    
     subroutine write_histday  
   
       use gr_phy_write_3d_m, only: gr_phy_write_3d  
       integer itau_w ! pas de temps ecriture  
   
       !------------------------------------------------  
   
       if (ok_journe) THEN  
          itau_w = itau_phy + itap  
          if (nqmx <= 4) then  
             call histwrite(nid_day, "Sigma_O3_Royer", itau_w, &  
                  gr_phy_write_3d(wo) * 1e3)  
             ! (convert "wo" from kDU to DU)  
          end if  
          if (ok_sync) then  
             call histsync(nid_day)  
          endif  
       ENDIF  
   
     End subroutine write_histday  
   
     !****************************  
   
     subroutine write_histhf  
   
       ! From phylmd/write_histhf.h, version 1.5 2005/05/25 13:10:09  
   
       !------------------------------------------------  
   
       call write_histhf3d  
   
       IF (ok_sync) THEN  
          call histsync(nid_hf)  
       ENDIF  
   
     end subroutine write_histhf  
   
     !***************************************************************  
   
1410      subroutine write_histins      subroutine write_histins
1411    
1412        ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09        ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09
1413    
1414        real zout        ! Ecriture des sorties
1415        integer itau_w ! pas de temps ecriture  
1416          use dimens_m, only: iim, jjm
1417          USE histsync_m, ONLY: histsync
1418          USE histwrite_m, ONLY: histwrite
1419    
1420          integer i, itau_w ! pas de temps ecriture
1421          REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1422    
1423        !--------------------------------------------------        !--------------------------------------------------
1424    
1425        IF (ok_instan) THEN        IF (ok_instan) THEN
1426           ! Champs 2D:           ! Champs 2D:
1427    
          zsto = dtphys * ecrit_ins  
          zout = dtphys * ecrit_ins  
1428           itau_w = itau_phy + itap           itau_w = itau_phy + itap
1429    
          i = NINT(zout/zsto)  
1430           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)
1431           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1432    
          i = NINT(zout/zsto)  
1433           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)
1434           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1435    
# Line 1794  contains Line 1499  contains
1499           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)
1500           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
1501    
1502           zx_tmp_fi2d(1:klon) = -1*sens(1:klon)           zx_tmp_fi2d(1:klon) = - sens(1:klon)
1503           ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)           ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)
1504           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1505           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
# Line 1856  contains Line 1561  contains
1561              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
1562                   zx_tmp_2d)                   zx_tmp_2d)
1563    
1564              zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)              zx_tmp_fi2d(1 : klon) = falbe(:, nsrf)
1565              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)              CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1566              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &              CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
1567                   zx_tmp_2d)                   zx_tmp_2d)
# Line 1864  contains Line 1569  contains
1569           END DO           END DO
1570           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)
1571           CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)
          CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)  
          CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)  
1572    
1573           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)           CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)
1574           CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)           CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
# Line 1925  contains Line 1628  contains
1628           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)           CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)
1629           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)           CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1630    
1631           if (ok_sync) then           call histsync(nid_ins)
             call histsync(nid_ins)  
          endif  
1632        ENDIF        ENDIF
1633    
1634      end subroutine write_histins      end subroutine write_histins
1635    
     !****************************************************  
   
     subroutine write_histhf3d  
   
       ! From phylmd/write_histhf3d.h, version 1.2 2005/05/25 13:10:09  
   
       integer itau_w ! pas de temps ecriture  
   
       !-------------------------------------------------------  
   
       itau_w = itau_phy + itap  
   
       ! Champs 3D:  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, qx(1, 1, ivap), zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d)  
   
       CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)  
       CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d)  
   
       if (nbtr >= 3) then  
          CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, tr_seri(1, 1, 3), &  
               zx_tmp_3d)  
          CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d)  
       end if  
   
       if (ok_sync) then  
          call histsync(nid_hf3d)  
       endif  
   
     end subroutine write_histhf3d  
   
1636    END SUBROUTINE physiq    END SUBROUTINE physiq
1637    
1638  end module physiq_m  end module physiq_m
Legend:
Removed from v.73  
changed lines
  Added in v.174
  ViewVC Help
Powered by ViewVC 1.1.21