/[lmdze]/trunk/phylmd/Interface_surf/pbl_surface.f
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/pbl_surface.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 186 by guez, Mon Mar 21 15:36:26 2016 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,&         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)
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 dimphy, ONLY: klev, klon, zmasq
37      ! dtime----input-R- interval du temps (secondes)      USE dimsoil, ONLY: nsoilmx
38      ! itap-----input-I- numero du pas de temps      use hbtm_m, only: hbtm
39      ! date0----input-R- jour initial      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
40      ! t--------input-R- temperature (K)      use stdlevvar_m, only: stdlevvar
41      ! q--------input-R- vapeur d'eau (kg/kg)      USE suphec_m, ONLY: rd, rg, rkappa
42      ! u--------input-R- vitesse u      use ustarhb_m, only: ustarhb
43      ! v--------input-R- vitesse v      use vdif_kcay_m, only: vdif_kcay
44      ! ts-------input-R- temperature du sol (en Kelvin)      use yamada4_m, only: yamada4
45      ! paprs----input-R- pression a intercouche (Pa)  
46      ! pplay----input-R- pression au milieu de couche (Pa)      REAL, INTENT(IN):: dtime ! interval du temps (secondes)
47      ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2      INTEGER, INTENT(IN):: itap ! numero du pas de temps
48      ! rlat-----input-R- latitude en degree      REAL, INTENT(inout):: pctsrf(klon, nbsrf)
49    
50        ! la nouvelle repartition des surfaces sortie de l'interface
51        REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
52    
53        REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
54        REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
55        REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
56        INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
57        REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal    
58        REAL, INTENT(IN):: ts(klon, nbsrf) ! temperature du sol (en Kelvin)
59        REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
60        REAL, INTENT(IN):: ksta, ksta_ter
61        LOGICAL, INTENT(IN):: ok_kzmin
62    
63        REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
64        ! soil temperature of surface fraction
65    
66        REAL, INTENT(inout):: qsol(klon)
67        ! column-density of water in soil, in kg m-2
68    
69        REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
70        REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
71        REAL snow(klon, nbsrf)
72        REAL qsurf(klon, nbsrf)
73        REAL evap(klon, nbsrf)
74        REAL, intent(inout):: falbe(klon, nbsrf)
75    
76        REAL fluxlat(klon, nbsrf)
77    
78        REAL, intent(in):: rain_fall(klon)
79        ! liquid water mass flux (kg/m2/s), positive down
80    
81        REAL, intent(in):: snow_f(klon)
82        ! solid water mass flux (kg/m2/s), positive down
83    
84        REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
85        REAL, intent(in):: fder(klon)
86        REAL, INTENT(IN):: rlat(klon) ! latitude en degr\'es
87    
88        REAL rugos(klon, nbsrf)
89      ! 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)  
90    
91        LOGICAL, INTENT(IN):: debut
92        real agesno(klon, nbsrf)
93        REAL, INTENT(IN):: rugoro(klon)
94    
95        REAL d_t(klon, klev), d_q(klon, klev)
96      ! d_t------output-R- le changement pour "t"      ! d_t------output-R- le changement pour "t"
97      ! d_q------output-R- le changement pour "q"      ! d_q------output-R- le changement pour "q"
98      ! d_u------output-R- le changement pour "u"  
99      ! d_v------output-R- le changement pour "v"      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
100      ! d_ts-----output-R- le changement pour "ts"      ! changement pour "u" et "v"
101    
102        REAL, intent(out):: d_ts(klon, nbsrf) ! le changement pour "ts"
103    
104        REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
105      ! 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)
106      !                    (orientation positive vers le bas)      !                    (orientation positive vers le bas)
107      ! 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)
108    
109        REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
110      ! 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
111      ! 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
112    
113        REAL, INTENT(out):: cdragh(klon), cdragm(klon)
114        real q2(klon, klev+1, nbsrf)
115    
116        REAL, INTENT(out):: dflux_t(klon), dflux_q(klon)
117      ! dflux_t derive du flux sensible      ! dflux_t derive du flux sensible
118      ! dflux_q derive du flux latent      ! dflux_q derive du flux latent
119      !IM "slab" ocean      !IM "slab" ocean
     ! flux_g---output-R-  flux glace (pour OCEAN='slab  ')  
     ! flux_o---output-R-  flux ocean (pour OCEAN='slab  ')  
