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

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

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

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