Changeset 407 for LMDZ.3.3/branches

Timestamp:

Oct 15, 2002, 4:59:04 PM (22 years ago)

Author:

lmdzadmin

Message:

Diagnostics supplementaires: T2m, q2m, u et v a 10m JP/IM
IM/LF

File:

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

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

@@ -16 +16 @@
      .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
      .                  dflux_t,dflux_q,
-     .                  zcoefh,zu1,zv1)
+     .                  zcoefh,zu1,zv1, t2m, q2m, u10m, v10m)
 cAA .                  itr, tr, flux_surf, d_tr)
 cAA REM:
@@ -35 +35 @@
       USE ioipsl
       USE interface_surf
+      USE stdlevvar_int
       IMPLICIT none
 c======================================================================
@@ -169 +170 @@
 c
 #include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
       REAL u1lay(klon), v1lay(klon)
       REAL delp(klon,klev)
@@ -203 +206 @@
       integer idayref
 #include "temps.h"
-      
+      REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
+      REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
+c
+      REAL yt2m(klon), yq2m(klon), yu10m(klon)
+c
+      REAL uzon(klon), vmer(klon)
+      REAL tair1(klon), qair1(klon), tairsol(klon)
+      REAL psfce(klon), patm(klon)
+c
+      REAL qairsol(klon), zgeo1(klon)
+      REAL rugo1(klon)
+c
+      LOGICAL zxli ! utiliser un jeu de fonctions simples
+      PARAMETER (zxli=.FALSE.)
+c
+      REAL zt, zqs, zdelta, zcor
+      REAL t_coup
+      PARAMETER(t_coup=273.15)
+C
+      PRINT*,'IMclmain klon=',klon
       IF (first_appel) THEN
           first_appel=.false.
@@ -276 +298 @@
       yv1 = 0.0
       yrads = 0.0
-      ypaprs = 0.0
       ypaprs = 0.0
       ypplay = 0.0
@@ -573 +594 @@
       ENDDO
 c
+c
+#undef T2m     
+#ifdef T2m
+ccc diagnostic t,q a 2m et u, v a 10m
+c
+      DO j=1, knon
+        i = ni(j)
+        uzon(j) = yu(j,1) + y_d_u(j,1)
+        vmer(j) = yv(j,1) + y_d_v(j,1)
+        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.EQ.is_oce) THEN
+         rugo1(j) = rugos(i,nsrf)
+        ENDIF
+        psfce(j)=ypaprs(j,1)
+        patm(j)=ypplay(j,1)
+c
+        IF (nsrf.EQ.1) THEN
+          qairsol(j) = yqsol(j)
+        ELSE IF(nsrf.GT.1) THEN
+         zt = ts(i,nsrf)
+         IF (thermcep) THEN
+           zdelta = MAX(0.,SIGN(1.,RTT-zt))
+           zqs = R2ES * FOEEW(zt,zdelta) / ypplay(j,1)
+           zqs = MIN(0.5,zqs)
+           zcor = 1./(1.-RETV*zqs)
+           zqs = zqs*zcor
+         ELSE
+           IF (zt .LT. t_coup) THEN
+             zqs = qsats(zt) / ypplay(j,1)
+           ELSE
+             zqs = qsatl(zt) / ypplay(j,1)
+           ENDIF
+         ENDIF   
+         qairsol(j) = zqs
+        ENDIF   
+      ENDDO
+c
+      IF(nsrf.EQ.3) THEN
+       j=1465
+       WRITE(*,*)' INstO',klon,knon,nsrf,zxli,uzon(j),vmer(j), 
+     & tair1(j),qair1(j),zgeo1(j),tairsol(j),qairsol(j),rugo1(j),
+     & psfce(j),patm(j)
+      ENDIF
+c
+      CALL stdlevvar(klon, knon, nsrf, zxli,
+     &               uzon, vmer, tair1, qair1, zgeo1,
+     &               tairsol, qairsol, rugo1, psfce, patm,
+     &               yt2m, yq2m, yu10m)
+
+c
+      IF(nsrf.EQ.3) THEN
+       j=1465
+       WRITE(*,*)' OUstd',klon,knon,nsrf,zxli,uzon(j),vmer(j), 
+     & tair1(j),qair1(j),zgeo1(j),tairsol(j),qairsol(j),rugo1(j),
+     & psfce(j),patm(j)
+       WRITE(*,*)' tqu',yt2m(j),yq2m(j),yu10m(j)
+      ENDIF
+c
+      DO j=1, knon
+       i = ni(j)
+       t2m(i,nsrf)=yt2m(j)
+
+        IF(nsrf.EQ.3) THEN
+         IF(j.EQ.1465) THEN
+         WRITE(*,*) 't2m APRES stdlev',j,i,tair1(j),t2m(i,nsrf),
+     $   tairsol(j),rlon(i),rlat(i)
+         ENDIF
+        ENDIF
+c
+       q2m(i,nsrf)=yq2m(j)
+c
+c 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)
+c
+      ENDDO
+#else
+       DO j=1, knon
+         t2m(i,nsrf)=0.
+         q2m(i,nsrf)=0.
+         u10m(i,nsrf)=0.
+         v10m(i,nsrf)=0.
+       ENDDO 
+#endif
 99999 CONTINUE
 c
@@ -1028 +1138 @@
      .                  u,v,t,q,
      .                  pcfm, pcfh)
+      USE clcdrag_int
       IMPLICIT none
 c======================================================================