120    
121      ! tslab-in/output-R temperature du slab ocean (en Kelvin)      REAL, intent(out):: ycoefh(klon, klev)
122      ! uniqmnt pour slab      REAL, intent(out):: zu1(klon)
123        REAL zv1(klon)
124        REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
125        REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
126    
127      ! seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      ! Ionela Musat cf. Anne Mathieu : pbl, hbtm (Comme les autres
128      !cc      ! diagnostics on cumule dans physiq ce qui permet de sortir les
129      ! ffonte----Flux thermique utilise pour fondre la neige      ! grandeurs par sous-surface)
130      ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      REAL pblh(klon, nbsrf)
131      !           hauteur de neige, en kg/m2/s      ! pblh------- HCL
132      ! on rajoute en output yu1 et yv1 qui sont les vents dans      REAL capcl(klon, nbsrf)
133      ! la premiere couche      REAL oliqcl(klon, nbsrf)
134      ! ces 4 variables sont maintenant traites dans phytrac      REAL cteicl(klon, nbsrf)
135      ! itr--------input-I- nombre de traceurs      REAL pblt(klon, nbsrf)
136      ! tr---------input-R- q. de traceurs      ! pblT------- T au nveau HCL
137      ! flux_surf--input-R- flux de traceurs a la surface      REAL therm(klon, nbsrf)
138      ! d_tr-------output-R tendance de traceurs      REAL trmb1(klon, nbsrf)
     !IM cf. AM : PBL  
139      ! trmb1-------deep_cape      ! trmb1-------deep_cape
140        REAL trmb2(klon, nbsrf)
141      ! trmb2--------inhibition      ! trmb2--------inhibition
142        REAL trmb3(klon, nbsrf)
143      ! trmb3-------Point Omega      ! trmb3-------Point Omega
144      ! 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)  
145      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)      REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
146      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      ! ffonte----Flux thermique utilise pour fondre la neige
147        ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
148      REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)      !           hauteur de neige, en kg/m2/s
149      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)  
   
     REAL rain_f(klon), snow_f(klon)  
     REAL fder(klon)  
   
     REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)  
     REAL rugos(klon, nbsrf)  
     ! la nouvelle repartition des surfaces sortie de l'interface  
     REAL pctsrf_new(klon, nbsrf)  
150    
151      REAL zcoefh(klon, klev)      ! Local:
     REAL zu1(klon)  
     REAL zv1(klon)  
152    
153      !$$$ PB ajout pour soil      REAL y_fqcalving(klon), y_ffonte(klon)
154      LOGICAL, INTENT (IN) :: soil_model      real y_run_off_lic_0(klon)
     !IM ajout seuils cdrm, cdrh  
     REAL cdmmax, cdhmax  
155    
156      REAL ksta, ksta_ter      REAL rugmer(klon)
     LOGICAL ok_kzmin  
157    
     REAL ftsoil(klon, nsoilmx, nbsrf)  
158      REAL ytsoil(klon, nsoilmx)      REAL ytsoil(klon, nsoilmx)
     REAL qsol(klon)  
   
     EXTERNAL clqh, clvent, coefkz, calbeta, cltrac  
159    
160      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
161      REAL yalb(klon)      REAL yalb(klon)
     REAL yalblw(klon)  
