Ignore:
Timestamp:
Aug 6, 2003, 4:50:49 PM (21 years ago)
Author:
lmdzadmin
Message:

Modifs sur les seuils (cdrag etc...), inclusion des diagnostics ISCCP par Ionela
LF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F

    r463 r467  
    191191      PARAMETER(klevp1=klev+1)
    192192#include "raddim.h"
    193       REAL*8 ZFSUP(KDLON,KFLEV+1)
    194       REAL*8 ZFSDN(KDLON,KFLEV+1)
    195       REAL*8 ZFSUP0(KDLON,KFLEV+1)
    196       REAL*8 ZFSDN0(KDLON,KFLEV+1)
     193cc      REAL*8 ZFSUP(KDLON,KFLEV+1)
     194cc      REAL*8 ZFSDN(KDLON,KFLEV+1)
     195cc      REAL*8 ZFSUP0(KDLON,KFLEV+1)
     196cc      REAL*8 ZFSDN0(KDLON,KFLEV+1)
     197c
     198      REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
     199      SAVE swdn0 , swdn, swup0, swup
    197200
    198201cccIM cf. FH
    199202      real u850(klon),v850(klon),u200(klon),v200(klon)
    200203      real u500(klon),v500(klon),phi500(klon),w500(klon)
     204cIM
     205      real prw(klon)
     206
     207cIM ISCCP - proprietes microphysiques des nuages convectifs
     208      REAL convliq(klon,klev)  ! eau liquide nuageuse convective
     209      REAL convfra(klon,klev)  ! fraction nuageuse convective
     210
     211      REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
     212      REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
     213      REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
     214      REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
     215
     216      INTEGER kinv, linv
     217
     218cIM ISCCP simulator BEGIN
     219      INTEGER igfi2D(iim,jjmp1)
     220cv3.4
     221      INTEGER debug, debugcol
     222      INTEGER npoints
     223      PARAMETER(npoints=klon)
     224      INTEGER sunlit(klon)
     225
     226      INTEGER ncol, seed(klon)
     227
     228cIM dans clesphys.h top_height, overlap
     229c     PARAMETER(ncol=100)
     230c     PARAMETER(ncol=625)
     231      PARAMETER(ncol=10)
     232      REAL tautab(0:255)
     233      INTEGER invtau(-20:45000)
     234      REAL emsfc_lw
     235      PARAMETER(emsfc_lw=0.99)
     236      REAL    ran0                      ! type for random number fuction
     237
     238      REAL pfull(klon,klev)
     239      REAL phalf(klon,klev+1)
     240      REAL cldtot(klon,klev)
     241      REAL dtau_s(klon,klev)
     242      REAL dtau_c(klon,klev)
     243      REAL dem_s(klon,klev)
     244      REAL dem_c(klon,klev)
     245cPLUS : variables de haut en bas pour le simulateur ISCCP
     246      REAL qv(klon,klev)
     247      REAL cc(klon,klev)
     248      REAL conv(klon,klev)
     249      REAL dtau_sH2B(klon,klev)
     250      REAL dtau_cH2B(klon,klev)
     251      REAL at(klon,klev)
     252      REAL dem_sH2B(klon,klev)
     253      REAL dem_cH2B(klon,klev)
     254
     255c output from ISCCP
     256      REAL fq_isccp(klon,7,7)
     257      REAL totalcldarea(klon)
     258      REAL meanptop(klon)
     259      REAL meantaucld(klon)
     260      REAL boxtau(klon,ncol)
     261      REAL boxptop(klon,ncol)
     262
     263c grille 4d physique
     264      INTEGER l, ni, nj, kmax, lmax, nrec
     265      INTEGER ni1, ni2, nj1, nj2
     266c     PARAMETER(kmax=7, lmax=7)
     267      PARAMETER(kmax=8, lmax=8)
     268      INTEGER kmaxm1, lmaxm1
     269      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
     270c     INTEGER iimx7, jjmx7, jjmp1x7
     271c     PARAMETER(iimx7=iim*7, jjmx7=jjm*7, jjmp1x7=jjmp1*7)
     272c     REAL fq4d(iim,jjmp1,7,7)
     273c     REAL fq3d(iimx7, jjmp1x7)
     274      INTEGER iimx8, jjmx8, jjmp1x8
     275      PARAMETER(iimx8=iim*8, jjmx8=jjm*8, jjmp1x8=jjmp1*8)
     276      REAL fq4d(iim,jjmp1,8,8)
     277      REAL fq3d(iimx8, jjmp1x8)
     278cIM180603     SAVE fq3d
     279
     280c     REAL maxfq3d, minfq3d
     281c
     282      INTEGER iw, iwmax
     283      REAL wmin, pas_w
     284c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
     285      PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
     286      REAL o500(klon)
     287      INTEGER nreg, nbreg
     288      PARAMETER(nbreg=5)
     289c     REAL histoW(iwmax,kmaxm1,lmaxm1)
     290      REAL histoW(kmaxm1,lmaxm1,iwmax,nbreg)
     291      REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbreg)
     292cIM180603     
     293c     SAVE histoW, nhistoW
     294c     SAVE nhistoW
     295      REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbreg)
     296      SAVE nhistoWt
     297
     298c     REAL histoWinv(kmaxm1,lmaxm1,iwmax)
     299c     REAL nhistoW(kmaxm1,lmaxm1,iwmax)
     300      INTEGER linv
     301c     LOGICAL pct_ocean(klon,nbreg)
     302      INTEGER pct_ocean(klon,nbreg)
     303      REAL rlonPOS(klon)
     304c     CHARACTER*4 pdirect
     305 
     306c sorties ISCCP
     307
     308      logical ok_isccp
     309      real ecrit_isccp
     310      integer nid_isccp
     311      save ok_isccp, ecrit_isccp, nid_isccp       
     312
     313#define histISCCP
     314#undef histISCCP
     315#ifdef histISCCP
     316c     data ok_isccp,ecrit_isccp/.true.,0.125/     
     317c     data ok_isccp,ecrit_isccp/.true.,1./     
     318      data ok_isccp/.true./     
     319#else
     320      data ok_isccp/.false./
     321#endif
     322
     323      REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
     324c     DATA zx_tau/0.1, 1.3, 3.6, 9.4, 23., 60./
     325c     DATA zx_pc/50., 180., 310., 440., 560., 680., 800., 1015./
     326c     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
     327cOK     DATA zx_tau/0.0, 0.1, 1.3, 3.6, 9.4, 23., 60./
     328cOK     DATA zx_pc/800., 680., 560., 440., 310., 180., 50./
     329
     330c tester l'alure
     331      DATA zx_tau/1., 2., 3., 4., 5., 6., 7./
     332c     DATA zx_pc/1., 2., 3., 4., 5., 6., 7./
     333      DATA zx_pc/7., 6., 5., 4., 3., 2., 1./
     334
     335      INTEGER komega, nhoriRD
     336
     337c statistiques regime dynamique END
     338
     339c     REAL del_lon(iim), del_lat(jjmp1)
     340      REAL del_lon, del_lat
     341c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
     342      REAL zx_lonx8(iimx8), zx_latx8(jjmp1x8)
     343c     INTEGER nhorix7
     344      INTEGER nhorix8
     345
     346cIM ISCCP simulator END
    201347
    202348      logical ok_hf
    … …  
    497643      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
    498644cccIM
    499       SAVE  ZFSUP,ZFSDN,ZFSUP0,ZFSDN0
    500645
    501646      INTEGER itaprad
    … …  
    753898         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
    754899     .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
    755      .       falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
     900     .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
    756901     .       dlw,radsol,frugs,agesno,clesphy0,
    757902     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
    … …  
    8801025c   Initialisation des sorties
    8811026c=============================================================
     1027
     1028#ifdef histISCCP
     1029#include "ini_histISCCP.h"
     1030#endif
     1031
    8821032#ifdef histhf
    8831033#include "ini_histhf.h"
    … …  
    9581108      ENDIF
    9591109C
    960       IF (if_ebil.ge.1) THEN
    9611110        DO i = 1, klon
    9621111          ztsol(i) = 0.
    … …  
    9671116          ENDDO
    9681117        ENDDO
     1118C
     1119      IF (if_ebil.ge.1) THEN
    9691120        ztit='after dynamic'
    9701121        CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
    … …  
    10721223      DO nsrf = 1, nbsrf
    10731224      DO i = 1, klon
    1074          frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
    1075 cccc        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
     1225c         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
     1226        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
    10761227      ENDDO
    10771228      ENDDO
    … …  
    10911242        rmu0 = -999.999
    10921243      ENDIF
    1093 C
     1244cIM BEG
     1245      DO i=1, klon
     1246       sunlit(i)=1
     1247       IF(rmu0(i).EQ.0.) sunlit(i)=0
     1248c      IF(rmu0(i).EQ.0.) THEN
     1249c       sunlit(i)=0
     1250c       PRINT*,' il fait nuit ',i,rlat(i),rlon(i)
     1251c      ENDIF
     1252      ENDDO
     1253cIM END
    10941254C     Calcul de l'abedo moyen par maille
    10951255      albsol(:)=0.
    … …  
    11031263C
    11041264C     Repartition sous maille des flux LW et SW
    1105       DO nsrf = 1, nbsrf
    1106       DO i = 1, klon
    1107         fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
    1108         fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
    1109       ENDDO
    1110       ENDDO
     1265C Modif OM+PASB+JLD
     1266C Repartition du longwave par sous-surface linearisee
     1267Cn
     1268       DO nsrf = 1, nbsrf
     1269       DO i = 1, klon
     1270c$$$        fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
     1271c$$$        fsollw(i,nsrf) = sollw(i)
     1272         fsollw(i,nsrf) = sollw(i)
     1273     $      + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i,nsrf))
     1274         fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
     1275       ENDDO
     1276       ENDDO
    11111277
    11121278      fder = dlw
    … …  
    11161282     e            julien, rmu0,
    11171283     e            ok_veget, ocean, npas, nexca, ftsol,
    1118      $            soil_model,ftsoil, qsol,
     1284     $            soil_model,cdmmax, cdhmax, ftsoil, qsol,
    11191285     $            paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
    11201286     $            fluxlat,
    … …  
    16161782      enddo
    16171783
     1784cIM ISCCP simulator BEGIN
     1785      IF (ok_isccp) THEN
     1786cIM calcul tau. emi nuages convectifs
     1787      convfra(:,:)=rnebcon(:,:)
     1788      convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
     1789c     CALL newmicro (paprs, pplay,ok_newmicro,
     1790c    .            t_seri, cldliq, cldfra, cldtau, cldemi,
     1791c    .            cldh, cldl, cldm, cldt, cldq)
     1792      CALL newmicro (paprs, pplay,ok_newmicro,
     1793     .            t_seri, convliq, convfra, dtau_c, dem_c,
     1794     .            cldh_c, cldl_c, cldm_c, cldt_c, cldq_c)
     1795
     1796cIM calcul tau. emi nuages startiformes
     1797      CALL newmicro (paprs, pplay,ok_newmicro,
     1798     .            t_seri, cldliq, cldfra, dtau_s, dem_s,
     1799     .            cldh_s, cldl_s, cldm_s, cldt_s, cldq_s)
     1800cIM calcul diagramme (PC, tau) cf. ISCCP D
     1801c     seed=50
     1802c     seed=ran0(klon)
     1803cT1O3     
     1804c     top_height=1
     1805cT3O3
     1806c     top_height=3
     1807c     overlap=3
     1808cIM cf GCM     
     1809      cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
     1810
     1811cIM inversion des niveaux de pression ==> de haut en bas
     1812      DO k=1,klev
     1813       kinv=klev-k+1
     1814       DO i=1,klon
     1815        pfull(i,k)=pplay(i,kinv)
     1816c on met toutes les variables de Haut 2 Bas
     1817        qv(i,k)=q_seri(i,kinv)
     1818        cc(i,k)=cldtot(i,kinv)
     1819        conv(i,k)=rnebcon(i,kinv)
     1820        dtau_sH2B(i,k)=dtau_s(i,kinv)
     1821        dtau_cH2B(i,k)=dtau_c(i,kinv)
     1822        at(i,k)=t_seri(i,kinv)
     1823        dem_sH2B(i,k)=dem_s(i,kinv)
     1824        dem_cH2B(i,k)=dem_c(i,kinv)
     1825
     1826       ENDDO
     1827      ENDDO
     1828
     1829      DO k=1,klev+1
     1830       kinv=klev-k+2
     1831       DO i=1,klon
     1832        phalf(i,k)=paprs(i,kinv)
     1833       ENDDO
     1834      ENDDO
     1835
     1836c     open(99,file='tautab.bin',access='sequential',
     1837c    $     form='unformatted',status='old')
     1838c     read(99) tautab
     1839
     1840cIM210503
     1841      IF (debut) THEN
     1842      open(99,file='tautab.formatted', FORM='FORMATTED')
     1843      read(99,'(f30.20)') tautab
     1844      close(99)
     1845c
     1846      open(99,file='invtau.formatted',form='FORMATTED')
     1847      read(99,'(i10)') invtau
     1848      close(99)
     1849c
     1850       nsrf=3
     1851       DO nreg=1, nbreg
     1852       DO i=1, klon
     1853
     1854c       IF (debut) THEN
     1855         IF(rlon(i).LT.0.) THEN
     1856           rlonPOS(i)=rlon(i)+360.
     1857         ELSE
     1858           rlonPOS(i)=rlon(i) 
     1859         ENDIF
     1860c       ENDIF
     1861
     1862c       pct_ocean(i,nreg)=.FALSE.
     1863        pct_ocean(i,nreg)=0
     1864
     1865c      DO nsrf = 1, nbsrf
     1866
     1867c test si c'est 1 point d'ocean
     1868        IF(pctsrf(i,nsrf).EQ.1.) THEN
     1869
     1870         IF(nreg.EQ.1) THEN
     1871
     1872c TROP
     1873          IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
     1874c          pct_ocean(i,nreg)=.TRUE.
     1875           pct_ocean(i,nreg)=1
     1876          ENDIF
     1877
     1878c PACIFIQUE NORD
     1879          ELSEIF(nreg.EQ.2) THEN
     1880           IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
     1881            IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN
     1882c            pct_ocean(i,nreg)=.TRUE.
     1883             pct_ocean(i,nreg)=1
     1884            ENDIF
     1885           ENDIF
     1886c CALIFORNIE ST-CU
     1887         ELSEIF(nreg.EQ.3) THEN
     1888          IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
     1889           IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
     1890c           pct_ocean(i,nreg)=.TRUE.
     1891            pct_ocean(i,nreg)=1
     1892           ENDIF
     1893          ENDIF
     1894c HAWAI
     1895        ELSEIF(nreg.EQ.4) THEN
     1896         IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
     1897          IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
     1898c          pct_ocean(i,nreg)=.TRUE.
     1899           pct_ocean(i,nreg)=1
     1900          ENDIF
     1901         ENDIF
     1902c WARM POOL
     1903        ELSEIF(nreg.EQ.5) THEN
     1904         IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
     1905          IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
     1906c          pct_ocean(i,nreg)=.TRUE.
     1907           pct_ocean(i,nreg)=1
     1908          ENDIF
     1909         ENDIF
     1910        ENDIF !nbreg
     1911c TROP
     1912c        IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
     1913c         pct_ocean(i)=.TRUE.
     1914c         WRITE(*,*) 'pct_ocean =',i, rlon(i), rlat(i)
     1915c          ENDIF !lon
     1916c         ENDIF !lat
     1917
     1918        ENDIF !pctsrf
     1919c      ENDDO
     1920       ENDDO !klon
     1921       ENDDO !nbreg
     1922 
     1923cIM somme de toutes les nhistoW BEG
     1924      DO nreg = 1, nbreg
     1925      DO k = 1, kmaxm1
     1926      DO l = 1, lmaxm1
     1927      DO iw = 1, iwmax
     1928       nhistoWt(k,l,iw,nreg)=0.
     1929      ENDDO
     1930      ENDDO
     1931      ENDDO
     1932      ENDDO
     1933cIM somme de toutes les nhistoW END
     1934      ENDIF
     1935
     1936
     1937c     CALL ISCCP_CLOUD_TYPES(nlev,ncol,seed,pfull,phalf,qv,
     1938c    &     cc,conv,dtau_s,dtau_c,top_height,overlap,
     1939c    &     tautab,invtau,skt,emsfc_lw,at,dem_s,dem_c,fq_isccp,
     1940c    &     totalcldarea,meanptop,meantaucld,boxtau,boxptop)
     1941
     1942c     DO i=1, klon
     1943c     i=1
     1944c1011  CONTINUE
     1945c
     1946cIM on verifie les donnees de INPUT en dehors du simulateur ISCCP
     1947cIM 1D non-vectorise (!) pour qu'on gagne du temps ...
     1948cIM
     1949c BEGIN find unpermittable data.....
     1950!     ---------------------------------------------------!
     1951!     find unpermittable data.....
     1952!
     1953c     do 13 k=1,klev
     1954c ca prend trop de temps ??
     1955c     cldtot(:,:) = min(max(cldtot(:,:),0.),1.)
     1956c     rnebcon(:,:) = min(max(rnebcon(:,:),0.),1.)
     1957c     dtau_s(:,:) = max(dtau_s(:,:),0.)
     1958c     dem_s(:,:) = min(max(dem_s(:,:),0.),1.)
     1959c     dtau_c(:,:) = max(dtau_c(:,:),0.)
     1960c     dem_c(:,:) = min(max(dem_c(:,:),0.),1.)
     1961c ca prend trop de temps ??
     1962
     1963c           if (cldtot(i,k) .lt. 0.) then
     1964c               print *, ' error = cloud fraction less than zero'
     1965c               STOP
     1966c           end if
     1967c           if (cldtot(i,k) .gt. 1.) then
     1968c               print *, ' error = cloud fraction greater than 1'
     1969c               STOP
     1970c           end if
     1971c           if (rnebcon(i,k) .lt. 0.) then
     1972c               print *,
     1973c    &           ' error = convective cloud fraction less than zero'
     1974c               STOP
     1975c           end if
     1976c           if (rnebcon(i,k) .gt. 1.) then
     1977c               print *,
     1978c    &           ' error = convective cloud fraction greater than 1'
     1979c               STOP
     1980c           end if
     1981
     1982c           if (dtau_s(i,k) .lt. 0.) then
     1983c               print *,
     1984c    &           ' error = stratiform cloud opt. depth less than zero'
     1985c               STOP
     1986c           end if
     1987c           if (dem_s(i,k) .lt. 0.) then
     1988c               print *,
     1989c    &           ' error = stratiform cloud emissivity less than zero'
     1990c               STOP
     1991c           end if
     1992c           if (dem_s(i,k) .gt. 1.) then
     1993c               print *,
     1994c    &           ' error = stratiform cloud emissivity greater than 1'
     1995c               STOP
     1996c           end if
     1997
     1998c           if (dtau_c(i,k) .lt. 0.) then
     1999c               print *,
     2000c    &           ' error = convective cloud opt. depth less than zero'
     2001c               STOP
     2002c           end if
     2003c           if (dem_c(i,k) .lt. 0.) then
     2004c               print *,
     2005c    &           ' error = convective cloud emissivity less than zero'
     2006c               STOP
     2007c           end if
     2008c           if (dem_c(i,k) .gt. 1.) then
     2009c               print *,
     2010c    &           ' error = convective cloud emissivity greater than 1'
     2011c               STOP
     2012c           end if
     2013c13    continue
     2014
     2015!     ---------------------------------------------------!
     2016c
     2017c END   find unpermittable data.....
     2018cv2.2.1.1     DO i=1, klon
     2019c     i=1
     2020c     seed=i
     2021c
     2022cv3.4
     2023      if (debut) then
     2024        DO i=1, klon
     2025          seed(i)=i+100
     2026c         seed(i)=i+50
     2027        ENDDO
     2028      endif
     2029c     seed=aint(ran0(klon))
     2030c     CALL ISCCP_CLOUD_TYPES(klev,ncol,seed,pfull(i,:),phalf(i,:)
     2031cv2.2.1.1
     2032c     CALL ISCCP_CLOUD_TYPES(klev,ncol,seed(i),pfull(i,:),phalf(i,:)
     2033c    &     ,q_seri(i,:),
     2034c    &     cldtot(i,:),rnebcon(i,:),dtau_s(i,:),dtau_c(i,:),
     2035c    &     top_height,overlap,
     2036c    &     tautab,invtau,ztsol,emsfc_lw,t_seri(i,:),dem_s(i,:),
     2037c    &     dem_c(i,:),
     2038c    &     fq_isccp(i,:,:),
     2039c    &     totalcldarea(i),meanptop(i),meantaucld(i),
     2040c    &     boxtau(i,:),boxptop(i,:))
     2041cv2.2.1.1
     2042cv3.4
     2043      debug=0
     2044      debugcol=0
     2045cIM260503
     2046c o500 ==> distribution nuage ftion du regime dynamique
     2047      DO i=1, klon
     2048       o500(i)=omega(i,8)*864.
     2049c      PRINT*,'pphi8 ',pphi(i,8),'zphi8,11,12',zphi(i,8),
     2050c    & zphi(i,11),zphi(i,12)
     2051      ENDDO
     2052
     2053c axe vertical pour les differents niveaux des histogrammes
     2054c     DO iw=1, iwmax
     2055c       zx_o500(iw)=wmin+(iw-1./2.)*pas_w
     2056c     ENDDO
     2057c     PRINT*,' phys AVANT seed(3361)=',seed(3361)
     2058      CALL ISCCP_CLOUD_TYPES(
     2059     &     debug,
     2060     &     debugcol,
     2061     &     klon,
     2062     &     sunlit,
     2063     &     klev,
     2064     &     ncol,
     2065     &     seed,
     2066     &     pfull,
     2067     &     phalf,
     2068c var de bas en haut ==> PB !
     2069c    &     q_seri,
     2070c    &     cldtot,
     2071c    &     rnebcon,
     2072c    &     dtau_s,
     2073c    &     dtau_c,
     2074c var de Haut en Bas BEG
     2075     &     qv, cc, conv, dtau_sH2B, dtau_cH2B,
     2076c var de Haut en Bas END
     2077     &     top_height,
     2078     &     overlap,
     2079     &     tautab,
     2080     &     invtau,
     2081     &     ztsol,
     2082     &     emsfc_lw,
     2083c var de bas en haut ==> PB !
     2084c    &     t_seri,
     2085c    &     dem_s,
     2086c    &     dem_c,
     2087c var de Haut en Bas BEG
     2088     &     at, dem_sH2B, dem_cH2B,
     2089cIM260503
     2090c    &     o500, pct_ocean,
     2091c var de Haut en Bas END
     2092     &     fq_isccp,
     2093     &     totalcldarea,
     2094     &     meanptop,
     2095     &     meantaucld,
     2096     &     boxtau,
     2097     &     boxptop)
     2098c    &     boxptop,
     2099cIM 260503
     2100c    &     histoW,
     2101c    &     nhistoW   
     2102c    &)
     2103
     2104cIM 200603
     2105c     PRINT*,'physiq fq_isccp(6,1,1)',fq_isccp(6,1,1)
     2106       
     2107cIM 200603
     2108cIM somme de toutes les nhistoW BEG
     2109c     DO k = 1, kmaxm1
     2110c     DO l = 1, lmaxm1
     2111c     DO iw = 1, iwmax
     2112c     nhistoWt(k,l,iw)=nhistoWt(k,l,iw)+nhistoW(k,l,iw)
     2113ccc      IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
     2114c      IF(nhistoWt(k,l,iw).NE.0.) THEN
     2115c       PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
     2116c      ENDIF
     2117c     ENDDO
     2118c     ENDDO
     2119c     ENDDO
     2120cIM somme de toutes les nhistoW END
     2121c     PRINT*,' phys APRES seed(3361)=',seed(3361)
     2122cv3.4
     2123c     i=i+1
     2124c     IF(i.LE.klon) THEN
     2125c      GOTO 1011
     2126c     ENDIF
     2127cv2.2.1.1     ENDDO
     2128
     2129c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
     2130c     minfq3d=100.
     2131c     maxfq3d=0.
     2132cIM calcul des correspondances entre la grille physique et
     2133cIM la grille dynamique
     2134c     DO i=1, klon
     2135c       grid_phys(i)=i
     2136c       PRINT*,'i, grid_phys',i,grid_phys(i)
     2137c     ENDDO
     2138c     CALL gr_fi_dyn(1,klon,iimp1,jjmp1,grid_phys,grid_dyn)
     2139c     DO j=1, jjmp1
     2140c       DO i=1, iimp1
     2141c        PRINT*,'i,j grid_dyn ',i,j,grid_dyn(i,j)
     2142c       ENDDO
     2143c     ENDDO
     2144c
     2145      DO l=1, lmax
     2146       DO k=1, kmax
     2147cIM grille physique ==> grille ecriture 2D (iim,jjmp1)
     2148c
     2149        DO i=1, iim
     2150          fq4d(i,1,k,l)=fq_isccp(1,k,l)
     2151cc         PRINT*,'first j=1 i =',i
     2152        ENDDO
     2153        DO j=2, jjm
     2154          DO i=1, iim
     2155cERROR ??         ig=i+iim*(j-1)
     2156          ig=i+1+(j-2)*iim
     2157cc         PRINT*,'i =',i,'j =',j,'ig=',ig
     2158          fq4d(i,j,k,l)=fq_isccp(ig,k,l)             
     2159         ENDDO
     2160        ENDDO
     2161        DO i=1, iim
     2162          fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l)
     2163cc         PRINT*,'last jjmp1 i =',i
     2164        ENDDO
     2165        IF(debut) THEN
     2166        DO j=1, jjmp1
     2167          DO i=1, iim
     2168            IF(j.GE.2.AND.j.LE.jjm) THEN
     2169              igfi2D(i,j)=i+1+(j-2)*iim
     2170c             PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
     2171            ELSEIF(j.EQ.1) THEN
     2172              igfi2D(i,j)=1
     2173c             PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
     2174            ELSEIF(j.EQ.jjmp1) THEN
     2175              igfi2D(i,j)=klon
     2176c             PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
     2177            ENDIF
     2178          ENDDO
     2179        ENDDO
     2180        ENDIF
     2181c       STOP
     2182c
     2183c       CALL gr_fi_ecrit(1,klon,iim,jjmp1,fq_isccp(:,k,l),
     2184c    $       fq4d(:,:,k,l))
     2185       ENDDO
     2186      ENDDO
     2187      DO l=1, lmax
     2188       fq4d(:,:,8,l)=-1.e+10
     2189       fq4d(:,:,l,8)=-1.e+10
     2190      ENDDO
     2191      DO l=1, lmax
     2192       DO k=1, kmax 
     2193        DO j=1, jjmp1
     2194         DO i=1, iim
     2195
     2196c inversion TAU ?!
     2197c         ni=(i-1)*lmax+l
     2198c         nj=(j-1)*kmax+kmax-k+1
     2199c
     2200c210503 inversion en PC ==> pas besoin !!!
     2201c         ni=(i-1)*lmax+lmax-l+1
     2202c         nj=(j-1)*kmax+k
     2203c
     2204c210503
     2205          ni=(i-1)*lmax+l
     2206          nj=(j-1)*kmax+k
     2207 
     2208c210503         if(k.EQ.8) then
     2209c          fq4d(i,j,8,l)=-1.e+10
     2210c         endif
     2211
     2212c210503         if(l.EQ.8) then
     2213c          fq4d(i,j,k,8)=-1.e+10
     2214c         endif
     2215
     2216          fq3d(ni,nj)=fq4d(i,j,k,l)
     2217
     2218c         if(fq3d(ni,nj).lt.0.) then
     2219c          print*,' fq3d LT ZERO ',ni,nj,fq3d(ni,nj)
     2220c         endif
     2221c         if(fq3d(ni,nj).gt.100.) then
     2222c          print*,' fq3d GT 100 ',ni,nj,fq3d(ni,nj)
     2223c         endif
     2224c max & min fq3d
     2225c         if(fq3d(ni,nj).gt.maxfq3d) maxfq3d=fq3d(ni,nj)
     2226c         if(fq3d(ni,nj).lt.minfq3d) minfq3d=fq3d(ni,nj)
     2227         
     2228         ENDDO
     2229        ENDDO
     2230c       fq4d(:,:,8,l)=-1.e+10
     2231c       fq4d(:,:,k,8)=-1.e+10
     2232c       k=k+1
     2233c       if(k.LE.kmax) then
     2234c        goto 1022
     2235c       endif
     2236       ENDDO
     2237c      l=l+1
     2238c      if(l.LE.lmax) then
     2239c       goto 1021
     2240c      endif
     2241      ENDDO
     2242
     2243c     print*,' minfq3d=',minfq3d,' maxfq3d=',maxfq3d
     2244c
     2245c calculs statistiques distribution nuage ftion du regime dynamique
     2246c     DO i=1, klon
     2247c!      o500(i)=omega(i,9)*864.
     2248c!      PRINT*,' o500=',o500(i),' pphi(9)=',pphi(i,9)
     2249c       o500(i)=omega(i,8)*864.
     2250cc      PRINT*,' pphi(8)',pphi(i,8),'pphi(11)',pphi(i,11),
     2251cc    .'pphi(12)',pphi(i,12)
     2252cc      PRINT*,' zphi8,11,12=',zphi(i,8),zphi(i,11),zphi(i,12)
     2253cc     PRINT*,' o500',o500(i),' w500',w500(i)
     2254c     ENDDO
     2255
     2256c axe vertical pour les differents niveaux des histogrammes
     2257c     DO iw=1, iwmax
     2258c       zx_o500(iw)=wmin+(iw-1./2.)*pas_w
     2259c     ENDDO
     2260
     2261
     2262c Ce calcul doit etre fait a partir de valeurs mensuelles ??
     2263cc     CALL histo_o500_pctau(o500,fq4d,histoW)
     2264cc     CALL histo_o500_pctau(paire,pctsrf,o500,fq4d,histoW)
     2265cc     CALL histo_o500_pctau(pct_ocean,rlat,o500,fq4d,histoW)
     2266ccOK ???     CALL histo_o500_pctau(pct_ocean,o500,fq4d,histoW)
     2267c     CALL histo_o500_pctau(klon,pct_ocean,o500,fq4d,histoW,nhistoW)
     2268c     CALL histo_o500_pctau(klon,pct_ocean,o500,fq_isccp,
     2269      CALL histo_o500_pctau(nbreg,pct_ocean,o500,fq_isccp,
     2270     &histoW,nhistoW)
     2271c
     2272cIM somme de toutes les nhistoW BEG
     2273      DO nreg=1, nbreg
     2274      DO k = 1, kmaxm1
     2275      DO l = 1, lmaxm1
     2276      DO iw = 1, iwmax
     2277       nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
     2278     & nhistoW(k,l,iw,nreg)
     2279ccc      IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
     2280c      IF(nhistoWt(k,l,iw).NE.0.) THEN
     2281c       PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
     2282c      ENDIF
     2283      ENDDO
     2284      ENDDO
     2285      ENDDO
     2286      ENDDO
     2287cIM somme de toutes les nhistoW END
     2288c
     2289c     IF(lafin) THEN   
     2290c     DO nreg=1, nbreg
     2291c     DO iw=1, iwmax
     2292c     DO l=1,lmaxm1
     2293c     DO k=1,kmaxm1
     2294c      IF(histoW(k,l,iw,nreg).NE.0.) then
     2295c        PRINT*,'physiq H nH',k,l,iw,
     2296c    &       histoW(k,l,iw,nreg),
     2297c    &       nhistoW(k,l,iw,nreg),nhistoWt(k,l,iw,nreg)
     2298c      ENDIF
     2299c     ENDDO
     2300c     ENDDO
     2301c     ENDDO
     2302c     ENDDO
     2303cIM verif fq_isccp, fq4d, fq3d
     2304c     DO l=1, lmaxm1
     2305c       DO k=1,kmaxm1
     2306c     i=74
     2307c     j=36
     2308c     DO j=1, jjmp1
     2309c      DO i=1, iim
     2310c       DO l=1, lmaxm1
     2311c         WRITE(*,'(a,3i4,7f10.4)')
     2312c    &    'fq_isccp,j,i,l=',j,i,l,
     2313c    &    (fq_isccp(igfi2D(i,j),k,l),k=1,kmaxm1)
     2314c         WRITE(*,'(a,3i4,7f10.4)')
     2315c    &    'fq4d,j,i,l=',j,i,l,(fq4d(i,j,k,l),k=1,kmaxm1)
     2316c       ENDDO
     2317c      ENDDO
     2318c     ENDDO
     2319c     ni1=(i-1)*8+1
     2320c     ni2=i*8
     2321c     nj1=(j-1)*8+1
     2322c     nj2=j*8
     2323c     DO ni=ni1,ni2
     2324c     WRITE(*,'(a,2i4,7f10.4)')
     2325c    &     'fq3d, ni,nj=',ni,nj,
     2326c    &      (fq3d(ni,nj),nj=nj1,nj2)
     2327c     ENDDO
     2328c     ENDIF
     2329
     2330c     DO iw=1, iwmax
     2331c      DO l=1,lmaxm1
     2332c       DO k=1,kmaxm1
     2333c        PRINT*,' iw,l,k,nhistoW=',iw,l,k,nhistoW(k,l,iw)
     2334c       ENDDO
     2335c      ENDDO
     2336c     ENDDO
     2337
     2338c       DO iw=1, iwmax
     2339c        DO l=1, lmaxm1
     2340c         linv=lmaxm1-l+1
     2341c         DO k=1, kmaxm1
     2342c         histoWinv(k,l,iw)=histoW(iw,k,l)
     2343c       ENDDO
     2344c      ENDDO
     2345c     ENDDO
     2346c
     2347c pb syncronisation ?? : 48 * 30 * 7 (jour1) + 48* 29 * 7 (jour suivant)
     2348c
     2349
     2350
     2351      ENDIF !ok_isccp
     2352cIM ISCCP simulator END
     2353
    16182354c   On prend la somme des fractions nuageuses et des contenus en eau
    16192355      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
    … …  
    17172453cccIMs             topsw0,toplw0,solsw0,sollw0)
    17182454     s             topsw0,toplw0,solsw0,sollw0,
    1719      s             ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
     2455     s             swdn0, swdn, swup0, swup     )
    17202456      itaprad = 0
    17212457      ENDIF
    … …  
    19682704cIM cf. FH     slp(:) = paprs(:,1)*exp(pphis(:)/(289.*t_seri(:,1)))
    19692705      slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
     2706c     PRINT*,' physiq slp ',slp(2185),paprs(2185,1),pphis(2185),
     2707c    .       RD,t_seri(2185,1)
     2708c
     2709ccc prw = eau precipitable
     2710      DO i = 1, klon
     2711       prw(i) = 0.
     2712      DO k = 1, klev
     2713       prw(i) = prw(i) +
     2714     .          q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
     2715      ENDDO
     2716c      PRINT*,' i ',i,' prw',prw(i)
     2717      ENDDO
    19702718c
    19712719
    … …  
    19732721c   Ecriture des sorties
    19742722c=============================================================
     2723
     2724#ifdef histISCCP
     2725#include "write_histISCCP.h"
     2726#endif
    19752727
    19762728#ifdef histhf
    … …  
    20242776         CALL phyredem ("restartphy.nc",dtime,radpas,
    20252777     .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
    2026      .      fsnow, falbe, fevap, rain_fall, snow_fall,
     2778     .      fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
    20272779     .      solsw, sollwdown,dlw,
    20282780     .      radsol,frugs,agesno,
Note: See TracChangeset for help on using the changeset viewer.