@@ -1070 +1181 @@
       REAL cepdu2, ckap, cb, cc, cd, clam
       PARAMETER (cepdu2 =(0.1)**2)
-      PARAMETER (ckap=0.35)
+      PARAMETER (CKAP=0.4)
       PARAMETER (cb=5.0)
       PARAMETER (cc=5.0)
@@ -1103 +1214 @@
       REAL zri(klon)
       REAL zl2(klon)
-      REAL zcfm1(klon), zcfm2(klon)
-      REAL zcfh1(klon), zcfh2(klon)
-      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
-      REAL zscf, zucf, zcr
+
+      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
+      REAL pcfm1(klon), pcfh1(klon)
+c
+      REAL zdphi, zdu2, ztvd, ztvu, zcdn
+      REAL zscf
       REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
       REAL z2geomf, zalh2, zalm2, zscfh, zscfm
@@ -1116 +1229 @@
 c essai qsurf
       real qsurf(klon)
-      real friv, frih
 c
       LOGICAL appel1er
@@ -1188 +1300 @@
 c Calculer le frottement au sol (Cdrag)
 c
-      CALL clcdrag(knon, nsrf, zxli, u, v, t, q, zgeop,
-     .             ts, qsurf, rugos, pcfm, pcfh, zcdn, zri)
+c     DO i = 1, knon
+      DO i = 1, klon
+       u1(i) = u(i,1)
+       v1(i) = v(i,1)
+       t1(i) = t(i,1)
+       q1(i) = q(i,1)
+       z1(i) = zgeop(i,1)
+      ENDDO
+c
+      CALL clcdrag(klon, knon, nsrf, zxli, 
+     $             u1, v1, t1, q1, z1,
+     $             ts, qsurf, rugos,
+     $             pcfm1, pcfh1) 
+C
+      DO i = 1, knon
+       pcfm(i,1)=pcfm1(i)
+       pcfh(i,1)=pcfh1(i)
+      ENDDO
 c
 c Calculer les coefficients turbulents dans l'atmosphere
@@ -1196 +1324 @@
          itop(i) = isommet
       ENDDO
+
+      PRINT*,' isommet=',isommet,' knon=',knon
 
       DO k = 2, isommet
@@ -1310 +1440 @@
       END
 
-      SUBROUTINE clcdrag(knon, nsrf, zxli, 
-     .                   u, v, t, q, zgeop, 
-     .                   ts, qsurf, rugos,
-     .                   pcfm, pcfh, zcdn, zri)
-c ================================================================= c
-c Objet : calcul cdrags pour le moment et les flux chaleur sensible, latente (pcfm,pcfh)
-c         et du nombre de Richardson zri
-c ================================================================= c
-      IMPLICIT NONE
-#include "dimensions.h"
-#include "dimphy.h"
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "indicesol.h"
-c
-      INTEGER knon, nsrf
-      REAL ts(klon), qsurf(klon)
-      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
-      REAL zgeop(klon,klev)
-      REAL rugos(klon), zri(klon)
-c
-      REAL zcdn
-      REAL pcfm(klon,klev), pcfh(klon,klev)
-c
-c Quelques constantes et options:
-c
-      REAL ckap, cb, cc, cd, cepdu2
-      PARAMETER (ckap=0.35)
-      PARAMETER (cb=5.0)
-      PARAMETER (cc=5.0)
-      PARAMETER (cd=5.0)
-      PARAMETER (cepdu2 =(0.1)**2)
-c
-c Variables locales
-      INTEGER i
-      REAL zdu2, zdphi, ztsolv, ztvd, zscf, zucf, zcr
-      REAL friv, frih
-      REAL zcfm1(klon), zcfm2(klon)
-      REAL zcfh1(klon), zcfh2(klon)
-c
-c Fonctions thermodynamiques et fonctions d'instabilite
-      REAL fsta, fins, x
-      LOGICAL zxli
-      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
-      fins(x) = SQRT(1.0-18.0*x)
-c
-c Calculer le frottement au sol (Cdrag)
-c
-      DO i = 1, knon
-         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
-         zdphi=zgeop(i,1)
-c        ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
-         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
-         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
-     .       *(1.+RETV*q(i,1))
-         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
-         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
-         IF (zri(i) .ge. 0.) THEN ! situation stable
-           IF (.NOT.zxli) THEN
-           zscf=SQRT(1.+cd*ABS(zri(i)))
-           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
-!           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf)
-           zcfm1(i) = zcdn * FRIV
-           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
-!           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
-           zcfh1(i) = zcdn * FRIH
-           pcfm(i,1) = zcfm1(i)
-           pcfh(i,1) = zcfh1(i)
-           ELSE
-           pcfm(i,1) = zcdn* fsta(zri(i))
-           pcfh(i,1) = zcdn* fsta(zri(i))
-           ENDIF
-         ELSE ! situation instable
-           IF (.NOT.zxli) THEN
-           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
-     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
-           zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
-           zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
-           pcfm(i,1) = zcfm2(i)
-           pcfh(i,1) = zcfh2(i)
-           ELSE
-           pcfm(i,1) = zcdn* fins(zri(i))
-           pcfh(i,1) = zcdn* fins(zri(i))
-           ENDIF
-           zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
-         ENDIF
-      END DO
-      RETURN
-      END
-
       SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
      .                  pcfm, pcfh)
@@ -1611 +1650 @@
       PARAMETER (isommet=klev)
       REAL vk
-      PARAMETER (vk=0.35)
+      PARAMETER (vk=0.40)
       REAL ricr
       PARAMETER (ricr=0.4)