162      REAL yu1(klon), yv1(klon)      REAL yu1(klon), yv1(klon)
163      REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      ! on rajoute en output yu1 et yv1 qui sont les vents dans
164      REAL yrain_f(klon), ysnow_f(klon)      ! la premiere couche
165      REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL ysnow(klon), yqsurf(klon), yagesno(klon)
166      REAL yfder(klon), ytaux(klon), ytauy(klon)  
167        real yqsol(klon)
168        ! column-density of water in soil, in kg m-2
169    
170        REAL yrain_f(klon)
171        ! liquid water mass flux (kg/m2/s), positive down
172    
173        REAL ysnow_f(klon)
174        ! solid water mass flux (kg/m2/s), positive down
175    
176        REAL yfder(klon)
177      REAL yrugm(klon), yrads(klon), yrugoro(klon)      REAL yrugm(klon), yrads(klon), yrugoro(klon)
178    
179      REAL yfluxlat(klon)      REAL yfluxlat(klon)
# Line 199  contains Line 184  contains
184      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)      REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
185      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)      REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
186      REAL y_dflux_t(klon), y_dflux_q(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
187      REAL ycoefh(klon, klev), ycoefm(klon, klev)      REAL coefh(klon, klev), coefm(klon, klev)
188      REAL yu(klon, klev), yv(klon, klev)      REAL yu(klon, klev), yv(klon, klev)
189      REAL yt(klon, klev), yq(klon, klev)      REAL yt(klon, klev), yq(klon, klev)
190      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)      REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
191    
     LOGICAL ok_nonloc  
     PARAMETER (ok_nonloc=.FALSE.)  
192      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)      REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
193    
     !IM 081204 hcl_Anne ? BEG  
194      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)      REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
195      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)      REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
196      REAL ykmq(klon, klev+1)      REAL ykmq(klon, klev+1)
197      REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)      REAL yq2(klon, klev+1)
198      REAL q2diag(klon, klev+1)      REAL q2diag(klon, klev+1)
     !IM 081204 hcl_Anne ? END  
199    
200      REAL u1lay(klon), v1lay(klon)      REAL u1lay(klon), v1lay(klon)
201      REAL delp(klon, klev)      REAL delp(klon, klev)
202      INTEGER i, k, nsrf      INTEGER i, k, nsrf
203    
204      INTEGER ni(klon), knon, j      INTEGER ni(klon), knon, j
205      ! Introduction d'une variable "pourcentage potentiel" pour tenir compte  
     ! des eventuelles apparitions et/ou disparitions de la glace de mer  
206      REAL pctsrf_pot(klon, nbsrf)      REAL pctsrf_pot(klon, nbsrf)
207        ! "pourcentage potentiel" pour tenir compte des \'eventuelles
208        ! apparitions ou disparitions de la glace de mer
209    
210      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
211    
     ! 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)  
   
212      REAL yt2m(klon), yq2m(klon), yu10m(klon)      REAL yt2m(klon), yq2m(klon), yu10m(klon)
213      REAL yustar(klon)      REAL yustar(klon)
     ! -- LOOP  
     REAL yu10mx(klon)  
     REAL yu10my(klon)  
     REAL ywindsp(klon)  
     ! -- LOOP  
214    
215      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)  
216      REAL ypblh(klon)      REAL ypblh(klon)
217      REAL ylcl(klon)      REAL ylcl(klon)
218      REAL ycapcl(klon)      REAL ycapcl(klon)
# Line 278  contains Line 223  contains
223      REAL ytrmb1(klon)      REAL ytrmb1(klon)
224      REAL ytrmb2(klon)      REAL ytrmb2(klon)
225      REAL ytrmb3(klon)      REAL ytrmb3(klon)
     REAL y_cd_h(klon), y_cd_m(klon)  
