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

Diff of /trunk/phylmd/clqh.f

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

trunk/libf/phylmd/clqh.f revision 12 by guez, Mon Jul 21 16:05:07 2008 UTC trunk/libf/phylmd/clqh.f90 revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC
# Line 1  Line 1 
1        SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin,  SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin, &
2       e                rlon, rlat, cufi, cvfi,       rlon, rlat, cufi, cvfi,  &
3       e                knon, nisurf, knindex, pctsrf,       knon, nisurf, knindex, pctsrf, &
4       $                soil_model,tsoil,qsol,       soil_model,tsoil,qsol, &
5       e                ok_veget, ocean, npas, nexca,       ok_veget, ocean, npas, nexca, &
6       e                rmu0, co2_ppm, rugos, rugoro,       rmu0, co2_ppm, rugos, rugoro, &
7       e                u1lay,v1lay,coef,       u1lay,v1lay,coef, &
8       e                t,q,ts,paprs,pplay,       t,q,ts,paprs,pplay, &
9       e                delp,radsol,albedo,alblw,snow,qsurf,       delp,radsol,albedo,alblw,snow,qsurf,  &
10       e                precip_rain, precip_snow, fder, taux, tauy,       precip_rain, precip_snow, fder, taux, tauy, ywindsp, &
11  c -- LOOP       sollw, sollwdown, swnet,fluxlat,  &
12       e                ywindsp,       pctsrf_new, agesno, &
13  c -- LOOP       d_t, d_q, d_ts, z0_new,  &
14       $                sollw, sollwdown, swnet,fluxlat,       flux_t, flux_q,dflux_s,dflux_l, &
15       s                pctsrf_new, agesno,       fqcalving,ffonte,run_off_lic_0, &
16       s                d_t, d_q, d_ts, z0_new,       flux_o,flux_g,tslab,seaice)
17       s                flux_t, flux_q,dflux_s,dflux_l,  
18       s                fqcalving,ffonte,run_off_lic_0,    use conf_phys_m
19  cIM "slab" ocean    use dimens_m
20       s                flux_o,flux_g,tslab,seaice)    use dimphy
21      use dimsoil
22        USE interface_surf    use fcttre
23      use indicesol
24        use dimens_m    USE interface_surf
25        use indicesol    use iniprint
26        use dimphy    use suphec_m, only: rcpd, rd, rg, rkappa
27        use dimsoil    use YOMCST
28        use iniprint    use yoethf
29        use YOMCST  
30        use yoethf    IMPLICIT none
31        use fcttre  
32        use conf_phys_m    ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
33        IMPLICIT none    ! Objet: diffusion verticale de "q" et de "h"
34  c======================================================================  
35  c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818    ! Arguments:
36  c Objet: diffusion verticale de "q" et de "h"    INTEGER knon
37  c======================================================================    REAL, intent(in):: dtime              ! intervalle du temps (s)
38      real date0
39  c Arguments:    REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
40        INTEGER knon    REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
41        REAL, intent(in):: dtime              ! intervalle du temps (s)    REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
42        real date0    !                               multiplie par le cisaillement du
43        REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)    !                               vent (dV/dz); la premiere valeur
44        REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)    !                               indique la valeur de Cdrag (sans unite)
45        REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)    REAL t(klon,klev)       ! temperature (K)
46  c                               multiplie par le cisaillement du    REAL q(klon,klev)       ! humidite specifique (kg/kg)
47  c                               vent (dV/dz); la premiere valeur    REAL ts(klon)           ! temperature du sol (K)
48  c                               indique la valeur de Cdrag (sans unite)    REAL evap(klon)         ! evaporation au sol
49        REAL t(klon,klev)       ! temperature (K)    REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
50        REAL q(klon,klev)       ! humidite specifique (kg/kg)    REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
51        REAL ts(klon)           ! temperature du sol (K)    REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
52        REAL evap(klon)         ! evaporation au sol    REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
53        REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)    REAL albedo(klon)       ! albedo de la surface
54        REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)    REAL alblw(klon)
55        REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)    REAL snow(klon)         ! hauteur de neige
56        REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2    REAL qsurf(klon)         ! humidite de l'air au dessus de la surface
57        REAL albedo(klon)       ! albedo de la surface    real precip_rain(klon), precip_snow(klon)
58        REAL alblw(klon)    REAL agesno(klon)
59        REAL snow(klon)         ! hauteur de neige    REAL rugoro(klon)
60        REAL qsurf(klon)         ! humidite de l'air au dessus de la surface    REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
61        real precip_rain(klon), precip_snow(klon)    integer jour            ! jour de l'annee en cours
62        REAL agesno(klon)    real rmu0(klon)         ! cosinus de l'angle solaire zenithal
63        REAL rugoro(klon)    real rugos(klon)        ! rugosite
64        REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent    integer knindex(klon)
65        integer jour            ! jour de l'annee en cours    real pctsrf(klon,nbsrf)
66        real rmu0(klon)         ! cosinus de l'angle solaire zenithal    real, intent(in):: rlon(klon), rlat(klon)
67        real rugos(klon)        ! rugosite    real cufi(klon), cvfi(klon)
68        integer knindex(klon)    logical ok_veget
69        real pctsrf(klon,nbsrf)    REAL co2_ppm            ! taux CO2 atmosphere
70        real, intent(in):: rlon(klon), rlat(klon)    character(len=*), intent(in):: ocean
71        real cufi(klon), cvfi(klon)    integer npas, nexca
72        logical ok_veget    ! -- LOOP
73        REAL co2_ppm            ! taux CO2 atmosphere    REAL yu10mx(klon)
74        character(len=*), intent(in):: ocean    REAL yu10my(klon)
75        integer npas, nexca    REAL ywindsp(klon)
76  c -- LOOP    ! -- LOOP
77         REAL yu10mx(klon)  
78         REAL yu10my(klon)  
79         REAL ywindsp(klon)    !
80  c -- LOOP    REAL d_t(klon,klev)     ! incrementation de "t"
81      REAL d_q(klon,klev)     ! incrementation de "q"
82      REAL d_ts(klon)         ! incrementation de "ts"
83  c    REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
84        REAL d_t(klon,klev)     ! incrementation de "t"    !                               sensible, flux de Cp*T, positif vers
85        REAL d_q(klon,klev)     ! incrementation de "q"    !                               le bas: j/(m**2 s) c.a.d.: W/m2
86        REAL d_ts(klon)         ! incrementation de "ts"    REAL flux_q(klon,klev)  ! flux de la vapeur d'eau:kg/(m**2 s)
87        REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur    REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
88  c                               sensible, flux de Cp*T, positif vers    REAL dflux_l(klon) ! derivee du flux latent dF/dTs
89  c                               le bas: j/(m**2 s) c.a.d.: W/m2    !IM cf JLD
90        REAL flux_q(klon,klev)  ! flux de la vapeur d'eau:kg/(m**2 s)    ! Flux thermique utiliser pour fondre la neige
91        REAL dflux_s(klon) ! derivee du flux sensible dF/dTs    REAL ffonte(klon)
92        REAL dflux_l(klon) ! derivee du flux latent dF/dTs    ! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
93  cIM cf JLD    ! hauteur de neige, en kg/m2/s
94  c Flux thermique utiliser pour fondre la neige    REAL fqcalving(klon)
95        REAL ffonte(klon)    !IM "slab" ocean
96  c Flux d'eau "perdue" par la surface et nécessaire pour que limiter la    REAL tslab(klon)  !temperature du slab ocean (K) (OCEAN='slab  ')
97  c hauteur de neige, en kg/m2/s    REAL seaice(klon) ! glace de mer en kg/m2
98        REAL fqcalving(klon)    REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2
99  cIM "slab" ocean    REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2
100        REAL tslab(klon)  !temperature du slab ocean (K) (OCEAN='slab  ')    !
101        REAL seaice(klon) ! glace de mer en kg/m2    !======================================================================
102        REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2    REAL t_grnd  ! temperature de rappel pour glace de mer
103        REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2    PARAMETER (t_grnd=271.35)
104  c    REAL t_coup
105  c======================================================================    PARAMETER(t_coup=273.15)
106        REAL t_grnd  ! temperature de rappel pour glace de mer    !======================================================================
107        PARAMETER (t_grnd=271.35)    INTEGER i, k
108        REAL t_coup    REAL zx_cq(klon,klev)
109        PARAMETER(t_coup=273.15)    REAL zx_dq(klon,klev)
110  c======================================================================    REAL zx_ch(klon,klev)
111        INTEGER i, k    REAL zx_dh(klon,klev)
112        REAL zx_cq(klon,klev)    REAL zx_buf1(klon)
113        REAL zx_dq(klon,klev)    REAL zx_buf2(klon)
114        REAL zx_ch(klon,klev)    REAL zx_coef(klon,klev)
115        REAL zx_dh(klon,klev)    REAL local_h(klon,klev) ! enthalpie potentielle
116        REAL zx_buf1(klon)    REAL local_q(klon,klev)
117        REAL zx_buf2(klon)    REAL local_ts(klon)
118        REAL zx_coef(klon,klev)    REAL psref(klon) ! pression de reference pour temperature potent.
119        REAL local_h(klon,klev) ! enthalpie potentielle    REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
120        REAL local_q(klon,klev)    !======================================================================
121        REAL local_ts(klon)    ! contre-gradient pour la vapeur d'eau: (kg/kg)/metre
122        REAL psref(klon) ! pression de reference pour temperature potent.    REAL gamq(klon,2:klev)
123        REAL zx_pkh(klon,klev), zx_pkf(klon,klev)    ! contre-gradient pour la chaleur sensible: Kelvin/metre
124  c======================================================================    REAL gamt(klon,2:klev)
125  c contre-gradient pour la vapeur d'eau: (kg/kg)/metre    REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
126        REAL gamq(klon,2:klev)    REAL zdelz
127  c contre-gradient pour la chaleur sensible: Kelvin/metre    !======================================================================
128        REAL gamt(klon,2:klev)    !======================================================================
129        REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)    ! Rajout pour l'interface
130        REAL zdelz    integer, intent(in):: itime
131  c======================================================================    integer nisurf
132  c======================================================================    logical, intent(in):: debut
133  c Rajout pour l'interface    logical, intent(in):: lafin
134        integer, intent(in):: itime    real zlev1(klon)
135        integer nisurf    real fder(klon), taux(klon), tauy(klon)
136        logical, intent(in):: debut    real temp_air(klon), spechum(klon)
137        logical, intent(in):: lafin    real epot_air(klon), ccanopy(klon)
138        real zlev1(klon)    real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
139        real fder(klon), taux(klon), tauy(klon)    real petBcoef(klon), peqBcoef(klon)
140        real temp_air(klon), spechum(klon)    real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
141        real epot_air(klon), ccanopy(klon)    real p1lay(klon)
142        real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)    !$$$C PB ajout pour soil
143        real petBcoef(klon), peqBcoef(klon)    LOGICAL, intent(in):: soil_model
144        real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)    REAL tsoil(klon, nsoilmx)
145        real p1lay(klon)    REAL qsol(klon)
146  c$$$C PB ajout pour soil  
147        LOGICAL, intent(in):: soil_model    ! Parametres de sortie
148        REAL tsoil(klon, nsoilmx)    real fluxsens(klon), fluxlat(klon)
149        REAL qsol(klon)    real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
150      real emis_new(klon), z0_new(klon)
151  ! Parametres de sortie    real pctsrf_new(klon,nbsrf)
152        real fluxsens(klon), fluxlat(klon)    ! JLD
153        real tsol_rad(klon), tsurf_new(klon), alb_new(klon)    real zzpk
154        real emis_new(klon), z0_new(klon)    !
155        real pctsrf_new(klon,nbsrf)    character (len = 20) :: modname = 'Debut clqh'
156  c JLD    LOGICAL check
157        real zzpk    PARAMETER (check=.false.)
158  C    !
159        character (len = 20) :: modname = 'Debut clqh'    if (check) THEN
160        LOGICAL check       write(*,*) modname,' nisurf=',nisurf
161        PARAMETER (check=.false.)       !C        call flush(6)
162  C    endif
163        if (check) THEN    !
164            write(*,*) modname,' nisurf=',nisurf    if (check) THEN
165  CC        call flush(6)       WRITE(*,*)' qsurf (min, max)' &
166        endif            , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
167  c       !C     call flush(6)
168        if (check) THEN    ENDIF
169         WRITE(*,*)' qsurf (min, max)'    !
170       $     , minval(qsurf(1:knon)), maxval(qsurf(1:knon))    !
171  CC     call flush(6)    if (iflag_pbl.eq.1) then
172        ENDIF       do k = 3, klev
 C  
 C  
       if (iflag_pbl.eq.1) then  
         do k = 3, klev  
           do i = 1, knon  
             gamq(i,k)= 0.0  
             gamt(i,k)=  -1.0e-03  
           enddo  
         enddo  
