/[lmdze]/trunk/Sources/phylmd/clmain.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/clmain.f

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

trunk/libf/phylmd/clmain.f90 revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC trunk/Sources/phylmd/clmain.f revision 222 by guez, Tue Apr 25 15:31:48 2017 UTC
# Line 4  module clmain_m Line 4  module clmain_m
4    
5  contains  contains
6    
7    SUBROUTINE clmain(dtime, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&    SUBROUTINE clmain(dtime, pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8         jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&         cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, pplay, fsnow, &
9         soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&         qsurf, evap, falbe, fluxlat, rain_fall, snow_f, fsolsw, fsollw, fder, &
10         qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&         frugs, agesno, rugoro, d_t, d_q, d_u, d_v, d_ts, flux_t, flux_q, &
11         rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&         flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, ycoefh, zu1, &
12         cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&         zv1, t2m, q2m, u10m, v10m, pblh, capcl, oliqcl, cteicl, pblt, therm, &
13         d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&         trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
14         dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&  
15         capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
16         fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)      ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
17        ! Objet : interface de couche limite (diffusion verticale)
18      ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19  
19      ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
20      ! Objet : interface de "couche limite" (diffusion verticale)      ! de la couche limite pour les traceurs se fait avec "cltrac" et
21        ! ne tient pas compte de la diff\'erentiation des sous-fractions
22      ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.      ! de sol.
23      ! Pour l'instant le calcul de la couche limite pour les traceurs  
24      ! se fait avec "cltrac" et ne tient pas compte de la différentiation      ! Pour pouvoir extraire les coefficients d'\'echanges et le vent
25      ! des sous-fractions de sol.      ! dans la premi\`ere couche, trois champs ont \'et\'e cr\'e\'es : "ycoefh",
26        ! "zu1" et "zv1". Nous avons moyenn\'e les valeurs de ces trois
27      ! Pour pouvoir extraire les coefficients d'échanges et le vent      ! champs sur les quatre sous-surfaces du mod\`ele.
     ! dans la première couche, trois champs supplémentaires ont été  
     ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons  
     ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces  
     ! du modèle. Dans l'avenir, si les informations des sous-surfaces  
     ! doivent être prises en compte, il faudra sortir ces mêmes champs  
     ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de  
     ! sous-surfaces).  
28    
     use calendar, ONLY : ymds2ju  
29      use clqh_m, only: clqh      use clqh_m, only: clqh
30        use clvent_m, only: clvent
31      use coefkz_m, only: coefkz      use coefkz_m, only: coefkz
32      use coefkzmin_m, only: coefkzmin      use coefkzmin_m, only: coefkzmin
33      USE conf_phys_m, ONLY : iflag_pbl      USE conf_gcm_m, ONLY: prt_level, lmt_pas
34      USE dimens_m, ONLY : iim, jjm      USE conf_phys_m, ONLY: iflag_pbl
35      USE dimphy, ONLY : klev, klon, zmasq      USE dimphy, ONLY: klev, klon, zmasq
36      USE dimsoil, ONLY : nsoilmx      USE dimsoil, ONLY: nsoilmx
     USE dynetat0_m, ONLY : day_ini  
     USE gath_cpl, ONLY : gath2cpl  
37      use hbtm_m, only: hbtm      use hbtm_m, only: hbtm
38      USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
39      use histwrite_m, only: histwrite      USE interfoce_lim_m, ONLY: interfoce_lim
40      USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf      use stdlevvar_m, only: stdlevvar
41      USE conf_gcm_m, ONLY : prt_level      USE suphec_m, ONLY: rd, rg, rkappa
42      USE suphec_m, ONLY : rd, rg, rkappa      use time_phylmdz, only: itap
43      USE temps, ONLY : annee_ref, itau_phy      use ustarhb_m, only: ustarhb
44        use vdif_kcay_m, only: vdif_kcay
45      use yamada4_m, only: yamada4      use yamada4_m, only: yamada4
46    
47      ! Arguments:      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
48    
49        REAL, INTENT(inout):: pctsrf(klon, nbsrf)
50        ! tableau des pourcentages de surface de chaque maille
51    
52        REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
53        REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
54        REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
55        INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
56        REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
57        REAL, INTENT(IN):: ftsol(:, :) ! (klon, nbsrf) temp\'erature du sol (en K)
58        REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
59        REAL, INTENT(IN):: ksta, ksta_ter
60        LOGICAL, INTENT(IN):: ok_kzmin
61    
62        REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
63        ! soil temperature of surface fraction
64    
65        REAL, INTENT(inout):: qsol(klon)
66        ! column-density of water in soil, in kg m-2
67    
68        REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
69        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
70        REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf) \'epaisseur neigeuse
71        REAL qsurf(klon, nbsrf)
72        REAL evap(klon, nbsrf)
73        REAL, intent(inout):: falbe(klon, nbsrf)
74        REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
75    
76        REAL, intent(in):: rain_fall(klon)
77        ! liquid water mass flux (kg/m2/s), positive down
78    
79        REAL, intent(in):: snow_f(klon)
80        ! solid water mass flux (kg/m2/s), positive down
81    
82        REAL, INTENT(IN):: fsolsw(klon, nbsrf), fsollw(klon, nbsrf)
83        REAL, intent(in):: fder(:) ! (klon)
84        REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
85        real agesno(klon, nbsrf)
86        REAL, INTENT(IN):: rugoro(klon)
87    
     REAL, INTENT (IN) :: dtime ! interval du temps (secondes)  
     REAL date0  
     ! date0----input-R- jour initial  
     INTEGER, INTENT (IN) :: itap  
     ! itap-----input-I- numero du pas de temps  
     REAL, INTENT(IN):: t(klon, klev), q(klon, klev)  
     ! t--------input-R- temperature (K)  
     ! q--------input-R- vapeur d'eau (kg/kg)  
     REAL, INTENT (IN):: u(klon, klev), v(klon, klev)  
     ! u--------input-R- vitesse u  
     ! v--------input-R- vitesse v  
     REAL, INTENT (IN):: paprs(klon, klev+1)  
     ! paprs----input-R- pression a intercouche (Pa)  
     REAL, INTENT (IN):: pplay(klon, klev)  
     ! pplay----input-R- pression au milieu de couche (Pa)  
     REAL, INTENT (IN):: rlon(klon), rlat(klon)  
     ! rlat-----input-R- latitude en degree  
     REAL cufi(klon), cvfi(klon)  
     ! cufi-----input-R- resolution des mailles en x (m)  
     ! cvfi-----input-R- resolution des mailles en y (m)  
88      REAL d_t(klon, klev), d_q(klon, klev)      REAL d_t(klon, klev), d_q(klon, klev)
89      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
90      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
91      REAL d_u(klon, klev), d_v(klon, klev)  
92      ! d_u------output-R- le changement pour "u"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
93      ! d_v------output-R- le changement pour "v"      ! changement pour "u" et "v"
94      REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)  
95      ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)      REAL, intent(out):: d_ts(:, :) ! (klon, nbsrf) variation of ftsol
96      !                    (orientation positive vers le bas)  
97      ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)      REAL, intent(out):: flux_t(klon, nbsrf)
98      REAL dflux_t(klon), dflux_q(klon)      ! flux de chaleur sensible (Cp T) (W/m2) (orientation positive vers
99        ! le bas) à la surface
100    
101        REAL, intent(out):: flux_q(klon, nbsrf)
102        ! flux de vapeur d'eau (kg/m2/s) à la surface
103    
104        REAL, intent(out):: flux_u(klon, nbsrf), flux_v(klon, nbsrf)
105        ! tension du vent à la surface, en Pa
106    
107        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
108        real q2(klon, klev+1, nbsrf)
109    
110        REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
111      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
112      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
113      !IM "slab" ocean      ! IM "slab" ocean
114      REAL flux_o(klon), flux_g(klon)  
115      !IM "slab" ocean      REAL, intent(out):: ycoefh(klon, klev)
116      ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')      REAL, intent(out):: zu1(klon)
117      ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')      REAL zv1(klon)
118      REAL y_flux_o(klon), y_flux_g(klon)      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
119      REAL tslab(klon), ytslab(klon)      REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
120      ! tslab-in/output-R temperature du slab ocean (en Kelvin)  
121      ! uniqmnt pour slab      ! Ionela Musat cf. Anne Mathieu : planetary boundary layer, hbtm
122      REAL seaice(klon), y_seaice(klon)      ! (Comme les autres diagnostics on cumule dans physiq ce qui
123      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      ! permet de sortir les grandeurs par sous-surface)
124      REAL y_fqcalving(klon), y_ffonte(klon)      REAL pblh(klon, nbsrf) ! height of planetary boundary layer
125        REAL capcl(klon, nbsrf)
126        REAL oliqcl(klon, nbsrf)
127        REAL cteicl(klon, nbsrf)
128        REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
129        REAL therm(klon, nbsrf)
130        REAL trmb1(klon, nbsrf)
131        ! trmb1-------deep_cape
132        REAL trmb2(klon, nbsrf)
133        ! trmb2--------inhibition
134        REAL trmb3(klon, nbsrf)
135        ! trmb3-------Point Omega
136        REAL plcl(klon, nbsrf)
137      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
138      ! ffonte----Flux thermique utilise pour fondre la neige      ! ffonte----Flux thermique utilise pour fondre la neige
139      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
140      !           hauteur de neige, en kg/m2/s      !           hauteur de neige, en kg/m2/s
141      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      REAL run_off_lic_0(klon)
142    
143      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      ! Local:
     ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal  
     REAL rugmer(klon), agesno(klon, nbsrf)  
     REAL, INTENT (IN) :: rugoro(klon)  
     REAL cdragh(klon), cdragm(klon)  
     ! jour de l'annee en cours                  
     INTEGER jour  
     REAL rmu0(klon) ! cosinus de l'angle solaire zenithal      
     ! taux CO2 atmosphere                      
     REAL co2_ppm  
     LOGICAL, INTENT (IN) :: debut  
     LOGICAL, INTENT (IN) :: lafin  
     LOGICAL ok_veget  
     CHARACTER (len=*), INTENT (IN) :: ocean  
     INTEGER npas, nexca  
   
     REAL pctsrf(klon, nbsrf)  
     REAL ts(klon, nbsrf)  
     ! ts-------input-R- temperature du sol (en Kelvin)  
     REAL d_ts(klon, nbsrf)  
     ! d_ts-----output-R- le changement pour "ts"  
     REAL snow(klon, nbsrf)  
     REAL qsurf(klon, nbsrf)  
     REAL evap(klon, nbsrf)  
     REAL albe(klon, nbsrf)  
     REAL alblw(klon, nbsrf)  
144    
145      REAL fluxlat(klon, nbsrf)      LOGICAL:: firstcal = .true.
146    
     REAL rain_f(klon), snow_f(klon)  
     REAL fder(klon)  
   
     REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)  
     REAL rugos(klon, nbsrf)  
     ! rugos----input-R- longeur de rugosite (en m)  
147      ! la nouvelle repartition des surfaces sortie de l'interface      ! la nouvelle repartition des surfaces sortie de l'interface
148      REAL pctsrf_new(klon, nbsrf)      REAL, save:: pctsrf_new_oce(klon)
149        REAL, save:: pctsrf_new_sic(klon)
150    
151      REAL zcoefh(klon, klev)      REAL y_fqcalving(klon), y_ffonte(klon)
152      REAL zu1(klon)      real y_run_off_lic_0(klon)
153      REAL zv1(klon)      REAL rugmer(klon)
154        REAL ytsoil(klon, nsoilmx)
155        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
156        REAL yalb(klon)
157    
158      !$$$ PB ajout pour soil      REAL yu1(klon), yv1(klon)
159      LOGICAL, INTENT (IN) :: soil_model      ! On ajoute en output yu1 et yv1 qui sont les vents dans
160      !IM ajout seuils cdrm, cdrh      ! la premi\`ere couche.
161      REAL cdmmax, cdhmax      
162        REAL snow(klon), yqsurf(klon), yagesno(klon)
163    
164      REAL ksta, ksta_ter      real yqsol(klon)
165      LOGICAL ok_kzmin      ! column-density of water in soil, in kg m-2
166    
167      REAL ftsoil(klon, nsoilmx, nbsrf)      REAL yrain_f(klon)
168      REAL ytsoil(klon, nsoilmx)      ! liquid water mass flux (kg/m2/s), positive down
     REAL qsol(klon)  
169    
170      EXTERNAL clvent, calbeta, cltrac      REAL ysnow_f(klon)
171        ! solid water mass flux (kg/m2/s), positive down
172    
173      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yfder(klon)
     REAL yalb(klon)  
     REAL yalblw(klon)  
     REAL yu1(klon), yv1(klon)  
     ! on rajoute en output yu1 et yv1 qui sont les vents dans  
     ! la premiere couche  
     REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)  
     REAL yrain_f(klon), ysnow_f(klon)  
     REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)  
     REAL yfder(klon), ytaux(klon), ytauy(klon)  
174      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
   
175      REAL yfluxlat(klon)      REAL yfluxlat(klon)
   
176      REAL y_d_ts(klon)      REAL y_d_ts(klon)
177      REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
178      REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
179      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon), y_flux_q(klon)
180      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon), y_flux_v(klon)
181      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
182      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
183      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
184      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
185      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
186    
     LOGICAL ok_nonloc  
     PARAMETER (ok_nonloc=.FALSE.)  
187      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
188    
189      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
190      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
191      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
192      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1)
193      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
194    
195      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
# Line 204  contains Line 199  contains
199      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
200    
201      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
202      ! "pourcentage potentiel" pour tenir compte des éventuelles      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
203      ! apparitions ou disparitions de la glace de mer      ! apparitions ou disparitions de la glace de mer
204    
205      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL zx_alf1, zx_alf2 ! valeur ambiante par extrapolation
   
     ! maf pour sorties IOISPL en cas de debugagage  
   
     CHARACTER (80) cldebug  
     SAVE cldebug  
     CHARACTER (8) cl_surf(nbsrf)  
     SAVE cl_surf  
     INTEGER nhoridbg, nidbg  
     SAVE nhoridbg, nidbg  
     INTEGER ndexbg(iim*(jjm+1))  
     REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian  
     REAL tabindx(klon)  
     REAL debugtab(iim, jjm+1)  
     LOGICAL first_appel  
     SAVE first_appel  
     DATA first_appel/ .TRUE./  
     LOGICAL :: debugindex = .FALSE.  
     INTEGER idayref  
     REAL t2m(klon, nbsrf), q2m(klon, nbsrf)  
     REAL u10m(klon, nbsrf), v10m(klon, nbsrf)  
206    
207      REAL yt2m(klon), yq2m(klon), yu10m(klon)      REAL yt2m(klon), yq2m(klon), yu10m(klon)
208      REAL yustar(klon)      REAL yustar(klon)
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
209    
210      REAL yt10m(klon), yq10m(klon)      REAL yt10m(klon), yq10m(klon)
     !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds  
     ! physiq ce qui permet de sortir les grdeurs par sous surface)  
     REAL pblh(klon, nbsrf)  
     ! pblh------- HCL  
     REAL plcl(klon, nbsrf)  
     REAL capcl(klon, nbsrf)  
     REAL oliqcl(klon, nbsrf)  
     REAL cteicl(klon, nbsrf)  
     REAL pblt(klon, nbsrf)  
     ! pblT------- T au nveau HCL  
     REAL therm(klon, nbsrf)  
     REAL trmb1(klon, nbsrf)  
     ! trmb1-------deep_cape  
     REAL trmb2(klon, nbsrf)  
     ! trmb2--------inhibition  
     REAL trmb3(klon, nbsrf)  
     ! trmb3-------Point Omega  
211      REAL ypblh(klon)      REAL ypblh(klon)
212      REAL ylcl(klon)      REAL ylcl(klon)
213      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 265  contains Line 218  contains
218      REAL ytrmb1(klon)      REAL ytrmb1(klon)
219      REAL ytrmb2(klon)      REAL ytrmb2(klon)
220      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
221      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
222      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
223      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 277  contains Line 229  contains
229      LOGICAL zxli      LOGICAL zxli
230      PARAMETER (zxli=.FALSE.)      PARAMETER (zxli=.FALSE.)
231    
     REAL zt, zqs, zdelta, zcor  
     REAL t_coup  
     PARAMETER (t_coup=273.15)  
   
     CHARACTER (len=20) :: modname = 'clmain'  
   
232      !------------------------------------------------------------      !------------------------------------------------------------
233    
234      ytherm = 0.      ytherm = 0.
235    
     IF (debugindex .AND. first_appel) THEN  
        first_appel = .FALSE.  
   
        ! initialisation sorties netcdf  
   
        idayref = day_ini  
        CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian)  
        CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)  
        DO i = 1, iim  
           zx_lon(i, 1) = rlon(i+1)  
           zx_lon(i, jjm+1) = rlon(i+1)  
        END DO  
        CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)  
        cldebug = 'sous_index'  
        CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &  
             iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)  
        ! no vertical axis  
        cl_surf(1) = 'ter'  
        cl_surf(2) = 'lic'  
        cl_surf(3) = 'oce'  
        cl_surf(4) = 'sic'  
        DO nsrf = 1, nbsrf  
           CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &  
                nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)  
        END DO  
        CALL histend(nidbg)  
        CALL histsync(nidbg)  
     END IF  
   
236      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
237         DO i = 1, klon         DO i = 1, klon
238            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k+1)
# Line 337  contains Line 254  contains
254      zu1 = 0.      zu1 = 0.
255      zv1 = 0.      zv1 = 0.
256      ypct = 0.      ypct = 0.
     yts = 0.  
     ysnow = 0.  
257      yqsurf = 0.      yqsurf = 0.
     yalb = 0.  
     yalblw = 0.  
258      yrain_f = 0.      yrain_f = 0.
259      ysnow_f = 0.      ysnow_f = 0.
     yfder = 0.  
     ytaux = 0.  
     ytauy = 0.  
     ysolsw = 0.  
     ysollw = 0.  
     ysollwdown = 0.  
260      yrugos = 0.      yrugos = 0.
261      yu1 = 0.      yu1 = 0.
262      yv1 = 0.      yv1 = 0.
     yrads = 0.  
263      ypaprs = 0.      ypaprs = 0.
264      ypplay = 0.      ypplay = 0.
265      ydelp = 0.      ydelp = 0.
# Line 361  contains Line 267  contains
267      yv = 0.      yv = 0.
268      yt = 0.      yt = 0.
269      yq = 0.      yq = 0.
     pctsrf_new = 0.  
     y_flux_u = 0.  
     y_flux_v = 0.  
     !$$ PB  
270      y_dflux_t = 0.      y_dflux_t = 0.
271      y_dflux_q = 0.      y_dflux_q = 0.
     ytsoil = 999999.  
272      yrugoro = 0.      yrugoro = 0.
     ! -- LOOP  
     yu10mx = 0.  
     yu10my = 0.  
     ywindsp = 0.  
     ! -- LOOP  
273      d_ts = 0.      d_ts = 0.
     !§§§ PB  
     yfluxlat = 0.  
274      flux_t = 0.      flux_t = 0.
275      flux_q = 0.      flux_q = 0.
276      flux_u = 0.      flux_u = 0.
277      flux_v = 0.      flux_v = 0.
278        fluxlat = 0.
279      d_t = 0.      d_t = 0.
280      d_q = 0.      d_q = 0.
281      d_u = 0.      d_u = 0.
282      d_v = 0.      d_v = 0.
283      zcoefh = 0.      ycoefh = 0.
284    
285      ! Boucler sur toutes les sous-fractions du sol:      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
286        ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
287      ! Initialisation des "pourcentages potentiels". On considère ici qu'on      ! (\`a affiner)
     ! peut avoir potentiellement de la glace sur tout le domaine océanique  
     ! (à affiner)  
288    
289      pctsrf_pot = pctsrf      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
290        pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
291      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
292      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
293    
294        ! Tester si c'est le moment de lire le fichier:
295        if (mod(itap - 1, lmt_pas) == 0) then
296           CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
297        endif
298    
299        ! Boucler sur toutes les sous-fractions du sol:
300    
301      loop_surface: DO nsrf = 1, nbsrf      loop_surface: DO nsrf = 1, nbsrf
302         ! Chercher les indices :         ! Chercher les indices :
303         ni = 0         ni = 0
304         knon = 0         knon = 0
305         DO i = 1, klon         DO i = 1, klon
306            ! Pour déterminer le domaine à traiter, on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
307            ! "potentielles"            ! "potentielles"
308            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
309               knon = knon + 1               knon = knon + 1
# Line 410  contains Line 311  contains
311            END IF            END IF
312         END DO         END DO
313    
314         ! variables pour avoir une sortie IOIPSL des INDEX         if_knon: IF (knon /= 0) then
        IF (debugindex) THEN  
           tabindx = 0.  
           DO i = 1, knon  
              tabindx(i) = real(i)  
           END DO  
           debugtab = 0.  
           ndexbg = 0  
           CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)  
           CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)  
        END IF  
   
        IF (knon == 0) CYCLE  
   
        DO j = 1, knon  
           i = ni(j)  
           ypct(j) = pctsrf(i, nsrf)  
           yts(j) = ts(i, nsrf)  
           ytslab(i) = tslab(i)  
           ysnow(j) = snow(i, nsrf)  
           yqsurf(j) = qsurf(i, nsrf)  
           yalb(j) = albe(i, nsrf)  
           yalblw(j) = alblw(i, nsrf)  
           yrain_f(j) = rain_f(i)  
           ysnow_f(j) = snow_f(i)  
           yagesno(j) = agesno(i, nsrf)  
           yfder(j) = fder(i)  
           ytaux(j) = flux_u(i, 1, nsrf)  
           ytauy(j) = flux_v(i, 1, nsrf)  
           ysolsw(j) = solsw(i, nsrf)  
           ysollw(j) = sollw(i, nsrf)  
           ysollwdown(j) = sollwdown(i)  
           yrugos(j) = rugos(i, nsrf)  
           yrugoro(j) = rugoro(i)  
           yu1(j) = u1lay(i)  
           yv1(j) = v1lay(i)  
           yrads(j) = ysolsw(j) + ysollw(j)  
           ypaprs(j, klev+1) = paprs(i, klev+1)  
           y_run_off_lic_0(j) = run_off_lic_0(i)  
           yu10mx(j) = u10m(i, nsrf)  
           yu10my(j) = v10m(i, nsrf)  
           ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))  
        END DO  
   
        ! IF bucket model for continent, copy soil water content  
        IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN  
315            DO j = 1, knon            DO j = 1, knon
316               i = ni(j)               i = ni(j)
317               yqsol(j) = qsol(i)               ypct(j) = pctsrf(i, nsrf)
318            END DO               yts(j) = ftsol(i, nsrf)
319         ELSE               snow(j) = fsnow(i, nsrf)
320            yqsol = 0.               yqsurf(j) = qsurf(i, nsrf)
321         END IF               yalb(j) = falbe(i, nsrf)
322         !$$$ PB ajour pour soil               yrain_f(j) = rain_fall(i)
323         DO k = 1, nsoilmx               ysnow_f(j) = snow_f(i)
324            DO j = 1, knon               yagesno(j) = agesno(i, nsrf)
325               i = ni(j)               yfder(j) = fder(i)
326               ytsoil(j, k) = ftsoil(i, k, nsrf)               yrugos(j) = frugs(i, nsrf)
327            END DO               yrugoro(j) = rugoro(i)
328         END DO               yu1(j) = u1lay(i)
329         DO k = 1, klev               yv1(j) = v1lay(i)
330            DO j = 1, knon               yrads(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
331               i = ni(j)               ypaprs(j, klev+1) = paprs(i, klev+1)
332               ypaprs(j, k) = paprs(i, k)               y_run_off_lic_0(j) = run_off_lic_0(i)
333               ypplay(j, k) = pplay(i, k)            END DO
334               ydelp(j, k) = delp(i, k)  
335               yu(j, k) = u(i, k)            ! For continent, copy soil water content
336               yv(j, k) = v(i, k)            IF (nsrf == is_ter) THEN
337               yt(j, k) = t(i, k)               yqsol(:knon) = qsol(ni(:knon))
338               yq(j, k) = q(i, k)            ELSE
339            END DO               yqsol = 0.
340         END DO            END IF
341    
342              ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
343    
        ! calculer Cdrag et les coefficients d'echange  
        CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&  
             yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)  
        IF (iflag_pbl == 1) THEN  
           CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)  
344            DO k = 1, klev            DO k = 1, klev
345               DO i = 1, knon               DO j = 1, knon
346                  ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))                  i = ni(j)
347                  ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))                  ypaprs(j, k) = paprs(i, k)
348                    ypplay(j, k) = pplay(i, k)
349                    ydelp(j, k) = delp(i, k)
350                    yu(j, k) = u(i, k)
351                    yv(j, k) = v(i, k)
352                    yt(j, k) = t(i, k)
353                    yq(j, k) = q(i, k)
354               END DO               END DO
355            END DO            END DO
        END IF  
356    
357         ! on seuille ycoefm et ycoefh            ! calculer Cdrag et les coefficients d'echange
358         IF (nsrf == is_oce) THEN            CALL coefkz(nsrf, ypaprs, ypplay, ksta, ksta_ter, yts(:knon), &
359            DO j = 1, knon                 yrugos, yu, yv, yt, yq, yqsurf(:knon), coefm(:knon, :), &
360               ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)                 coefh(:knon, :))
361               ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)            IF (iflag_pbl == 1) THEN
362            END DO               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
363         END IF               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
364                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
365              END IF
366    
367         IF (ok_kzmin) THEN            ! on met un seuil pour coefm et coefh
368            ! Calcul d'une diffusion minimale pour les conditions tres stables            IF (nsrf == is_oce) THEN
369            CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm(:, 1), &               coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
370                 ycoefm0, ycoefh0)               coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
371              END IF
372    
373            DO k = 1, klev            IF (ok_kzmin) THEN
374               DO i = 1, knon               ! Calcul d'une diffusion minimale pour les conditions tres stables
375                  ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))               CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
376                  ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))                    coefm(:knon, 1), ycoefm0, ycoefh0)
377               END DO               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
378            END DO               coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
379         END IF            END IF
380    
381         IF (iflag_pbl >= 3) THEN            IF (iflag_pbl >= 3) THEN
382            ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
383            yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &               ! Fr\'ed\'eric Hourdin
384                 1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
385            DO k = 2, klev                    + ypplay(:knon, 1))) &
386               yzlay(1:knon, k) = yzlay(1:knon, k-1) &                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
387                    + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &               DO k = 2, klev
388                    / ypaprs(1:knon, k) &                  yzlay(1:knon, k) = yzlay(1:knon, k-1) &
389                    * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg                       + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
390            END DO                       / ypaprs(1:knon, k) &
391            DO k = 1, klev                       * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
392               yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &               END DO
393                    / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))               DO k = 1, klev
394            END DO                  yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
395            yzlev(1:knon, 1) = 0.                       / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
396            yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)               END DO
397            DO k = 2, klev               yzlev(1:knon, 1) = 0.
398               yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))               yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
399            END DO                    - yzlay(:knon, klev - 1)
400            DO k = 1, klev + 1               DO k = 2, klev
401               DO j = 1, knon                  yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
402                  i = ni(j)               END DO
403                  yq2(j, k) = q2(i, k, nsrf)               DO k = 1, klev + 1
404                    DO j = 1, knon
405                       i = ni(j)
406                       yq2(j, k) = q2(i, k, nsrf)
407                    END DO
408               END DO               END DO
           END DO  
409    
410            y_cd_m(1:knon) = ycoefm(1:knon, 1)               CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
411            y_cd_h(1:knon) = ycoefh(1:knon, 1)               IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
412    
413            IF (prt_level>9) THEN               ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
              PRINT *, 'USTAR = ', yustar  
           END IF  
414    
415            ! iflag_pbl peut être utilisé comme longueur de mélange               IF (iflag_pbl >= 11) THEN
416                    CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
417                         yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
418                         iflag_pbl)
419                 ELSE
420                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
421                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
422                 END IF
423    
424            IF (iflag_pbl >= 11) THEN               coefm(:knon, 2:) = ykmm(:knon, 2:klev)
425               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               coefh(:knon, 2:) = ykmn(:knon, 2:klev)
                   yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &  
                   iflag_pbl)  
           ELSE  
              CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &  
                   y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)  
426            END IF            END IF
427    
428            ycoefm(1:knon, 1) = y_cd_m(1:knon)            ! calculer la diffusion des vitesses "u" et "v"
429            ycoefh(1:knon, 1) = y_cd_h(1:knon)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
430            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay, ydelp, y_d_u, y_flux_u(:knon))
431            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
432         END IF                 ypplay, ydelp, y_d_v, y_flux_v(:knon))
433    
434         ! calculer la diffusion des vitesses "u" et "v"            ! calculer la diffusion de "q" et de "h"
435         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &            CALL clqh(dtime, julien, firstcal, nsrf, ni(:knon), &
436              ydelp, y_d_u, y_flux_u)                 ytsoil(:knon, :), yqsol, mu0, yrugos, yrugoro, yu1, yv1, &
437         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                 coefh(:knon, :), yt, yq, yts(:knon), ypaprs, ypplay, ydelp, &
438              ydelp, y_d_v, y_flux_v)                 yrads(:knon), yalb(:knon), snow(:knon), yqsurf, yrain_f, &
439                   ysnow_f, yfder(:knon), yfluxlat(:knon), pctsrf_new_sic, &
440         ! pour le couplage                 yagesno(:knon), y_d_t, y_d_q, y_d_ts(:knon), yz0_new, &
441         ytaux = y_flux_u(:, 1)                 y_flux_t(:knon), y_flux_q(:knon), y_dflux_t(:knon), &
442         ytauy = y_flux_v(:, 1)                 y_dflux_q(:knon), y_fqcalving, y_ffonte, y_run_off_lic_0)
   
        ! calculer la diffusion de "q" et de "h"  
        CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&  
             cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&  
             yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&  
             yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&  
             ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &  
             yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&  
             yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&  
             yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&  
             y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&  
             ytslab, y_seaice)  
   
        ! calculer la longueur de rugosite sur ocean  
        yrugm = 0.  
        IF (nsrf == is_oce) THEN  
           DO j = 1, knon  
              yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &  
                   0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))  
              yrugm(j) = max(1.5E-05, yrugm(j))  
           END DO  
        END IF  
        DO j = 1, knon  
           y_dflux_t(j) = y_dflux_t(j)*ypct(j)  
           y_dflux_q(j) = y_dflux_q(j)*ypct(j)  
           yu1(j) = yu1(j)*ypct(j)  
           yv1(j) = yv1(j)*ypct(j)  
        END DO  
443    
444         DO k = 1, klev            ! calculer la longueur de rugosite sur ocean
445              yrugm = 0.
446              IF (nsrf == is_oce) THEN
447                 DO j = 1, knon
448                    yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
449                         0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
450                    yrugm(j) = max(1.5E-05, yrugm(j))
451                 END DO
452              END IF
453            DO j = 1, knon            DO j = 1, knon
454               i = ni(j)               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
455               ycoefh(j, k) = ycoefh(j, k)*ypct(j)               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
456               ycoefm(j, k) = ycoefm(j, k)*ypct(j)               yu1(j) = yu1(j)*ypct(j)
457               y_d_t(j, k) = y_d_t(j, k)*ypct(j)               yv1(j) = yv1(j)*ypct(j)
              y_d_q(j, k) = y_d_q(j, k)*ypct(j)  
              flux_t(i, k, nsrf) = y_flux_t(j, k)  
              flux_q(i, k, nsrf) = y_flux_q(j, k)  
              flux_u(i, k, nsrf) = y_flux_u(j, k)  
              flux_v(i, k, nsrf) = y_flux_v(j, k)  
              y_d_u(j, k) = y_d_u(j, k)*ypct(j)  
              y_d_v(j, k) = y_d_v(j, k)*ypct(j)  
458            END DO            END DO
        END DO  
459    
460         evap(:, nsrf) = -flux_q(:, 1, nsrf)            DO k = 1, klev
461                 DO j = 1, knon
462                    i = ni(j)
463                    coefh(j, k) = coefh(j, k)*ypct(j)
464                    coefm(j, k) = coefm(j, k)*ypct(j)
465                    y_d_t(j, k) = y_d_t(j, k)*ypct(j)
466                    y_d_q(j, k) = y_d_q(j, k)*ypct(j)
467                    y_d_u(j, k) = y_d_u(j, k)*ypct(j)
468                    y_d_v(j, k) = y_d_v(j, k)*ypct(j)
469                 END DO
470              END DO
471    
472         albe(:, nsrf) = 0.            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
473         alblw(:, nsrf) = 0.            flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
474         snow(:, nsrf) = 0.            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
475         qsurf(:, nsrf) = 0.            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
476         rugos(:, nsrf) = 0.  
477         fluxlat(:, nsrf) = 0.            evap(:, nsrf) = -flux_q(:, nsrf)
478         DO j = 1, knon  
479            i = ni(j)            falbe(:, nsrf) = 0.
480            d_ts(i, nsrf) = y_d_ts(j)            fsnow(:, nsrf) = 0.
481            albe(i, nsrf) = yalb(j)            qsurf(:, nsrf) = 0.
482            alblw(i, nsrf) = yalblw(j)            frugs(:, nsrf) = 0.
           snow(i, nsrf) = ysnow(j)  
           qsurf(i, nsrf) = yqsurf(j)  
           rugos(i, nsrf) = yz0_new(j)  
           fluxlat(i, nsrf) = yfluxlat(j)  
           IF (nsrf == is_oce) THEN  
              rugmer(i) = yrugm(j)  
              rugos(i, nsrf) = yrugm(j)  
           END IF  
           agesno(i, nsrf) = yagesno(j)  
           fqcalving(i, nsrf) = y_fqcalving(j)  
           ffonte(i, nsrf) = y_ffonte(j)  
           cdragh(i) = cdragh(i) + ycoefh(j, 1)  
           cdragm(i) = cdragm(i) + ycoefm(j, 1)  
           dflux_t(i) = dflux_t(i) + y_dflux_t(j)  
           dflux_q(i) = dflux_q(i) + y_dflux_q(j)  
           zu1(i) = zu1(i) + yu1(j)  
           zv1(i) = zv1(i) + yv1(j)  
        END DO  
        IF (nsrf == is_ter) THEN  
483            DO j = 1, knon            DO j = 1, knon
484               i = ni(j)               i = ni(j)
485               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
486            END DO               falbe(i, nsrf) = yalb(j)
487         END IF               fsnow(i, nsrf) = snow(j)
488         IF (nsrf == is_lic) THEN               qsurf(i, nsrf) = yqsurf(j)
489                 frugs(i, nsrf) = yz0_new(j)
490                 fluxlat(i, nsrf) = yfluxlat(j)
491                 IF (nsrf == is_oce) THEN
492                    rugmer(i) = yrugm(j)
493                    frugs(i, nsrf) = yrugm(j)
494                 END IF
495                 agesno(i, nsrf) = yagesno(j)
496                 fqcalving(i, nsrf) = y_fqcalving(j)
497                 ffonte(i, nsrf) = y_ffonte(j)
498                 cdragh(i) = cdragh(i) + coefh(j, 1)
499                 cdragm(i) = cdragm(i) + coefm(j, 1)
500                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
501                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
502                 zu1(i) = zu1(i) + yu1(j)
503                 zv1(i) = zv1(i) + yv1(j)
504              END DO
505              IF (nsrf == is_ter) THEN
506                 qsol(ni(:knon)) = yqsol(:knon)
507              else IF (nsrf == is_lic) THEN
508                 DO j = 1, knon
509                    i = ni(j)
510                    run_off_lic_0(i) = y_run_off_lic_0(j)
511                 END DO
512              END IF
513    
514              ftsoil(:, :, nsrf) = 0.
515              ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
516    
517            DO j = 1, knon            DO j = 1, knon
518               i = ni(j)               i = ni(j)
519               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
520                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
521                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
522                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
523                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
524                    ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
525                 END DO
526            END DO            END DO
527         END IF  
528         !$$$ PB ajout pour soil            ! diagnostic t, q a 2m et u, v a 10m
529         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
530            DO j = 1, knon            DO j = 1, knon
531               i = ni(j)               i = ni(j)
532               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
533            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
534         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
535                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
536                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
537                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
538                 tairsol(j) = yts(j) + y_d_ts(j)
539                 rugo1(j) = yrugos(j)
540                 IF (nsrf == is_oce) THEN
541                    rugo1(j) = frugs(i, nsrf)
542                 END IF
543                 psfce(j) = ypaprs(j, 1)
544                 patm(j) = ypplay(j, 1)
545    
546         DO j = 1, knon               qairsol(j) = yqsurf(j)
           i = ni(j)  
           DO k = 1, klev  
              d_t(i, k) = d_t(i, k) + y_d_t(j, k)  
              d_q(i, k) = d_q(i, k) + y_d_q(j, k)  
              d_u(i, k) = d_u(i, k) + y_d_u(j, k)  
              d_v(i, k) = d_v(i, k) + y_d_v(j, k)  
              zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)  
547            END DO            END DO
        END DO  
   
        !cc diagnostic t, q a 2m et u, v a 10m  
   
        DO j = 1, knon  
           i = ni(j)  
           uzon(j) = yu(j, 1) + y_d_u(j, 1)  
           vmer(j) = yv(j, 1) + y_d_v(j, 1)  
           tair1(j) = yt(j, 1) + y_d_t(j, 1)  
           qair1(j) = yq(j, 1) + y_d_q(j, 1)  
           zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &  
                1)))*(ypaprs(j, 1)-ypplay(j, 1))  
           tairsol(j) = yts(j) + y_d_ts(j)  
           rugo1(j) = yrugos(j)  
           IF (nsrf == is_oce) THEN  
              rugo1(j) = rugos(i, nsrf)  
           END IF  
           psfce(j) = ypaprs(j, 1)  
           patm(j) = ypplay(j, 1)  
548    
549            qairsol(j) = yqsurf(j)            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
550         END DO                 zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
551                   yt10m, yq10m, yu10m, yustar)
552    
553         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &            DO j = 1, knon
554              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               i = ni(j)
555              yu10m, yustar)               t2m(i, nsrf) = yt2m(j)
556                 q2m(i, nsrf) = yq2m(j)
        DO j = 1, knon  
           i = ni(j)  
           t2m(i, nsrf) = yt2m(j)  
           q2m(i, nsrf) = yq2m(j)  
   
           ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman  
           u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)  
           v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)  
557    
558         END DO               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
559                 u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
560                 v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
561              END DO
562    
563         DO i = 1, knon            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t(:knon), &
564            y_cd_h(i) = ycoefh(i, 1)                 y_flux_q(:knon), yu, yv, yt, yq, ypblh(:knon), ycapcl, &
565            y_cd_m(i) = ycoefm(i, 1)                 yoliqcl, ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
        END DO  
        CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &  
             y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &  
             ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)  
   
        DO j = 1, knon  
           i = ni(j)  
           pblh(i, nsrf) = ypblh(j)  
           plcl(i, nsrf) = ylcl(j)  
           capcl(i, nsrf) = ycapcl(j)  
           oliqcl(i, nsrf) = yoliqcl(j)  
           cteicl(i, nsrf) = ycteicl(j)  
           pblt(i, nsrf) = ypblt(j)  
           therm(i, nsrf) = ytherm(j)  
           trmb1(i, nsrf) = ytrmb1(j)  
           trmb2(i, nsrf) = ytrmb2(j)  
           trmb3(i, nsrf) = ytrmb3(j)  
        END DO  
566    
        DO j = 1, knon  
           DO k = 1, klev + 1  
              i = ni(j)  
              q2(i, k, nsrf) = yq2(j, k)  
           END DO  
        END DO  
        !IM "slab" ocean  
        IF (nsrf == is_oce) THEN  
567            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
568               i = ni(j)               i = ni(j)
569               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
570                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
571               ELSE               capcl(i, nsrf) = ycapcl(j)
572                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
573               END IF               cteicl(i, nsrf) = ycteicl(j)
574                 pblt(i, nsrf) = ypblt(j)
575                 therm(i, nsrf) = ytherm(j)
576                 trmb1(i, nsrf) = ytrmb1(j)
577                 trmb2(i, nsrf) = ytrmb2(j)
578                 trmb3(i, nsrf) = ytrmb3(j)
579            END DO            END DO
        END IF  
580    
        IF (nsrf == is_sic) THEN  
581            DO j = 1, knon            DO j = 1, knon
582               i = ni(j)               DO k = 1, klev + 1
583               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
584               IF (pctsrf_new(i, is_sic)>epsfra) THEN                  q2(i, k, nsrf) = yq2(j, k)
585                  flux_g(i) = y_flux_g(j)               END DO
              ELSE  
                 flux_g(i) = 0.  
              END IF  
586            END DO            END DO
587           else
588         END IF            fsnow(:, nsrf) = 0.
589         IF (ocean == 'slab  ') THEN         end IF if_knon
           IF (nsrf == is_oce) THEN  
              tslab(1:klon) = ytslab(1:klon)  
              seaice(1:klon) = y_seaice(1:klon)  
           END IF  
        END IF  
590      END DO loop_surface      END DO loop_surface
591    
592      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
593        frugs(:, is_oce) = rugmer
594        pctsrf(:, is_oce) = pctsrf_new_oce
595        pctsrf(:, is_sic) = pctsrf_new_sic
596    
597      rugos(:, is_oce) = rugmer      firstcal = .false.
     pctsrf = pctsrf_new  
598    
599    END SUBROUTINE clmain    END SUBROUTINE clmain
600    

Legend:
Removed from v.57  
changed lines
  Added in v.222

  ViewVC Help
Powered by ViewVC 1.1.21