226      REAL uzon(klon), vmer(klon)      REAL uzon(klon), vmer(klon)
227      REAL tair1(klon), qair1(klon), tairsol(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
228      REAL psfce(klon), patm(klon)      REAL psfce(klon), patm(klon)
# Line 290  contains Line 234  contains
234      LOGICAL zxli      LOGICAL zxli
235      PARAMETER (zxli=.FALSE.)      PARAMETER (zxli=.FALSE.)
236    
     REAL zt, zqs, zdelta, zcor  
     REAL t_coup  
     PARAMETER (t_coup=273.15)  
   
     CHARACTER (len=20) :: modname = 'clmain'  
   
237      !------------------------------------------------------------      !------------------------------------------------------------
238    
     ! initialisation Anne  
239      ytherm = 0.      ytherm = 0.
240    
     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  
   
241      DO k = 1, klev ! epaisseur de couche      DO k = 1, klev ! epaisseur de couche
242         DO i = 1, klon         DO i = 1, klon
243            delp(i, k) = paprs(i, k) - paprs(i, k+1)            delp(i, k) = paprs(i, k) - paprs(i, k+1)
# Line 342  contains Line 250  contains
250         v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2         v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
251      END DO      END DO
252    
253      ! initialisation:      ! Initialization:
254        rugmer = 0.
255      DO i = 1, klon      cdragh = 0.
256         rugmer(i) = 0.0      cdragm = 0.
257         cdragh(i) = 0.0      dflux_t = 0.
258         cdragm(i) = 0.0      dflux_q = 0.
259         dflux_t(i) = 0.0      zu1 = 0.
260         dflux_q(i) = 0.0      zv1 = 0.
261         zu1(i) = 0.0      ypct = 0.
262         zv1(i) = 0.0      yts = 0.
263      END DO      ysnow = 0.
264      ypct = 0.0      yqsurf = 0.
265      yts = 0.0      yrain_f = 0.
266      ysnow = 0.0      ysnow_f = 0.
267      yqsurf = 0.0      yfder = 0.
268      yalb = 0.0      yrugos = 0.
269      yalblw = 0.0      yu1 = 0.
270      yrain_f = 0.0      yv1 = 0.
271      ysnow_f = 0.0      yrads = 0.
272      yfder = 0.0      ypaprs = 0.
273      ytaux = 0.0      ypplay = 0.
274      ytauy = 0.0      ydelp = 0.
275      ysolsw = 0.0      yu = 0.
276      ysollw = 0.0      yv = 0.
277      ysollwdown = 0.0      yt = 0.
278      yrugos = 0.0      yq = 0.
279      yu1 = 0.0      pctsrf_new = 0.
280      yv1 = 0.0      y_flux_u = 0.
281      yrads = 0.0      y_flux_v = 0.
282      ypaprs = 0.0      y_dflux_t = 0.
283      ypplay = 0.0      y_dflux_q = 0.
     ydelp = 0.0  
     yu = 0.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  
284      ytsoil = 999999.      ytsoil = 999999.
285      yrugoro = 0.      yrugoro = 0.
286      ! -- LOOP      d_ts = 0.
     yu10mx = 0.0  
     yu10my = 0.0  
     ywindsp = 0.0  
     ! -- LOOP  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           d_ts(i, nsrf) = 0.0  
        END DO  
     END DO  
     !§§§ PB  
287      yfluxlat = 0.      yfluxlat = 0.
288      flux_t = 0.      flux_t = 0.
289      flux_q = 0.      flux_q = 0.
290      flux_u = 0.      flux_u = 0.
291      flux_v = 0.      flux_v = 0.
292      DO k = 1, klev      d_t = 0.
293         DO i = 1, klon      d_q = 0.
294            d_t(i, k) = 0.0      d_u = 0.
295            d_q(i, k) = 0.0      d_v = 0.
296            d_u(i, k) = 0.0      ycoefh = 0.
297            d_v(i, k) = 0.0  
298            zcoefh(i, k) = 0.0      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
299         END DO      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
300      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)  
301    
302      pctsrf_pot = pctsrf      pctsrf_pot = pctsrf
303      pctsrf_pot(:, is_oce) = 1. - zmasq      pctsrf_pot(:, is_oce) = 1. - zmasq
304      pctsrf_pot(:, is_sic) = 1. - zmasq      pctsrf_pot(:, is_sic) = 1. - zmasq
305    
306      DO nsrf = 1, nbsrf      ! Boucler sur toutes les sous-fractions du sol:
307         ! chercher les indices:  
308        loop_surface: DO nsrf = 1, nbsrf
309           ! Chercher les indices :
310         ni = 0         ni = 0
311         knon = 0         knon = 0
312         DO i = 1, klon         DO i = 1, klon
313            ! pour determiner le domaine a traiter on utilise les surfaces            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
314            ! "potentielles"            ! "potentielles"
315            IF (pctsrf_pot(i, nsrf) > epsfra) THEN            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
316               knon = knon + 1               knon = knon + 1
# Line 435  contains Line 318  contains
318            END IF            END IF
319         END DO         END DO
320    
321         ! 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  
322            DO j = 1, knon            DO j = 1, knon
323               i = ni(j)               i = ni(j)
324               ytsoil(j, k) = ftsoil(i, k, nsrf)               ypct(j) = pctsrf(i, nsrf)
325            END DO               yts(j) = ts(i, nsrf)
326         END DO               ysnow(j) = snow(i, nsrf)
327         DO k = 1, klev               yqsurf(j) = qsurf(i, nsrf)
328            DO j = 1, knon               yalb(j) = falbe(i, nsrf)
329               i = ni(j)               yrain_f(j) = rain_fall(i)
330               ypaprs(j, k) = paprs(i, k)               ysnow_f(j) = snow_f(i)
331               ypplay(j, k) = pplay(i, k)               yagesno(j) = agesno(i, nsrf)
332               ydelp(j, k) = delp(i, k)               yfder(j) = fder(i)
333               yu(j, k) = u(i, k)               yrugos(j) = rugos(i, nsrf)
334               yv(j, k) = v(i, k)               yrugoro(j) = rugoro(i)
335               yt(j, k) = t(i, k)               yu1(j) = u1lay(i)
336               yq(j, k) = q(i, k)               yv1(j) = v1lay(i)
337            END DO               yrads(j) = solsw(i, nsrf) + sollw(i, nsrf)
338         END DO               ypaprs(j, klev+1) = paprs(i, klev+1)
339                 y_run_off_lic_0(j) = run_off_lic_0(i)
340              END DO
341    
342              ! For continent, copy soil water content
343              IF (nsrf == is_ter) THEN
344                 yqsol(:knon) = qsol(ni(:knon))
345              ELSE
346                 yqsol = 0.
347              END IF
348    
349         ! calculer Cdrag et les coefficients d'echange            DO k = 1, nsoilmx
350         CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&               DO j = 1, knon
351              yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)                  i = ni(j)
352         !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))  
353               END DO               END DO
354            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)  
355    
           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  