173          do i = 1, knon          do i = 1, knon
174            gamq(i,2) = 0.0             gamq(i,k)= 0.0
175            gamt(i,2) = -2.5e-03             gamt(i,k)=  -1.0e-03
176          enddo          enddo
177        else       enddo
178          do k = 2, klev       do i = 1, knon
179            do i = 1, knon          gamq(i,2) = 0.0
180              gamq(i,k) = 0.0          gamt(i,2) = -2.5e-03
181              gamt(i,k) = 0.0       enddo
182            enddo    else
183         do k = 2, klev
184            do i = 1, knon
185               gamq(i,k) = 0.0
186               gamt(i,k) = 0.0
187          enddo          enddo
188        endif       enddo
189      endif
190    
191        DO i = 1, knon    DO i = 1, knon
192           psref(i) = paprs(i,1) !pression de reference est celle au sol       psref(i) = paprs(i,1) !pression de reference est celle au sol
193           local_ts(i) = ts(i)       local_ts(i) = ts(i)
194        ENDDO    ENDDO
195        DO k = 1, klev    DO k = 1, klev
196        DO i = 1, knon       DO i = 1, knon
197           zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA          zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
198           zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA          zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
199           local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)          local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
200           local_q(i,k) = q(i,k)          local_q(i,k) = q(i,k)
201        ENDDO       ENDDO
202        ENDDO    ENDDO
203  c    !
204  c Convertir les coefficients en variables convenables au calcul:    ! Convertir les coefficients en variables convenables au calcul:
205  c    !
206  c    !
207        DO k = 2, klev    DO k = 2, klev
208        DO i = 1, knon       DO i = 1, knon
209           zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))          zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
210       .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2               *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
211           zx_coef(i,k) = zx_coef(i,k) * dtime*RG          zx_coef(i,k) = zx_coef(i,k) * dtime*RG
212        ENDDO       ENDDO
213        ENDDO    ENDDO
214  c    !
215  c Preparer les flux lies aux contre-gardients    ! Preparer les flux lies aux contre-gardients
216  c    !
217        DO k = 2, klev    DO k = 2, klev
218        DO i = 1, knon       DO i = 1, knon
219           zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)          zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k) &
220       .              *(pplay(i,k-1)-pplay(i,k))               *(pplay(i,k-1)-pplay(i,k))
221           z_gamaq(i,k) = gamq(i,k) * zdelz          z_gamaq(i,k) = gamq(i,k) * zdelz
222           z_gamah(i,k) = gamt(i,k) * zdelz *RCPD * zx_pkh(i,k)          z_gamah(i,k) = gamt(i,k) * zdelz *RCPD * zx_pkh(i,k)
223        ENDDO       ENDDO
224        ENDDO    ENDDO
225        DO i = 1, knon    DO i = 1, knon
226           zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)       zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
227           zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)       zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev) &
228       .                   -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)            -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
229           zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)       zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
230  c       !
231           zzpk=(pplay(i,klev)/psref(i))**RKAPPA       zzpk=(pplay(i,klev)/psref(i))**RKAPPA
232           zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)       zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
233           zx_ch(i,klev) = (local_h(i,klev)*zzpk*delp(i,klev)       zx_ch(i,klev) = (local_h(i,klev)*zzpk*delp(i,klev) &
234       .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)            -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
235           zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)       zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
236        ENDDO    ENDDO
237        DO k = klev-1, 2 , -1    DO k = klev-1, 2 , -1
238        DO i = 1, knon       DO i = 1, knon
239           zx_buf1(i) = delp(i,k)+zx_coef(i,k)          zx_buf1(i) = delp(i,k)+zx_coef(i,k) &
240       .               +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))               +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
241           zx_cq(i,k) = (local_q(i,k)*delp(i,k)          zx_cq(i,k) = (local_q(i,k)*delp(i,k) &
242       .                 +zx_coef(i,k+1)*zx_cq(i,k+1)               +zx_coef(i,k+1)*zx_cq(i,k+1) &
243       .                 +zx_coef(i,k+1)*z_gamaq(i,k+1)               +zx_coef(i,k+1)*z_gamaq(i,k+1) &
244       .                 -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)               -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
245           zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)          zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
246  c          !
247           zzpk=(pplay(i,k)/psref(i))**RKAPPA          zzpk=(pplay(i,k)/psref(i))**RKAPPA
248           zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k)          zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k) &
249       .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
250           zx_ch(i,k) = (local_h(i,k)*zzpk*delp(i,k)          zx_ch(i,k) = (local_h(i,k)*zzpk*delp(i,k) &
251       .                 +zx_coef(i,k+1)*zx_ch(i,k+1)               +zx_coef(i,k+1)*zx_ch(i,k+1) &
252       .                 +zx_coef(i,k+1)*z_gamah(i,k+1)               +zx_coef(i,k+1)*z_gamah(i,k+1) &
253       .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)               -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
254           zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)          zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
255        ENDDO       ENDDO
256        ENDDO    ENDDO
257  C    !
258  C nouvelle formulation JL Dufresne    ! nouvelle formulation JL Dufresne
259  C    !
260  C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt    ! q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
261  C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt    ! h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
262  C    !
263        DO i = 1, knon    DO i = 1, knon
264           zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))       zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
265           zx_cq(i,1) = (local_q(i,1)*delp(i,1)       zx_cq(i,1) = (local_q(i,1)*delp(i,1) &
266       .                 +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))            +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2))) &
267       .                /zx_buf1(i)            /zx_buf1(i)
268           zx_dq(i,1) = -1. * RG / zx_buf1(i)       zx_dq(i,1) = -1. * RG / zx_buf1(i)
269  c       !
270           zzpk=(pplay(i,1)/psref(i))**RKAPPA       zzpk=(pplay(i,1)/psref(i))**RKAPPA
271           zx_buf2(i) = zzpk*delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))       zx_buf2(i) = zzpk*delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
272           zx_ch(i,1) = (local_h(i,1)*zzpk*delp(i,1)       zx_ch(i,1) = (local_h(i,1)*zzpk*delp(i,1) &
273       .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))            +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2))) &
274       .                /zx_buf2(i)            /zx_buf2(i)
275           zx_dh(i,1) = -1. * RG / zx_buf2(i)       zx_dh(i,1) = -1. * RG / zx_buf2(i)
276        ENDDO    ENDDO
277    
278  C Appel a interfsurf (appel generique) routine d'interface avec la surface    ! Appel a interfsurf (appel generique) routine d'interface avec la surface
279    
280  c initialisation    ! initialisation
281         petAcoef =0.    petAcoef =0.
282          peqAcoef = 0.    peqAcoef = 0.
283          petBcoef =0.    petBcoef =0.
284          peqBcoef = 0.    peqBcoef = 0.
285          p1lay =0.    p1lay =0.
286            
287  c      do i = 1, knon    !      do i = 1, knon
288          petAcoef(1:knon) = zx_ch(1:knon,1)    petAcoef(1:knon) = zx_ch(1:knon,1)
289          peqAcoef(1:knon) = zx_cq(1:knon,1)    peqAcoef(1:knon) = zx_cq(1:knon,1)
290          petBcoef(1:knon) =  zx_dh(1:knon,1)    petBcoef(1:knon) =  zx_dh(1:knon,1)
291          peqBcoef(1:knon) = zx_dq(1:knon,1)    peqBcoef(1:knon) = zx_dq(1:knon,1)
292          tq_cdrag(1:knon) =coef(1:knon,1)    tq_cdrag(1:knon) =coef(1:knon,1)
293          temp_air(1:knon) =t(1:knon,1)    temp_air(1:knon) =t(1:knon,1)
294          epot_air(1:knon) =local_h(1:knon,1)    epot_air(1:knon) =local_h(1:knon,1)
295          spechum(1:knon)=q(1:knon,1)    spechum(1:knon)=q(1:knon,1)
296          p1lay(1:knon) = pplay(1:knon,1)    p1lay(1:knon) = pplay(1:knon,1)
297          zlev1(1:knon) = delp(1:knon,1)    zlev1(1:knon) = delp(1:knon,1)
298  c        swnet = swdown * (1. - albedo)    !        swnet = swdown * (1. - albedo)
299  c    !
300  cIM swdown=flux SW incident sur terres    !IM swdown=flux SW incident sur terres
301  cIM swdown=flux SW net sur les autres surfaces    !IM swdown=flux SW net sur les autres surfaces
302  cIM     swdown(1:knon) = swnet(1:knon)    !IM     swdown(1:knon) = swnet(1:knon)
303          if(nisurf.eq.is_ter) THEN    if(nisurf.eq.is_ter) THEN
304           swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon))       swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon))
305          else    else
306           swdown(1:knon) = swnet(1:knon)       swdown(1:knon) = swnet(1:knon)
307          endif    endif
308  c      enddo    !      enddo
309        ccanopy = co2_ppm    ccanopy = co2_ppm
310    
311        CALL interfsurf_hq(itime, dtime, date0, jour, rmu0,    CALL interfsurf_hq(itime, dtime, date0, jour, rmu0, &
312       e klon, iim, jjm, nisurf, knon, knindex, pctsrf,         klon, iim, jjm, nisurf, knon, knindex, pctsrf,  &
313       e rlon, rlat, cufi, cvfi,         rlon, rlat, cufi, cvfi,  &
314       e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol,         debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol, &
315       e zlev1,  u1lay, v1lay, temp_air, spechum, epot_air, ccanopy,         zlev1,  u1lay, v1lay, temp_air, spechum, epot_air, ccanopy,  &
316       e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,         tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
317       e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,         precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
318       e fder, taux, tauy,         fder, taux, tauy, &
319  c -- LOOP         ywindsp, rugos, rugoro, &
320       e ywindsp,         albedo, snow, qsurf, &
321  c -- LOOP         ts, p1lay, psref, radsol, &
322       e rugos, rugoro,         ocean, npas, nexca, zmasq, &
323       e albedo, snow, qsurf,         evap, fluxsens, fluxlat, dflux_l, dflux_s,               &
324       e ts, p1lay, psref, radsol,         tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new,  &
325       e ocean, npas, nexca, zmasq,         pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0, &
326       s evap, fluxsens, fluxlat, dflux_l, dflux_s,                       flux_o, flux_g, tslab, seaice)
327       s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new,  
328       s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,  
329  cIM "slab" ocean    do i = 1, knon
330       s flux_o, flux_g, tslab, seaice)       flux_t(i,1) = fluxsens(i)
331         flux_q(i,1) = - evap(i)
332         d_ts(i) = tsurf_new(i) - ts(i)
333        do i = 1, knon       albedo(i) = alb_new(i)
334          flux_t(i,1) = fluxsens(i)    enddo
335          flux_q(i,1) = - evap(i)  
336          d_ts(i) = tsurf_new(i) - ts(i)    !==== une fois on a zx_h_ts, on peut faire l'iteration ========
337          albedo(i) = alb_new(i)    DO i = 1, knon
338        enddo       local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
339         local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i,1)*dtime
340  c==== une fois on a zx_h_ts, on peut faire l'iteration ========    ENDDO
341        DO i = 1, knon    DO k = 2, klev
342           local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime       DO i = 1, knon
          local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i,1)*dtime  
       ENDDO  
       DO k = 2, klev  
       DO i = 1, knon  
