Changeset 86 for LMDZ.3.3

Timestamp:

May 10, 2000, 12:30:30 PM (24 years ago)

Author:

lmdzadmin

Message:

nouvelle version de la routine incluant l'appel a la routine d'interface
avec la surface
LF

File:

• LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F (modified) (7 diffs)

LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F

@@ -143 +143 @@
       REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
 c======================================================================
+
+      write(*,*)'CLMAIN.NEW'
+
       DO k = 1, klev   ! epaisseur de couche
       DO i = 1, klon
@@ -322 +325 @@
 c
 c calculer la diffusion de "q" et de "h"
-      CALL clqh(knon, dtime, yu1, yv1,
+      CALL clqh(knon, dtime, nsrf,yu1, yv1,
      e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
      e          ycal,ybeta,ydif,ygamt,ygamq,
@@ -429 +432 @@
       RETURN
       END
-      SUBROUTINE clqh(knon,dtime,u1lay,v1lay,coef,t,q,ts,paprs,pplay,
+      SUBROUTINE clqh(knon,dtime,nisurf,u1lay,v1lay,coef,
+     e                t,q,ts,paprs,pplay,
      e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
      s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
+
+      USE interface_surf
+
       IMPLICIT none
 c======================================================================
@@ -517 +524 @@
 #include "YOETHF.h"
 #include "FCTTRE.h"
-c======================================================================
+#include "indicesol.h"
+c======================================================================
+c Rajout pour l'interface
+      integer itime
+      integer jour
+      integer nisurf
+      integer knindex(klon)
+      logical debut, lafin, ok_veget
+      real rlon(klon), rlat(klon)
+      real zlev(klon), zlflu(klon)
+      real temp_air(klon), spechum(klon)
+      real hum_air(klon), ccanopy(klon)
+      real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
+      real petBcoef(klon), peqBcoef(klon)
+      real precip_rain(klon), precip_snow(klon)
+      real lwdown(klon), swnet(klon), swdown(klon), ps(klon)
+      real p1lay(klon)
+      real coef1lay(klon)
+      character*6 ocean
+
+! Parametres de sortie
+      real evap(klon), fluxsens(klon), fluxlat(klon)
+      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
+      real emis_new(klon), z0_new(klon)
+      real dflux_l(klon), dflux_s(klon)
+      real pctsrf_new(klon,nbsrf)
+c
+
       DO i = 1, knon
          psref(i) = paprs(i,1) !pression de reference est celle au sol
@@ -533 +567 @@
 c Convertir les coefficients en variables convenables au calcul:
 c
-      DO i = 1, knon
-         zx_coef(i,1) = coef(i,1)
-     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
-     .                 * pplay(i,1)/(RD*t(i,1))
-      ENDDO
 c
       DO k = 2, klev
@@ -557 +586 @@
       ENDDO
       ENDDO
-c======================================================================
-c
-c zx_qs = qsat en kg/kg
-c
-      DO i = 1, knon
-         IF (thermcep) THEN
-           zdelta=MAX(0.,SIGN(1.,rtt-local_ts(i)))
-           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
-           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
-           zx_qs= r2es * FOEEW(local_ts(i),zdelta)/paprs(i,1)
-           zx_qs=MIN(0.5,zx_qs)
-           zcor=1./(1.-retv*zx_qs)
-           zx_qs=zx_qs*zcor
-           zx_dq_s_dh = FOEDE(local_ts(i),zdelta,zcvm5,zx_qs,zcor)
-     .                 /RLVTT / zx_pkh(i,1)
-         ELSE
-           IF (local_ts(i).LT.t_coup) THEN
-              zx_qs = qsats(local_ts(i)) / paprs(i,1)
-              zx_dq_s_dh = dqsats(local_ts(i),zx_qs)/RLVTT
-     .                    / zx_pkh(i,1)
-           ELSE
-              zx_qs = qsatl(local_ts(i)) / paprs(i,1)
-              zx_dq_s_dh = dqsatl(local_ts(i),zx_qs)/RLVTT
-     .                    / zx_pkh(i,1)
-           ENDIF
-         ENDIF
 
-c        zx_dq_s_dh = 0.
-         zx_dq_s_dt(i) = RCPD * zx_pkh(i,1) * zx_dq_s_dh
-         zx_qsat(i) = zx_qs
-c
-      ENDDO
       DO i = 1, knon
          zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
@@ -639 +637 @@
       ENDDO
 
-C === Calcul de la temperature de surface ===
-C 
-C zx_sl = chaleur latente d'evaporation ou de sublimation
-C
-      do i = 1, knon
-        zx_sl(i) = RLVTT
-        if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT
-        zx_k1(i) = zx_coef(i,1)
-      enddo
-      do i = 1, knon
-C Q
-        zx_oq(i) = 1. - (beta(i) * zx_k1(i) * zx_dq(i,1) * dtime)
-        zx_mq(i) = beta(i) * zx_k1(i) * 
-     .              (zx_cq(i,1) - zx_qsat(i) 
-     .                             + zx_dq_s_dt(i) * local_ts(i))
-     .              / zx_oq(i)
-        zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i))
-     .                                                / zx_oq(i)
-c H
-        zx_oh(i) = 1. - (zx_k1(i) * zx_dh(i,1) * dtime)
-        zx_mh(i) = zx_k1(i) * zx_ch(i,1) / zx_oh(i)
-        zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i,1))/ zx_oh(i)
-c Tsurface
-        local_ts(i) = (ts(i) + cal(i)/(RCPD * zx_pkh(i,1)) * dtime *
-     .                   (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) 
-     .                       + dif_grnd(i) * t_grnd * dtime)/
-     .                ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i,1)) * (
-     .                       zx_nh(i) + zx_sl(i) * zx_nq(i))  
-     .                     + dtime * dif_grnd(i))
-        zx_h_ts(i) = local_ts(i) * RCPD * zx_pkh(i,1)
-        d_ts(i) = local_ts(i) - ts(i)
-        zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
-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)
-        flux_q(i,1) = zx_mq(i) + zx_nq(i) * local_ts(i) 
-        flux_t(i,1) = zx_mh(i) + zx_nh(i) * local_ts(i)
-c Derives des flux dF/dTs (W m-2 K-1):
-        dflux_s(i) = zx_nh(i)
-        dflux_l(i) = (zx_sl(i) * zx_nq(i))
-      ENDDO
+C Appel a interfsurf (appel generique) routine d'interface avec la surface
 
+      ok_veget = .false.
+      ocean = 'force '
+
+      petAcoef=zx_ch(:,1)
+      peqAcoef=zx_cq(:,1)
+      petBcoef=zx_dh(:,1)
+      peqBcoef=zx_dq(:,1)
+      coef1lay=coef(:,1)
+      temp_air=t(:,1)
+      spechum=q(:,1)
+      p1lay = pplay(:,1)
+
+      CALL interfsurf(itime, dtime, jour,
+     . klon, nisurf, knon, knindex, rlon, rlat,
+     . debut, lafin, ok_veget,
+     . zlev, zlflu, u1lay, v1lay, temp_air, spechum, hum_air, ccanopy, 
+     . tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
+     . precip_rain, precip_snow, lwdown, swnet, swdown,
+     . ts, p1lay, cal, beta, coef1lay, psref, radsol, dif_grnd,
+     . ocean,
+     . evap, fluxsens, fluxlat, dflux_l, dflux_s,              
+     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
+
+      flux_t(:,1) = fluxsens
+      flux_q(:,1) = fluxlat
+      d_ts = tsurf_new - ts
 
 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)*dtime
-         local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i)*dtime
+         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