356            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  
357               DO j = 1, knon               DO j = 1, knon
358                  i = ni(j)                  i = ni(j)
359                  yq2(j, k) = q2(i, k, nsrf)                  ypaprs(j, k) = paprs(i, k)
360                    ypplay(j, k) = pplay(i, k)
361                    ydelp(j, k) = delp(i, k)
362                    yu(j, k) = u(i, k)
363                    yv(j, k) = v(i, k)
364                    yt(j, k) = t(i, k)
365                    yq(j, k) = q(i, k)
366               END DO               END DO
367            END DO            END DO
368    
369            !   Bug introduit volontairement pour converger avec les resultats            ! calculer Cdrag et les coefficients d'echange
370            !  du papier sur les thermiques.            CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
371            IF (1==1) THEN                 yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
372               y_cd_m(1:knon) = ycoefm(1:knon, 1)            IF (iflag_pbl == 1) THEN
373               y_cd_h(1:knon) = ycoefh(1:knon, 1)               CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
374            ELSE               coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
375               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)  
376            END IF            END IF
           CALL ustarhb(knon, yu, yv, y_cd_m, yustar)  
377    
378            IF (prt_level>9) THEN            ! on met un seuil pour coefm et coefh
379               PRINT *, 'USTAR = ', yustar            IF (nsrf == is_oce) THEN
380                 coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
381                 coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
382            END IF            END IF
383    
384            !   iflag_pbl peut etre utilise comme longuer de melange            IF (ok_kzmin) THEN
385                 ! Calcul d'une diffusion minimale pour les conditions tres stables
386                 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
387                      coefm(:knon, 1), ycoefm0, ycoefh0)
388                 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
389                 coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
390              END IF
391    
392            IF (iflag_pbl>=11) THEN            IF (iflag_pbl >= 3) THEN
393               CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &               ! Mellor et Yamada adapt\'e \`a Mars, Richard Fournier et
394                    yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &               ! Fr\'ed\'eric Hourdin
395                    iflag_pbl)               yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
396            ELSE                    + ypplay(:knon, 1))) &
397               CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &                    * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
398                    yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)               DO k = 2, klev
399                    yzlay(1:knon, k) = yzlay(1:knon, k-1) &
400                         + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
401                         / ypaprs(1:knon, k) &
402                         * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
403                 END DO
404                 DO k = 1, klev
405                    yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
406                         / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
407                 END DO
408                 yzlev(1:knon, 1) = 0.
409                 yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
410                      - yzlay(:knon, klev - 1)
411                 DO k = 2, klev
412                    yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
413                 END DO
414                 DO k = 1, klev + 1
415                    DO j = 1, knon
416                       i = ni(j)
417                       yq2(j, k) = q2(i, k, nsrf)
418                    END DO
419                 END DO
420    
421                 CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
422                 IF (prt_level > 9) PRINT *, 'USTAR = ', yustar
423    
424                 ! iflag_pbl peut \^etre utilis\'e comme longueur de m\'elange
425    
426                 IF (iflag_pbl >= 11) THEN
427                    CALL vdif_kcay(knon, dtime, rg, ypaprs, yzlev, yzlay, yu, yv, &
428                         yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, yustar, &
429                         iflag_pbl)
430                 ELSE
431                    CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
432                         coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
433                 END IF
434    
435                 coefm(:knon, 2:) = ykmm(:knon, 2:klev)
436                 coefh(:knon, 2:) = ykmn(:knon, 2:klev)
437            END IF            END IF
438    
439            ycoefm(1:knon, 1) = y_cd_m(1:knon)            ! calculer la diffusion des vitesses "u" et "v"
440            ycoefh(1:knon, 1) = y_cd_h(1:knon)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
441            ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)                 ypplay, ydelp, y_d_u, y_flux_u)
442            ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)            CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
443         END IF                 ypplay, ydelp, y_d_v, y_flux_v)
444    
445         ! calculer la diffusion des vitesses "u" et "v"            ! calculer la diffusion de "q" et de "h"
446         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &            CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni(:knon), &
447              ydelp, y_d_u, y_flux_u)                 pctsrf, ytsoil, yqsol, rmu0, yrugos, yrugoro, yu1, &
448         CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &                 yv1, coefh(:knon, :), yt, yq, yts, ypaprs, ypplay, ydelp, &
449              ydelp, y_d_v, y_flux_v)                 yrads, yalb(:knon), ysnow, yqsurf, yrain_f, ysnow_f, yfder, &
450                   yfluxlat, pctsrf_new, yagesno(:knon), y_d_t, y_d_q, &
451         ! pour le couplage                 y_d_ts(:knon), yz0_new, y_flux_t, y_flux_q, y_dflux_t, &
452         ytaux = y_flux_u(:, 1)                 y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0)
453         ytauy = y_flux_v(:, 1)  
454              ! calculer la longueur de rugosite sur ocean
455         ! FH modif sur le cdrag temperature            yrugm = 0.
456         !$$$PB : déplace dans clcdrag            IF (nsrf == is_oce) THEN
457         !$$$      do i=1, knon               DO j = 1, knon
458         !$$$         ycoefh(i, 1)=ycoefm(i, 1)*0.8                  yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
459         !$$$      enddo                       0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
460                    yrugm(j) = max(1.5E-05, yrugm(j))
461         ! calculer la diffusion de "q" et de "h"               END DO
462         CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&            END IF
             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  