343          local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)          local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
344          local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)          local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
345        ENDDO       ENDDO
346        ENDDO    ENDDO
347  c======================================================================    !======================================================================
348  c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas    !== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
349  c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)    !== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
350        DO k = 2, klev    DO k = 2, klev
351        DO i = 1, knon       DO i = 1, knon
352          flux_q(i,k) = (zx_coef(i,k)/RG/dtime)          flux_q(i,k) = (zx_coef(i,k)/RG/dtime) &
353       .                * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))               * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
354          flux_t(i,k) = (zx_coef(i,k)/RG/dtime)          flux_t(i,k) = (zx_coef(i,k)/RG/dtime) &
355       .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))               * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k)) &
356       .                / zx_pkh(i,k)               / zx_pkh(i,k)
357        ENDDO       ENDDO
358        ENDDO    ENDDO
359  c======================================================================    !======================================================================
360  C Calcul tendances    ! Calcul tendances
361        DO k = 1, klev    DO k = 1, klev
362        DO i = 1, knon       DO i = 1, knon
363           d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)          d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
364           d_q(i,k) = local_q(i,k) - q(i,k)          d_q(i,k) = local_q(i,k) - q(i,k)
365        ENDDO       ENDDO
366        ENDDO    ENDDO
 c  
367    
368        RETURN  END SUBROUTINE clqh
       END  
Legend:
Removed from v.12  
changed lines
  Added in v.38
  ViewVC Help
Powered by ViewVC 1.1.21