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

Legend:
Removed from v.38  
changed lines
  Added in v.155

  ViewVC Help
Powered by ViewVC 1.1.21