463            DO j = 1, knon            DO j = 1, knon
464               yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &               y_dflux_t(j) = y_dflux_t(j)*ypct(j)
465                    0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))               y_dflux_q(j) = y_dflux_q(j)*ypct(j)
466               yrugm(j) = max(1.5E-05, yrugm(j))               yu1(j) = yu1(j)*ypct(j)
467            END DO               yv1(j) = yv1(j)*ypct(j)
468         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  
469    
470         DO k = 1, klev            DO k = 1, klev
471            DO j = 1, knon               DO j = 1, knon
472               i = ni(j)                  i = ni(j)
473               ycoefh(j, k) = ycoefh(j, k)*ypct(j)                  coefh(j, k) = coefh(j, k)*ypct(j)
474               ycoefm(j, k) = ycoefm(j, k)*ypct(j)                  coefm(j, k) = coefm(j, k)*ypct(j)
475               y_d_t(j, k) = y_d_t(j, k)*ypct(j)                  y_d_t(j, k) = y_d_t(j, k)*ypct(j)
476               y_d_q(j, k) = y_d_q(j, k)*ypct(j)                  y_d_q(j, k) = y_d_q(j, k)*ypct(j)
477               !§§§ PB                  flux_t(i, k, nsrf) = y_flux_t(j, k)
478               flux_t(i, k, nsrf) = y_flux_t(j, k)                  flux_q(i, k, nsrf) = y_flux_q(j, k)
479               flux_q(i, k, nsrf) = y_flux_q(j, k)                  flux_u(i, k, nsrf) = y_flux_u(j, k)
480               flux_u(i, k, nsrf) = y_flux_u(j, k)                  flux_v(i, k, nsrf) = y_flux_v(j, k)
481               flux_v(i, k, nsrf) = y_flux_v(j, k)                  y_d_u(j, k) = y_d_u(j, k)*ypct(j)
482               !$$$ 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)
483               !$$$ 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)  
484            END DO            END DO
        END DO  
485    
486         evap(:, nsrf) = -flux_q(:, 1, nsrf)            evap(:, nsrf) = -flux_q(:, 1, nsrf)
487    
488         albe(:, nsrf) = 0.            falbe(:, nsrf) = 0.
489         alblw(:, nsrf) = 0.            snow(:, nsrf) = 0.
490         snow(:, nsrf) = 0.            qsurf(:, nsrf) = 0.
491         qsurf(:, nsrf) = 0.            rugos(:, nsrf) = 0.
492         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  
493            DO j = 1, knon            DO j = 1, knon
494               i = ni(j)               i = ni(j)
495               qsol(i) = yqsol(j)               d_ts(i, nsrf) = y_d_ts(j)
496                 falbe(i, nsrf) = yalb(j)
497                 snow(i, nsrf) = ysnow(j)
498                 qsurf(i, nsrf) = yqsurf(j)
499                 rugos(i, nsrf) = yz0_new(j)
500                 fluxlat(i, nsrf) = yfluxlat(j)
501                 IF (nsrf == is_oce) THEN
502                    rugmer(i) = yrugm(j)
503                    rugos(i, nsrf) = yrugm(j)
504                 END IF
505                 agesno(i, nsrf) = yagesno(j)
506                 fqcalving(i, nsrf) = y_fqcalving(j)
507                 ffonte(i, nsrf) = y_ffonte(j)
508                 cdragh(i) = cdragh(i) + coefh(j, 1)
509                 cdragm(i) = cdragm(i) + coefm(j, 1)
510                 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
511                 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
512                 zu1(i) = zu1(i) + yu1(j)
513                 zv1(i) = zv1(i) + yv1(j)
514              END DO
515              IF (nsrf == is_ter) THEN
516                 qsol(ni(:knon)) = yqsol(:knon)
517              else IF (nsrf == is_lic) THEN
518                 DO j = 1, knon
519                    i = ni(j)
520                    run_off_lic_0(i) = y_run_off_lic_0(j)
521                 END DO
522              END IF
523    
524              ftsoil(:, :, nsrf) = 0.
525              DO k = 1, nsoilmx
526                 DO j = 1, knon
527                    i = ni(j)
528                    ftsoil(i, k, nsrf) = ytsoil(j, k)
529                 END DO
530            END DO            END DO
531         END IF  
        IF (nsrf==is_lic) THEN  
532            DO j = 1, knon            DO j = 1, knon
533               i = ni(j)               i = ni(j)
534               run_off_lic_0(i) = y_run_off_lic_0(j)               DO k = 1, klev
535                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
536                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
537                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
538                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
539                    ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
540                 END DO
541            END DO            END DO
542         END IF  
543         !$$$ PB ajout pour soil            ! diagnostic t, q a 2m et u, v a 10m
544         ftsoil(:, :, nsrf) = 0.  
        DO k = 1, nsoilmx  
545            DO j = 1, knon            DO j = 1, knon
546               i = ni(j)               i = ni(j)
547               ftsoil(i, k, nsrf) = ytsoil(j, k)               uzon(j) = yu(j, 1) + y_d_u(j, 1)
548            END DO               vmer(j) = yv(j, 1) + y_d_v(j, 1)
549         END DO               tair1(j) = yt(j, 1) + y_d_t(j, 1)
550                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
551                 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
552                      1)))*(ypaprs(j, 1)-ypplay(j, 1))
553                 tairsol(j) = yts(j) + y_d_ts(j)
554                 rugo1(j) = yrugos(j)
555                 IF (nsrf == is_oce) THEN
556                    rugo1(j) = rugos(i, nsrf)
557                 END IF
558                 psfce(j) = ypaprs(j, 1)
559                 patm(j) = ypplay(j, 1)
560    
561         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)  
562            END DO            END DO
        END DO  
563    
564         !cc diagnostic t, q a 2m et u, v a 10m            CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
565                   zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
566                   yt10m, yq10m, yu10m, yustar)
567    
568         DO j = 1, knon            DO j = 1, knon
569            i = ni(j)               i = ni(j)
570            uzon(j) = yu(j, 1) + y_d_u(j, 1)               t2m(i, nsrf) = yt2m(j)
571            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)  
   
           qairsol(j) = yqsurf(j)  
        END DO  
572    
573         CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &               ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
574              tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &               u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
575              yu10m, yustar)               v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
        !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)  
576    
577         END DO            END DO
578    
579         DO i = 1, knon            CALL hbtm(knon, ypaprs, ypplay, yt2m, yq2m, yustar, y_flux_t, &
580            y_cd_h(i) = ycoefh(i, 1)                 y_flux_q, yu, yv, yt, yq, ypblh(:knon), ycapcl, yoliqcl, &
581            y_cd_m(i) = ycoefm(i, 1)                 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  
582    
        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  
583            DO j = 1, knon            DO j = 1, knon
              ! on projette sur la grille globale  
584               i = ni(j)               i = ni(j)
585               IF (pctsrf_new(i, is_oce)>epsfra) THEN               pblh(i, nsrf) = ypblh(j)
586                  flux_o(i) = y_flux_o(j)               plcl(i, nsrf) = ylcl(j)
587               ELSE               capcl(i, nsrf) = ycapcl(j)
588                  flux_o(i) = 0.               oliqcl(i, nsrf) = yoliqcl(j)
589               END IF               cteicl(i, nsrf) = ycteicl(j)
590                 pblt(i, nsrf) = ypblt(j)
591                 therm(i, nsrf) = ytherm(j)
592                 trmb1(i, nsrf) = ytrmb1(j)
593                 trmb2(i, nsrf) = ytrmb2(j)
594                 trmb3(i, nsrf) = ytrmb3(j)
595            END DO            END DO
        END IF  
596    
        IF (nsrf==is_sic) THEN  
597            DO j = 1, knon            DO j = 1, knon
598               i = ni(j)               DO k = 1, klev + 1
599               ! On pondère lorsque l'on fait le bilan au sol :                  i = ni(j)
600               ! flux_g(i) = y_flux_g(j)*ypct(j)                  q2(i, k, nsrf) = yq2(j, k)
601               IF (pctsrf_new(i, is_sic)>epsfra) THEN               END DO
                 flux_g(i) = y_flux_g(j)  
              ELSE  
                 flux_g(i) = 0.  
              END IF  
602            END DO            END DO
603           end IF if_knon
604         END IF      END DO loop_surface
        !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                                                        
           END IF  
           !OCEAN                                                        
        END IF  
     END DO  
605    
606      ! On utilise les nouvelles surfaces      ! On utilise les nouvelles surfaces
     ! A rajouter: conservation de l'albedo  
607    
608      rugos(:, is_oce) = rugmer      rugos(:, is_oce) = rugmer
609      pctsrf = pctsrf_new      pctsrf = pctsrf_new

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

  ViewVC Help
Powered by ViewVC 1.1.21