Changeset 175 for trunk/LMDZ.TITAN/libf


Ignore:
Timestamp:
Jun 28, 2011, 12:56:50 PM (13 years ago)
Author:
slebonnois
Message:

S.LEBONNOIS:

  • Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
  • Quelques ajouts dans la doc.
Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
22 added
2 deleted
23 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/clesphys.h

    r104 r175  
    1111       INTEGER nbapp_rad, nbapp_chim, iflag_con, iflag_ajs
    1212       REAL    ecritphy
     13       INTEGER lev_histmth, lev_histday
    1314       REAL    solaire
     15
     16! Parametres pour PBL:
    1417       REAL    z0, lmixmin
    15        REAL    ksta, inertie
     18       REAL    ksta
    1619       LOGICAL ok_kzmin
    17        INTEGER lev_histmth, lev_histday
     20
     21! Parametres surface:
     22       REAL    inertie
     23
     24! Parametres Chimie:
    1825       logical chimi,ylellouch,hcnrad
    1926       integer vchim,aerprod,htoh2
    20        integer microfi,cutoff
    21        real    tx,tcorrect
     27       
     28! Parametres Microphysique:
     29       integer microfi,cutoff,clouds
     30       real    tx,tcorrect,p_prodaer
     31       real    xnuf
    2232
    2333
    24        COMMON/clesphys/cycle_diurne, soil_model,                        &
    25      &     ok_orodr, ok_orolf, ok_gw_nonoro, nbapp_rad, nbapp_chim      &
    26      &     , ecritphy                                                   &
    27      &     , iflag_con, iflag_ajs, solaire, z0, lmixmin, ksta           &
    28      &     , ok_kzmin, lev_histmth, lev_histday                         &
    29      &     , inertie, chimi,vchim,aerprod,htoh2,ylellouch,hcnrad        &
    30      &     , microfi,cutoff,tx,tcorrect
     34       COMMON/clesphys_i/                                               &
     35     &     nbapp_rad, nbapp_chim, iflag_con, iflag_ajs,                 &
     36     &     lev_histmth, lev_histday, vchim,aerprod,htoh2,               &
     37     &     microfi,cutoff,clouds
    3138
     39       COMMON/clesphys_r/                                               &
     40     &     ecritphy, solaire, z0, lmixmin, ksta, inertie,               &
     41     &     tx,tcorrect,p_prodaer,xnuf
     42
     43       COMMON/clesphys_l/cycle_diurne, soil_model,                      &
     44     &     ok_orodr, ok_orolf, ok_gw_nonoro, ok_kzmin,                  &
     45     &     chimi,ylellouch,hcnrad
     46

  • trunk/LMDZ.TITAN/libf/phytitan/conf_phys.F90

    r97 r175  
    367367
    368368!
     369!Config Key  = p_prodaer
     370!Config Desc = pressure level for aerosol production (in Pa)
     371!Config Def  = 1.
     372!Config Help =
     373!
     374  p_prodaer = 1.
     375  call getin('p_prodaer',p_prodaer)
     376!
    369377!Config Key  = cutoff
    370378!Config Desc =
    … …  
    375383  call getin('cutoff',cutoff)
    376384
     385!
     386!Config Key  = clouds
     387!Config Desc = activation des nuages
     388!Config Def  = 1
     389!Config Help =
     390!
     391  clouds = 1
     392  call getin('clouds',clouds)
     393  if (microfi.ne.1) clouds = 0      ! On  ne fait pas de nuages sans microphysique !
     394
     395!
     396!Config Key  = xnuf
     397!Config Desc = fraction nuageuse
     398!Config Def  = 0.5
     399!Config Help =
     400!
     401  xnuf = 0.5
     402  call getin('xnuf',xnuf)
     403  xnuf = amax1(xnuf,0.1)          ! On garde au minimum 10% de nuages.
     404  if (clouds.eq.0) xnuf = 0.      ! Si il n'y pas de nuages, on ne met pas de fraction
     405                                  ! nuageuse -> permet de retomber sur le TR habituel.
    377406!
    378407!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    448477  write(numout,*)' tx = ', tx
    449478  write(numout,*)' tcorrect = ', tcorrect
     479  write(numout,*)' p_prodaer = ', p_prodaer
    450480  write(numout,*)' cutoff = ', cutoff
     481  write(numout,*)' clouds = ', clouds
     482  write(numout,*)' xnuf = ', xnuf
    451483  write(numout,*)' lev_histmth = ',lev_histmth
    452484  write(numout,*)' lev_histday = ',lev_histday

  • trunk/LMDZ.TITAN/libf/phytitan/cooling.F

    r121 r175  
    1       SUBROUTINE COOLING(NL,PRESS,TEMP,Z,Q0,lwnet,pfluxi)
     1      SUBROUTINE COOLING(NG,NL,PRESS,TEMP,Z,Q0,lwnet,pfluxi,icld)
    22
    33c=======================================================================
    … …  
    6666c   ----------
    6767
    68       INTEGER NL
    69       REAL PRESS(ngrid,nl),TEMP(ngrid,nl)
    70       REAL Z(ngrid,nl),Q0(ngrid,nl-1)
    71       REAL lwnet(ngrid,nl),UBARI2
    72       real pfluxi(ngrid)
     68      INTEGER NG,NL,icld
     69      REAL PRESS(NG,NL),TEMP(NG,NL)
     70      REAL Z(NG,NL),Q0(NG,NL-1)
     71      REAL lwnet(NG,NL),UBARI2
     72      real pfluxi(NG)
    7373
    7474
    … …  
    7777
    7878C DTAU IS PASSED EN-MASS, SO ITS DEMENSIONS ARE CRITICAL
    79       COMMON /IRTAUS/ DTAUI(ngrid,NLAYER,NSPECI)
    80       REAL dtaui
     79      REAL dtaui(ngrid,NLAYER,NSPECI)
     80      REAL dtauip(ngrid,NLAYER,NSPECI)
     81      COMMON /IRTAUS/ dtaui,dtauip
    8182
    8283c   Local:
    … …  
    8889      REAL DW,WAVEN,TJ,BSURF,QOUT,QIN,eff_g,COLDEN
    8990
    90       INTEGER ig,K,J,NLEVEL,I,L
     91      INTEGER ig,K,J,I,L
    9192
    9293c     EXTERNAL PLNCK
    … …  
    139140      UBARI2=1./1.66
    140141      UBARI2=UBARI
    141       NLEVEL=NL
    142142
    143143C ZERO THE NET FLUXES
    … …  
    161161          B0 = 0.
    162162
    163           DO J=1,NLEVEL-1
    164              DO ig=1,ngrid
     163          DO J=1,NL-1
     164             DO ig=1,NG
    165165                TJ=TEMP(ig,J)
    166166
    … …  
    182182             ENDDO
    183183
    184              DO ig=1,ngrid
    185              zz4=EXP(-DTAUI(ig,J,K)/UBARI2)
    186                 EM(ig,J)=zz4
    187             ENDDO
     184             IF (ICLD.EQ.1) THEN
     185               DO ig=1,NG
     186                 zz4=EXP(-DTAUI(ig,J,K)/UBARI2)
     187                 EM(ig,J)=zz4
     188               ENDDO
     189             ELSE
     190               DO ig=1,NG
     191                 zz4=EXP(-DTAUIP(ig,J,K)/UBARI2)
     192                 EM(ig,J)=zz4
     193               ENDDO
     194             ENDIF
    188195          ENDDO
    189196
    … …  
    197204           FUPIS=0.
    198205
    199         DO 2220 J=1,NLEVEL-1
    200            DO 2230 ig=1,ngrid
     206        DO 2220 J=1,NL-1
     207           DO 2230 ig=1,NG
    201208              FDI(ig,J+1) = FDI(ig,J)*EM(ig,J) + 2.*RPI*UBARI*
    202209     &        B0(ig,J)*(1.-EM(ig,J))
    … …  
    218225C UPWARD FLUXES: SURFACE EMISSIONS
    219226
    220         DO 2310 ig=1,ngrid
     227        DO 2310 ig=1,NG
    221228          PLNCK=0.
    222229          DO I=-2,2,1
    223230             WAVNUM=WAVEN + I*zz1
    224              zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,NLEVEL))
     231             zz2=EXP(-1.4388 * WAVNUM/TEMP(ig,NL))
    225232             zz3=WAVNUM*WAVNUM*WAVNUM
    226233             PLNCK=PLNCK+1.191E-5* zz3*zz2/(1.-zz2)
    227234          ENDDO
    228 c          BSURF=PLNCK( WAVEN, TEMP(ig,NLEVEL), DW)
     235c          BSURF=PLNCK( WAVEN, TEMP(ig,NL), DW)
    229236           BSURF=.2*PLNCK
    230         FUPI(ig,NLEVEL)=BSURF * 2.*RPI*UBARI + RSFI*FDI(ig,NLEVEL)
    231         FUPIS(ig,NLEVEL,K)=BSURF*2.*RPI*UBARI+RSFI*FDIS(ig,NLEVEL,K)
     237        FUPI(ig,NL)=BSURF * 2.*RPI*UBARI + RSFI*FDI(ig,NL)
     238        FUPIS(ig,NL,K)=BSURF*2.*RPI*UBARI+RSFI*FDIS(ig,NL,K)
    2322392310    CONTINUE
    233240c     write(*,*)
    234 c     write(*,*) 'cooling : FUPI/NLEVEL =' ,
    235 c    & ((FUPI(i,l),l=nl,nl),i=1,ngrid)
    236 c     write(*,*)
    237 c     write(*,*) 'cooling : FDI/NLEVEL =' ,
    238 c    & ((FDI(i,l),l=nl,nl),i=1,ngrid)
    239 
    240         DO 2320 J=NLEVEL-1,1,-1
    241            DO 2330 ig=1,ngrid
     241c     write(*,*) 'cooling : FUPI/NL =' ,
     242c    & ((FUPI(i,l),l=nl,nl),i=1,NG)
     243c     write(*,*)
     244c     write(*,*) 'cooling : FDI/NL =' ,
     245c    & ((FDI(i,l),l=nl,nl),i=1,NG)
     246
     247        DO 2320 J=NL-1,1,-1
     248           DO 2330 ig=1,NG
    242249              FUPI(ig,J) = FUPI(ig,J+1)*EM(ig,J) + 2.*RPI*UBARI*
    243250     &        B0(ig,J)*(1.-EM(ig,J))
    … …  
    258265c compute the downward IR flux at the surface:
    259266c
    260           DO 3520 ig=1,ngrid
    261              pfluxi(ig)=pfluxi(ig)+ DWNI(K)*FDI(ig,NLEVEL)
     267          DO 3520 ig=1,NG
     268             pfluxi(ig)=pfluxi(ig)+ DWNI(K)*FDI(ig,NL)
    2622693520      CONTINUE
    263270
    264271c compute the net IR flux, (+) upward:
    265272c
    266           DO J=1,NLEVEL
    267           DO ig=1,ngrid
     273          DO J=1,NL
     274          DO ig=1,NG
    268275             lwnet(ig,J)= lwnet(ig,J)+ DWNI(K)*(FUPI(ig,J)-FDI(ig,J))
    269276          ENDDO
    270277          ENDDO
    271278         
    272           DO 3210 J=1,NLEVEL-1
    273              DO 3220 ig=1,ngrid
     279          DO 3210 J=1,NL-1
     280             DO 3220 ig=1,NG
    274281                QOUT=FUPI(ig,J) + FDI(ig,J+1)   ! OUT OF LAYER
    275282                QIN =FDI(ig,J)  + FUPI(ig,J+1)  ! INTO LAYER
    … …  
    297304
    298305c   convertion erg/cm2 -> J/m2
    299       DO 3550 ig=1,ngrid
     306      DO 3550 ig=1,NG
    300307         pfluxi(ig)  = 1.e-3*pfluxi(ig)
    301308         lwnet(ig,:) = 1.e-3*lwnet(ig,:)
    … …  
    307314C COMPUTE THE BASELINE COOLING RATE
    308315
    309        DO 3000 J=1,NLEVEL-1
     316       DO 3000 J=1,NL-1
    310317C TURN THE Q'S INTO TIMESCALES.....
    311           DO 3300 ig=1,ngrid
     318          DO 3300 ig=1,NG
    312319             eff_g = RG*(RA/(RA+Z(ig,J)))**2 ! 10% DIFF AT 1 MBAR
    313320             COLDEN = RHOP*(PRESS(ig,J+1)-PRESS(ig,J))/eff_g

  • trunk/LMDZ.TITAN/libf/phytitan/effg.F

    r3 r175  
    11      FUNCTION EFFG(Z)
    2       DATA G/1.35/
    3       DATA RPLANT/2575./
    4       EFFG = G * (RPLANT/(RPLANT + Z ) )**2
     2#include "YOMCST.h"
     3      EFFG = RG * (RA/(RA + Z ) )**2
    54      RETURN
    65      END

  • trunk/LMDZ.TITAN/libf/phytitan/heating.F

    r104 r175  
    1        SUBROUTINE heating(dist,rmu0,fract,sol_htg,swnet)
     1       SUBROUTINE heating(dist,rmu0,fract,sol_htg,swnet,icld)
    22
    33
    … …  
    1616c rmu0-----input-R- cosinus de l'angle zenithal
    1717c fract----input-R- duree d'ensoleillement normalisee
     18c icld-----input-I- calcul avec nuages.
    1819c        p(klon,nl)    pressure (level)
    1920c
    … …  
    4647
    4748      real dist, rmu0(klon), fract(klon)
     49      integer icld
    4850
    4951      real sol_htg(klon,klev)
    … …  
    9698               ubar0=rmu0(ig)
    9799
    98                CALL sfluxv(iprint,ig,dist)           ! #3
     100               CALL sfluxv(iprint,ig,dist,icld)           ! #3
    99101
    100102               fnetv(ig,:) = fnetv(ig,:) *fract(ig)   ! >0 vers le haut

  • trunk/LMDZ.TITAN/libf/phytitan/ini_histday.h

    r106 r175  
    1 !
    2 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ini_histday.h,v 1.2 2005/01/27 10:06:12 fairhead Exp $
    3 !
    41      IF (ok_journe) THEN
    52c
    … …  
    104101     .                "ave(X)", zsto,zout)
    105102c
    106       ENDIF !lev_histday.GE.2
    107 c
    108 c-------------------------------------------------------
    109       IF(lev_histday.GE.3) THEN
    110 c
    111103cccccccccccccccccc  Tracers
    112104c
    113105         if (iflag_trac.eq.1) THEN
    114           if (microfi.eq.1) then
    115            DO iq=1,nmicro
    116          CALL histdef(nid_day, tname(iq), ttext(iq), "n/m2",
    117      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
    118      .                "ave(X)", zsto,zout)
    119            ENDDO
     106          if (microfi.ge.1) then
     107c           DO iq=1,nmicro
     108c             CALL histdef(nid_day, tname(iq), ttext(iq), "n/m2",
     109c     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     110c     .                "ave(X)", zsto,zout)
     111c           ENDDO
     112             CALL histdef(nid_day, "qaer","nb tot aer" , "n/m2",
     113     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     114     .                    "ave(X)", zsto,zout)
     115             CALL histdef(nid_day, "qnoy","nb tot noy" , "n/m2",
     116     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     117     .                    "ave(X)", zsto,zout)
     118             CALL histdef(nid_day, "qgl1","V tot gl1" , "m3/m2",
     119     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     120     .                    "ave(X)", zsto,zout)
     121             CALL histdef(nid_day, "qgl2","V tot gl2" , "m3/m2",
     122     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     123     .                    "ave(X)", zsto,zout)
     124             CALL histdef(nid_day, "qgl3","V tot gl3" , "m3/m2",
     125     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     126     .                    "ave(X)", zsto,zout)
     127c--------------
     128c ----- SATURATION ESP NUAGES
     129             if (clouds.eq.1) then
     130               CALL histdef(nid_day,"ch4sat", "saturation CH4", "--",
     131     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     132     .                "ave(X)", zsto,zout)
     133               CALL histdef(nid_day,"c2h6sat", "saturation C2H6", "--",
     134     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     135     .                "ave(X)", zsto,zout)
     136               CALL histdef(nid_day,"c2h2sat", "saturation C2H2", "--",
     137     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     138     .                "ave(X)", zsto,zout)
     139c --------------
     140c ----- RESERVOIR DE SURFACE
     141               CALL histdef(nid_day, "reserv", "Reservoir surface","m",
     142     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     143     .                "ave(X)", zsto,zout)
     144c --------------
     145c ----- PRECIPITATIONS (precipitations cumulatives)
     146               CALL histdef(nid_day,"prech4","Precip CH4","m",
     147     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     148     .                "ave(X)", zsto,zout)
     149               CALL histdef(nid_day,"prec2h6","Precip C2H6",
     150     .                "m",iim,jjmp1,nhori, 1,1,1, -99, 32,
     151     .                "ave(X)", zsto,zout)
     152               CALL histdef(nid_day,"prec2h2","Precip C2H2",
     153     .                "m",iim,jjmp1,nhori, 1,1,1, -99, 32,
     154     .                "ave(X)", zsto,zout)
     155c --------------
     156c ----- FLUX GLACE
     157               CALL histdef(nid_day,"flxgl1", "flux gl CH4",
     158     .              "kg/m2/s",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     159     .              "ave(X)", zsto,zout)
     160               CALL histdef(nid_day,"flxgl2", "flux gl C2H6",
     161     .              "kg/m2/s",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     162     .              "ave(X)", zsto,zout)
     163               CALL histdef(nid_day,"flxgl3", "flux gl C2H2",
     164     .              "kg/m2/s",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     165     .              "ave(X)", zsto,zout)
     166c --------------
     167c ----- RAYON DES GOUTTES
     168               CALL histdef(nid_day,"rcldbar", "rayon moyen goutte",
     169     .                "m",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     170     .                "ave(X)", zsto,zout)
     171             endif
    120172          endif
     173c --------------
     174c ----- TRACEURS CHIMIQUES
    121175          if (nmicro.lt.nqmax) then
    122176           DO iq=nmicro+1,nqmax
    … …  
    128182         endif
    129183c
     184      ENDIF !lev_histday.GE.2
     185c
     186c-------------------------------------------------------
     187      IF(lev_histday.GE.3) THEN
     188c
    130189cccccccccccccccccc  Radiative transfer
    131190c
    … …  
    154213     .                32, "ave(X)", zsto1,zout)
    155214c
     215c --------------
     216c ----- OPACITE BRUME
    156217         DO k=7,NSPECV,10
    157218           write(str1,'(i2.2)') k
    … …  
    161222         ENDDO
    162223c
     224         DO k=8,NSPECI,10
     225           write(str1,'(i2.2)') k
     226         CALL histdef(nid_day,"thi"//str1,"Haze Opa IR",
     227     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     228     .                "ave(X)",zsto1,zout)
     229         ENDDO
     230c
     231c --------------
     232c ----- EXTINCTION BRUME
    163233         DO k=7,NSPECV,10
    164234           write(str1,'(i2.2)') k
    … …  
    168238         ENDDO
    169239c
     240         DO k=8,NSPECI,10
     241           write(str1,'(i2.2)') k
     242         CALL histdef(nid_day,"khi"//str1,"Haze ext IR ",
     243     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     244     .                "ave(X)",zsto1,zout)
     245         ENDDO
     246c
     247c --------------
     248c ----- OPACITE GAZ
    170249         DO k=7,NSPECV,10
    171250           write(str1,'(i2.2)') k
    172          CALL histdef(nid_day,"tgv"//str1,"Haze Opa Vis",
    173      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    174      .                "ave(X)",zsto1,zout)
    175          ENDDO
    176 c
     251         CALL histdef(nid_day,"tgv"//str1,"Gas Opa Vis",
     252     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     253     .                "ave(X)",zsto1,zout)
     254         ENDDO
     255c
     256         DO k=8,NSPECI,10
     257           write(str1,'(i2.2)') k
     258         CALL histdef(nid_day,"tgi"//str1,"Gas Opa IR",
     259     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     260     .                "ave(X)",zsto1,zout)
     261         ENDDO
     262c
     263c --------------
     264c ----- EXTINCTION GAZ
    177265         DO k=7,NSPECV,10
    178266           write(str1,'(i2.2)') k
    179          CALL histdef(nid_day,"kgv"//str1,"Haze ext Vis ",
     267         CALL histdef(nid_day,"kgv"//str1,"Gas ext Vis ",
    180268     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    181269     .                "ave(X)",zsto1,zout)
    … …  
    184272         DO k=8,NSPECI,10
    185273           write(str1,'(i2.2)') k
    186          CALL histdef(nid_day,"thi"//str1,"Haze Opa IR",
    187      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    188      .                "ave(X)",zsto1,zout)
    189          ENDDO
    190 c
    191          DO k=8,NSPECI,10
    192            write(str1,'(i2.2)') k
    193          CALL histdef(nid_day,"khi"//str1,"Haze ext IR ",
    194      .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    195      .                "ave(X)",zsto1,zout)
    196          ENDDO
    197 c
    198          DO k=8,NSPECI,10
    199            write(str1,'(i2.2)') k
    200          CALL histdef(nid_day,"tgi"//str1,"Haze Opa IR",
    201      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    202      .                "ave(X)",zsto1,zout)
    203          ENDDO
    204 c
    205          DO k=8,NSPECI,10
    206            write(str1,'(i2.2)') k
    207          CALL histdef(nid_day,"kgi"//str1,"Haze ext IR ",
    208      .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    209      .                "ave(X)",zsto1,zout)
    210          ENDDO
     274         CALL histdef(nid_day,"kgi"//str1,"Gas ext IR ",
     275     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     276     .                "ave(X)",zsto1,zout)
     277         ENDDO
     278c
     279c --------------
     280c ----- OPACITE NUAGES
     281         if (clouds.eq.1) then
     282           CALL histdef(nid_day,"tcld","Cld Opa proxy",
     283     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     284     .                "ave(X)",zsto,zout)
     285c
     286c --------------
     287c ----- EXTINCTION NUAGES
     288           CALL histdef(nid_day,"kcld","Cld Ext proxy",
     289     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     290     .                "ave(X)",zsto,zout)
     291         endif
    211292c
    212293      ENDIF !lev_histday.GE.3

  • trunk/LMDZ.TITAN/libf/phytitan/ini_histins.h

    r107 r175  
    1 !
    2 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ini_histins.h,v 1.1.1.1 2004/05/19 12:53:08 lmdzadmin Exp $
    3 !
    41      IF (ok_instan) THEN
    52c
    … …  
    112109c
    113110         if (iflag_trac.eq.1) THEN
    114           if (microfi.eq.1) then
     111          if (microfi.ge.1) then
    115112           DO iq=1,nmicro
    116113         CALL histdef(nid_ins, tname(iq), ttext(iq), "n/m2",
    117114     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
    118      .                "ins(X)", zsto,zout)
     115     .                "inst(X)", zsto,zout)
    119116           ENDDO
    120117          endif
    … …  
    123120         CALL histdef(nid_ins, tname(iq), ttext(iq), "ppm",
    124121     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
    125      .                "ins(X)", zsto,zout)
     122     .                "inst(X)", zsto,zout)
    126123           ENDDO
    127124          endif
    … …  
    154151     .                32, "inst(X)", zsto,zout)
    155152c
     153c --------------
     154c ----- OPACITE BRUME
    156155         DO k=7,NSPECV,10
    157156           write(str1,'(i2.2)') k
    … …  
    161160         ENDDO
    162161c
     162         DO k=8,NSPECI,10
     163           write(str1,'(i2.2)') k
     164         CALL histdef(nid_ins,"thi"//str1,"Haze Opa IR",
     165     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     166     .                "ins(X)",zsto,zout)
     167         ENDDO
     168c
     169c --------------
     170c ----- EXTINCTION BRUME
    163171         DO k=7,NSPECV,10
    164172           write(str1,'(i2.2)') k
    … …  
    168176         ENDDO
    169177c
     178         DO k=8,NSPECI,10
     179           write(str1,'(i2.2)') k
     180         CALL histdef(nid_ins,"khi"//str1,"Haze ext IR ",
     181     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     182     .                "ins(X)",zsto,zout)
     183         ENDDO
     184c
     185c --------------
     186c ----- OPACITE GAZ
    170187         DO k=7,NSPECV,10
    171188           write(str1,'(i2.2)') k
    … …  
    175192         ENDDO
    176193c
     194         DO k=8,NSPECI,10
     195           write(str1,'(i2.2)') k
     196         CALL histdef(nid_ins,"tgi"//str1,"Haze Opa IR",
     197     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     198     .                "ins(X)",zsto,zout)
     199         ENDDO
     200c
     201c --------------
     202c ----- EXTINCTION GAZ
    177203         DO k=7,NSPECV,10
    178204           write(str1,'(i2.2)') k
    179205         CALL histdef(nid_ins,"kgv"//str1,"Haze ext Vis ",
    180206     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    181      .                "ins(X)",zsto,zout)
    182          ENDDO
    183 c
    184          DO k=8,NSPECI,10
    185            write(str1,'(i2.2)') k
    186          CALL histdef(nid_ins,"thi"//str1,"Haze Opa IR",
    187      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    188      .                "ins(X)",zsto,zout)
    189          ENDDO
    190 c
    191          DO k=8,NSPECI,10
    192            write(str1,'(i2.2)') k
    193          CALL histdef(nid_ins,"khi"//str1,"Haze ext IR ",
    194      .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    195      .                "ins(X)",zsto,zout)
    196          ENDDO
    197 c
    198          DO k=8,NSPECI,10
    199            write(str1,'(i2.2)') k
    200          CALL histdef(nid_ins,"tgi"//str1,"Haze Opa IR",
    201      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    202207     .                "ins(X)",zsto,zout)
    203208         ENDDO

  • trunk/LMDZ.TITAN/libf/phytitan/ini_histmth.h

    r106 r175  
    105105     .                "ave(X)", zsto1,zout)
    106106c
     107cccccccccccccccccc  Tracers
     108c
    107109         if (iflag_trac.eq.1) THEN
    108 c         if (microfi.eq.1) then
    109 c          DO iq=1,nmicro
    110 c        CALL histdef(nid_mth, tname(iq), ttext(iq), "n/m2",
    111 c    .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
    112 c    .                "ave(X)", zsto,zout)
    113 c          ENDDO
    114 c         endif
    115           if (nmicro.lt.nqmax) then
     110          if (microfi.ge.1) then
     111c           DO iq=1,nmicro
     112c             CALL histdef(nid_mth, tname(iq), ttext(iq), "n/m2",
     113c     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     114c     .                "ave(X)", zsto,zout)
     115c           ENDDO
     116             CALL histdef(nid_mth, "qaer","nb tot aer" , "n/m2",
     117     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     118     .                    "ave(X)", zsto,zout)
     119             CALL histdef(nid_mth, "qnoy","nb tot noy" , "n/m2",
     120     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     121     .                    "ave(X)", zsto,zout)
     122             CALL histdef(nid_mth, "qgl1","V tot gl1" , "m3/m2",
     123     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     124     .                    "ave(X)", zsto,zout)
     125             CALL histdef(nid_mth, "qgl2","V tot gl2" , "m3/m2",
     126     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     127     .                    "ave(X)", zsto,zout)
     128             CALL histdef(nid_mth, "qgl3","V tot gl3" , "m3/m2",
     129     .                    iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     130     .                    "ave(X)", zsto,zout)
     131c--------------
     132c ----- SATURATION ESP NUAGES
     133             if (clouds.eq.1) then
     134               CALL histdef(nid_mth,"ch4sat", "saturation CH4", "--",
     135     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     136     .                "ave(X)", zsto,zout)
     137               CALL histdef(nid_mth,"c2h6sat", "saturation C2H6", "--",
     138     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     139     .                "ave(X)", zsto,zout)
     140               CALL histdef(nid_mth,"c2h2sat", "saturation C2H2", "--",
     141     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     142     .                "ave(X)", zsto,zout)
     143c --------------
     144c ----- RESERVOIR DE SURFACE
     145               CALL histdef(nid_mth, "reserv", "Reservoir surface","m",
     146     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     147     .                "ave(X)", zsto,zout)
     148c --------------
     149c ----- PRECIPITATIONS (precipitations cumulatives)
     150               CALL histdef(nid_mth,"prech4","Precip CH4","m",
     151     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     152     .                "ave(X)", zsto,zout)
     153               CALL histdef(nid_mth,"prec2h6","Precip C2H6",
     154     .                "m",iim,jjmp1,nhori, 1,1,1, -99, 32,
     155     .                "ave(X)", zsto,zout)
     156               CALL histdef(nid_mth,"prec2h2","Precip C2H2",
     157     .                "m",iim,jjmp1,nhori, 1,1,1, -99, 32,
     158     .                "ave(X)", zsto,zout)
     159               CALL histdef(nid_mth,"prenoy","Precip NOY",
     160     .                "m",iim,jjmp1,nhori, 1,1,1, -99, 32,
     161     .                "ave(X)", zsto,zout)
     162               CALL histdef(nid_mth,"preaer","Precip AER",
     163     .                "m",iim,jjmp1,nhori, 1,1,1, -99, 32,
     164     .                "ave(X)", zsto,zout)
     165c --------------
     166c ----- FLUX GLACE
     167               CALL histdef(nid_mth,"flxgl1", "flux gl CH4",
     168     .              "kg/m2/s",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     169     .              "ave(X)", zsto,zout)
     170               CALL histdef(nid_mth,"flxgl2", "flux gl C2H6",
     171     .              "kg/m2/s",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     172     .              "ave(X)", zsto,zout)
     173               CALL histdef(nid_mth,"flxgl3", "flux gl C2H2",
     174     .              "kg/m2/s",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     175     .              "ave(X)", zsto,zout)
     176c --------------
     177c ----- RAYON DES GOUTTES
     178               CALL histdef(nid_mth,"rcldbar", "rayon moyen goutte",
     179     .                "m",iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     180     .                "ave(X)", zsto,zout)
     181             endif
     182          endif
     183c --------------
     184c ----- TRACEURS CHIMIQUES
     185          if (nmicro.lt.nqmax) then
    116186           DO iq=nmicro+1,nqmax
    117187         CALL histdef(nid_mth, tname(iq), ttext(iq), "ppm",
    … …  
    125195c    .                "ave(X)", zsto,zout)
    126196c          ENDDO
    127           endif
     197          endif
    128198         endif
    129199c
    … …  
    159229     .                32, "ave(X)", zsto1,zout)
    160230c
     231c --------------
     232c ----- OPACITE BRUME
    161233         DO k=7,NSPECV,10
    162234           write(str1,'(i2.2)') k
    … …  
    166238         ENDDO
    167239c
     240         DO k=8,NSPECI,10
     241           write(str1,'(i2.2)') k
     242         CALL histdef(nid_mth,"thi"//str1,"Haze Opa IR",
     243     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     244     .                "ave(X)",zsto1,zout)
     245         ENDDO
     246c
     247c --------------
     248c ----- EXTINCTION BRUME
    168249         DO k=7,NSPECV,10
    169250           write(str1,'(i2.2)') k
    … …  
    173254         ENDDO
    174255c
     256         DO k=8,NSPECI,10
     257           write(str1,'(i2.2)') k
     258         CALL histdef(nid_mth,"khi"//str1,"Haze ext IR ",
     259     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     260     .                "ave(X)",zsto1,zout)
     261         ENDDO
     262c
     263c --------------
     264c ----- OPACITE GAZ
    175265         DO k=7,NSPECV,10
    176266           write(str1,'(i2.2)') k
    177          CALL histdef(nid_mth,"tgv"//str1,"Haze Opa Vis",
    178      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    179      .                "ave(X)",zsto1,zout)
    180          ENDDO
    181 c
     267         CALL histdef(nid_mth,"tgv"//str1,"Gas Opa Vis",
     268     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     269     .                "ave(X)",zsto1,zout)
     270         ENDDO
     271c
     272         DO k=8,NSPECI,10
     273           write(str1,'(i2.2)') k
     274         CALL histdef(nid_mth,"tgi"//str1,"Haze Opa IR",
     275     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     276     .                "ave(X)",zsto1,zout)
     277         ENDDO
     278c
     279c --------------
     280c ----- EXTINCTION GAZ
    182281         DO k=7,NSPECV,10
    183282           write(str1,'(i2.2)') k
    184          CALL histdef(nid_mth,"kgv"//str1,"Haze ext Vis ",
     283         CALL histdef(nid_mth,"kgv"//str1,"Gas ext Vis ",
    185284     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    186285     .                "ave(X)",zsto1,zout)
    … …  
    189288         DO k=8,NSPECI,10
    190289           write(str1,'(i2.2)') k
    191          CALL histdef(nid_mth,"thi"//str1,"Haze Opa IR",
    192      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    193      .                "ave(X)",zsto1,zout)
    194          ENDDO
    195 c
    196          DO k=8,NSPECI,10
    197            write(str1,'(i2.2)') k
    198          CALL histdef(nid_mth,"khi"//str1,"Haze ext IR ",
    199      .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    200      .                "ave(X)",zsto1,zout)
    201          ENDDO
    202 c
    203          DO k=8,NSPECI,10
    204            write(str1,'(i2.2)') k
    205          CALL histdef(nid_mth,"tgi"//str1,"Haze Opa IR",
    206      .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    207      .                "ave(X)",zsto1,zout)
    208          ENDDO
    209 c
    210          DO k=8,NSPECI,10
    211            write(str1,'(i2.2)') k
    212          CALL histdef(nid_mth,"kgi"//str1,"Haze ext IR ",
    213      .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
    214      .                "ave(X)",zsto1,zout)
    215          ENDDO
     290         CALL histdef(nid_mth,"kgi"//str1,"Gas ext IR ",
     291     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     292     .                "ave(X)",zsto1,zout)
     293         ENDDO
     294c
     295c --------------
     296c ----- OPACITE NUAGES
     297         if (clouds.eq.1) then
     298           CALL histdef(nid_mth,"tcld","Cld Opa proxy",
     299     .                "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     300     .                "ave(X)",zsto,zout)
     301c
     302c --------------
     303c ----- EXTINCTION NUAGES
     304           CALL histdef(nid_mth,"kcld","Cld Ext proxy",
     305     .                "m-1",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     306     .                "ave(X)",zsto,zout)
     307         endif
     308c
     309c --------------
     310c ----- OCCURENCE NUAGES
     311           do k=1,12
     312             write(str1,'(i2.2)') k
     313             CALL histdef(nid_mth,"occcld"//str1,"occ cld",
     314     .       "--",iim,jjmp1,nhori,klev,1,klev,nvert,32,
     315     .       "ave(X)",zsto,zout)
     316           enddo
    216317c
    217318      ENDIF !lev_histmth.GE.3

  • trunk/LMDZ.TITAN/libf/phytitan/microtab.h

    r3 r175  
    1 *-------------------------------------------------------------------
    2 * INCLUDE microtab.h
    3 *
     1!-------------------------------------------------------------------
     2! INCLUDE microtab.h
     3!
    44         REAL wco,df_GP
    5          INTEGER nz,nrad,nztop
     5         INTEGER nz,nrad,imono,nztop,ntype
    66
    7          parameter(nz=llm,nrad=10,nztop=1)  !VERSION X
     7         parameter(nz=llm,ntype=5,nrad=10,imono=5,nztop=1)  !VERSION X
     8!        parameter(nz=llm,ntype=1,nrad=10,imono=5,nztop=1)  !VERSION X
    89
    910         parameter(wco=177.,df_GP=2.)      !FOR FRACTAL PARTCICLES
    10 c        parameter(wco=1.E+6,df_GP=3.)     !FOR SPHERE PARTICLES
     11!        parameter(wco=1.E+6,df_GP=3.)     !FOR SPHERE PARTICLES
    1112
    1213      real rf(nrad),df(nrad),zf,aknc
    1314      common/frac/rf,df,zf,aknc
    1415
    15 *********************************************************************
    16 * tcorrect, tx, microfi, cutoff: definis dans physiq.def (clesphys.h)
    17 *------------
     16!********************************************************************
     17! tcorrect, tx, microfi, cutoff: definis dans physiq.def (clesphys.h)
     18!------------
    1819         ! WARNING: tx=production rate
    1920         !          tcorrect is readjustment factor: =1 is continiuty
    2021         !                                           =X is q()*X
    21 *------------
    22 *(*1): si microfi=1, optcv et optci sont appeles a chaque appels de la
    23 *       physique  pour reactualiser les TAU's. De meme, pg2.F est
    24 *       active a chaque appel de la physique....
    25 *      si microfi=0., optcv et optci, ainsi que pg2, ne sont appele qu'une
    26 *       fois au debut, comme dans la version originale....
    27 *------------
    28 *       dans optci et optcv:
    29 *      si cutoff=1, brume coupee facon Pascal -> T ok au sol et dans la strato
    30 *                                             -> T tropopause mauvaise
    31 *                                             -> albedo ok
    32 *      si cutoff=2, brume coupee sous 100mbar -> T ok sol/tropopause/strato
    33 *                                             -> mais albedo mauvais
     22!------------
     23!(*1): si microfi=1, optcv et optci sont appeles a chaque appels de la
     24!       physique  pour reactualiser les TAU's. De meme, pg2.F est
     25!       active a chaque appel de la physique....
     26!      si microfi=0., optcv et optci, ainsi que pg2, ne sont appele qu'une
     27!       fois au debut, comme dans la version originale....
     28!------------
     29!       dans optci et optcv:
     30!      si cutoff=1, brume coupee facon Pascal -> T ok au sol et dans la strato
     31!                                             -> T tropopause mauvaise
     32!                                             -> albedo ok
     33!      si cutoff=2, brume coupee sous 100mbar -> T ok sol/tropopause/strato
     34!                                             -> mais albedo mauvais
    3435

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F

    r130 r175  
    3232     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
    3333
    34       COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
    35      & , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
     34      COMMON /CLOUD/
     35     &               RCLDI(NSPECI), XICLDI(NSPECI)
     36     &             , RCLDV(NSPECV), XICLDV(NSPECV)
     37     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
     38     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
    3639
    3740      COMMON /TAUS/   TAUHI(ngrid,NSPECI),TAUCI(ngrid,NSPECI),
    … …  
    4144
    4245      COMMON /TAUD/   TAUHID(ngrid,NLAYER,NSPECI)
     46     &               ,TAUCID(ngrid,NLAYER,NSPECI)
    4347     &               ,TAUGID(ngrid,NLAYER,NSPECI)
    4448     &               ,TAUHVD(ngrid,NLAYER,NSPECV)
     49     &               ,TAUCVD(ngrid,NLAYER,NSPECV)
    4550     &               ,TAUGVD(ngrid,NLAYER,NSPECV)
    4651
    … …  
    5055     &               ,WBARI(ngrid,NLAYER,NSPECI)
    5156     &               ,COSBI(ngrid,NLAYER,NSPECI)
     57     &               ,DTAUIP(ngrid,NLAYER,NSPECI)
     58     &               ,TAUIP(ngrid,NLEVEL,NSPECI)
     59     &               ,WBARIP(ngrid,NLAYER,NSPECI)
     60     &               ,COSBIP(ngrid,NLAYER,NSPECI)
    5261
    5362      COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
    5463     &                DWNI(NSPECI), WLNI(NSPECI)
     64
     65      REAL DTAUP(ngrid,NLAYER,NSPECI),DTAUPP(ngrid,NLAYER,NSPECI)
     66      COMMON /IRTAUS/ DTAUP,DTAUPP
    5567
    5668      COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
    … …  
    5870      COMMON /CONST/RGAS,RHOP,PI,SIGMA
    5971      COMMON /part/v,rayon,vrat,dr,dv
     72
     73c-----Rayons nuages et "composition" de la goutte
     74c     sur la grille ...
     75      integer ncount(ngrid,NLAYER)
     76      real    rmcbar(ngrid,NLAYER)
     77      real    xfbar(ngrid,NLAYER,4)
     78      COMMON/rnuabar/ncount,rmcbar,xfbar
    6079
    6180      DIMENSION PROD(NLEVEL)
    … …  
    7897      data iopti,iwarning,seulmtunpt/0,0,0/
    7998
    80       real   zqaer_1pt(NLAYER,nrad)
    81       real   TAUHID_1pt(NLAYER,NSPECI)
    82       real   TAUGID_1pt(NLAYER,NSPECI)
    83       real   TAUHI_1pt(NSPECI),TAUCI_1pt(NSPECI)
    84       real   TAUGI_1pt(NSPECI)
    85       real   DTAUI_1pt(NLAYER,NSPECI),TAUI_1pt(NLEVEL,NSPECI)
    86       real   WBARI_1pt(NLAYER,NSPECI)
    87       real   COSBI_1pt(NLAYER,NSPECI)
     99      real   zqaer_1pt(NLAYER,2*nrad)
     100#include "optci_1pt.h"
     101
    88102      character*100 dummy
    89103      real   dummy2,dummy3
    … …  
    128142c il faut quand meme qu on lise la look-up table de dim nrad=10
    129143c et si microfi=1, on doit avoir nmicro=nrad (dans microtab.h)
    130        if ((nmicro.ne.nrad).and.(microfi.eq.1)) then
    131           print*,"nmicro.ne.nrad",nmicro,nrad
    132           print*,"PROBLEME pour zqaer_1pt dans optci !!"
    133           stop
     144c
     145c Nouvelle verif pour nuages :
     146c La condition ci-dessus n'est plus realisable !
     147c nmicro comprend maintenant aussi des glaces
     148c Donc on teste juste que nmicro soit > 2*nrad (ou nrad si on ne fait pas de nuages)
     149       if (microfi.ge.1) then
     150         if ((clouds.eq.1).and.(nmicro.lt.2*nrad)) then
     151           print*,"OPTCI :"
     152           print*,"clouds = 1 MAIS nmicro < 2*nrad"
     153           print*,"Probleme pour zqaer_1pt dans optci."
     154           stop
     155         endif
     156         if ((clouds.eq.0).and.(nmicro.lt.nrad)) then
     157           print*,"OPTCI :"
     158           print*,"nmicro < nrad"
     159           print*,"Probleme pour zqaer_1pt dans optci."
     160           stop
     161         endif
    134162       endif
    135163
    136      
    137164      DO 420 K=1,NSPECI
    138165C LETS USE THE THOLIN OPTICAL CONSTANTS FOR THE HAZE.
    … …  
    141168          XIMGI(K)=TNI*FHIR
    142169C SET UP THE OPTICAL CONSTANTS FOR THE CLOUD
    143           RCLDI(K)=1.27
    144           XICLDI(K)=REFLIQ(WNOI(K))
     170          CALL LIQCH4(WLNI(K),TNR,TNI)
     171          RCLDI(K)=TNR
     172          XICLDI(K)=TNI
     173          CALL LIQC2H6(WLNI(K),TNR,TNI)
     174          RCLDI2(K)=TNR
     175          XICLDI2(K)=TNI
    145176 420  CONTINUE
    146177
    … …  
    197228c     endif   
    198229
    199         if (microfi.eq.1) then
    200            do iq=1,nrad
    201              do j=1,NLAYER
    202                  zqaer_1pt(j,iq)=qaer(ig,j,iq)
    203              enddo
    204            enddo
     230        if (microfi.ge.1) then
     231          do iq=1,2*nrad
     232c           si on ne fait pas de nuages on ne copie que les nrad premieres valeurs.
     233            if (clouds.eq.0.and.iq.gt.nrad) then
     234              zqaer_1pt(:,iq)=0.
     235            else
     236              do j=1,NLAYER
     237                zqaer_1pt(j,iq)=qaer(ig,j,iq)
     238              enddo
     239            endif
     240          enddo
    205241        else
    206242         if (ig.eq.1)  then
     243           zqaer_1pt = 0.
    207244c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
    208245c boucle sur nrad=10 (dans microtab.h)
    … …  
    217254c ici, les tableaux definissant la structure des aerosols sont
    218255c remplis: rf,df(nq),rayon(nq,)v(nq)......
    219        call rdf()
     256           call rdf()
    220257         endif
    221258        endif
    … …  
    229266       
    230267        if (seulmtunpt.eq.0) then
    231            call optci_1pt(zqaer_1pt,iopti,
    232      .            COSBI_1pt,DTAUI_1pt,TAUHI_1pt,TAUHID_1pt,TAUCI_1pt,
    233      .            TAUGI_1pt,TAUGID_1pt,WBARI_1pt,TAUI_1pt,IPRINT)
     268          call optci_1pt2(zqaer_1pt,rmcbar(ig,:),xfbar(ig,:,:),
     269     &                   iopti,IPRINT)
    234270           iopti = 1
    235271        endif
    … …  
    241277        endif
    242278       
    243         COSBI(ig,:,:)  = COSBI_1pt(:,:)
    244         WBARI(ig,:,:)  = WBARI_1pt(:,:)
    245         DTAUI(ig,:,:)  = DTAUI_1pt(:,:)
    246         TAUHI(ig,:)    = TAUHI_1pt(:)
    247         TAUCI(ig,:)    = TAUCI_1pt(:)
    248         TAUGI(ig,:)    = TAUGI_1pt(:)
    249         TAUI(ig,:,:)   = TAUI_1pt(:,:)
    250         TAUHID(ig,:,:) = TAUHID_1pt(:,:)
    251         TAUGID(ig,:,:) = TAUGID_1pt(:,:)
     279        COSBI(ig,:,:)  = COSBI_1pt(:,:)
     280        WBARI(ig,:,:)  = WBARI_1pt(:,:)
     281        DTAUI(ig,:,:)  = DTAUI_1pt(:,:)
     282        TAUI(ig,:,:)   = TAUI_1pt(:,:)
     283
     284        COSBIP(ig,:,:)  = COSBIP_1pt(:,:)
     285        WBARIP(ig,:,:)  = WBARIP_1pt(:,:)
     286        DTAUIP(ig,:,:)  = DTAUIP_1pt(:,:)
     287        TAUIP(ig,:,:)   = TAUIP_1pt(:,:)
     288
     289        TAUHI(ig,:)    = TAUHI_1pt(:)
     290        TAUCI(ig,:)    = TAUCI_1pt(:)
     291        TAUGI(ig,:)    = TAUGI_1pt(:)
     292
     293        TAUHID(ig,:,:) = TAUHID_1pt(:,:)
     294        TAUCID(ig,:,:) = TAUCID_1pt(:,:)
     295        TAUGID(ig,:,:) = TAUGID_1pt(:,:)
    252296
    253297c************************************************************************
    … …  
    256300c************************************************************************
    257301c************************************************************************
     302C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED
     303        DO 225 IG=1,klon
     304         DO 220 J=1,NLAYER
     305          DO 230 K=1,NSPECI
     306              DTAUP(IG,J,K)=DTAUI(IG,J,K)
     307              DTAUPP(IG,J,K)=DTAUIP(IG,J,K)
     308230       CONTINUE
     309220      CONTINUE
     310225     CONTINUE
    258311
    259312      print*, 'FIN OPTCI'

  • trunk/LMDZ.TITAN/libf/phytitan/optci_1pt.F

    r130 r175  
    1       SUBROUTINE optci_1pt(zqaer_1pt,iopti,
    2      .            COSBI_1pt,DTAUI_1pt,TAUHI_1pt,TAUHID_1pt,TAUCI_1pt,
    3      .            TAUGI_1pt,TAUGID_1pt,WBARI_1pt,TAUI_1pt,IPRINT)
     1      SUBROUTINE optci_1pt(zqaer_1pt,rcdb,xfrb,iopti,IPRINT)
    42      use dimphy
    53#include "dimensions.h"
    … …  
    1614C iopti: premier appel, on ne calcule qu'une fois les QM et QF
    1715* nrad dans microtab.h
    18       real   zqaer_1pt(NLAYER,nrad)
    19       real   TAUHID_1pt(NLAYER,NSPECI)
    20       real   TAUGID_1pt(NLAYER,NSPECI)
    21       real   TAUHI_1pt(NSPECI),TAUCI_1pt(NSPECI)
    22       real   TAUGI_1pt(NSPECI)
    23       real   DTAUI_1pt(NLAYER,NSPECI),TAUI_1pt(NLEVEL,NSPECI)
    24       real   WBARI_1pt(NLAYER,NSPECI)
    25       real   COSBI_1pt(NLAYER,NSPECI)
     16      real   zqaer_1pt(NLAYER,2*nrad)
     17#include "optci_1pt.h"
    2618c   ---------
    2719
    … …  
    3729     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
    3830
    39       COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
    40      & , RCLDI(NSPECI), XICLDI(NSPECI), RCLDV(NSPECV), XICLDV(NSPECV)
     31      COMMON /CLOUD/
     32     &               RCLDI(NSPECI), XICLDI(NSPECI)
     33     &             , RCLDV(NSPECV), XICLDV(NSPECV)
     34     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
     35     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
    4136
    4237      COMMON /SPECTI/ BWNI(NSPC1I), WNOI(NSPECI),
    … …  
    5752      REAL QM1(nrad,NSPECI),QM2(nrad,NSPECI)
    5853      REAL QM3(nrad,NSPECI),QM4(nrad,NSPECI)
     54      REAL QC1(nrad,NSPECI),QC2(nrad,NSPECI)
     55      REAL QC3(nrad,NSPECI),QC4(nrad,NSPECI)
    5956      real emu
    6057      REAL TAEROSM1(NSPECI),TAEROSCATM1(NSPECI),DELTAZM1(NSPECI)
    6158     
     59c ---- nuages     
     60      REAL TNUAGE,TNUAGESCAT
     61      REAL rcdb(nlayer),xfrb(nlayer,4)
     62
    6263      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
    6364
    … …  
    192193             if (iopti.eq.0) then
    193194
    194        CALL OPTFRAC(XMONO,10000./WNOI(K)
    195      &                         ,QEXT,QSCT,QABS,QBAR)
    196 
    197 
    198 c      CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
    199 c    &                ,XMONO,QSCT,QEXT,QABS,QBAR)
     195c       CALL OPTFRAC(XMONO,10000./WNOI(K)
     196c     &                         ,QEXT,QSCT,QABS,QBAR)
     197
     198
     199       CALL CFFFV11(1.e-2/WNOI(K),REALI(K),XIMGI(K),RF(inq),2.
     200     &                ,XMONO,QSCT,QEXT,QABS,QBAR)
    200201
    201202
    … …  
    213214
    214215               IF(TAEROS.LT.1.e-10) TAEROS=1.e-10                                             
    215 
    216216         ENDIF
    217217       ENDDO             ! FIN DE LA BOUCLE SUR nrad
    … …  
    279279
    280280C NEXT COMPUTE TAU CLOUD
    281       TAUCLD=0.0
    282       CBARC =0.0
    283        IF ( XNCLD(J) .GT. 0..and .taufac.gt.0.) THEN
    284          print*,'Appel a xmie avec radcld=',radcld(j)
    285                 CALL XMIE(RADCLD(J),RCLDI(K),XICLDI(K),
    286      &                         QEXTC,QSCTC,QABSC,CBARC,WNOI(K))
    287                 TAUCLD=QEXTC*XNCLD(J)
     281      IF (clouds.eq.0) THEN
     282        TNUAGE=0.
     283        TNUAGESCAT=0.
     284        CNBAR=0.
     285      ELSE
     286        TNUAGE=0.
     287        TNUAGESCAT=0.
     288        CNBAR=0.
     289        QEXTC=0.
     290        QSCTC=0.
     291        QABSC=0.
     292        CBARC=0.
     293
     294        DO inq=1,nrad         !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
     295          QC1(inq,K)=0.
     296          QC2(inq,K)=0.
     297          QC3(inq,K)=0.
     298          QC4(inq,K)=0.
     299        ENDDO
     300
     301** OPTICAL CONSTANT : MIXING RULES
     302
     303        IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
     304
     305          XNR=xfrb(nlayer+1-J,1) *REALI(K)
     306     &    +xfrb(nlayer+1-J,2) *RCLDI(K)
     307     &    +xfrb(nlayer+1-J,3) *RCLDI2(K)
     308     &    +xfrb(nlayer+1-J,4) *RCLDI2(K)
     309
     310          XNI=xfrb(nlayer+1-J,1) *XIMGI(K)
     311     &    +xfrb(nlayer+1-J,2) *XICLDI(K)
     312     &    +xfrb(nlayer+1-J,3) *XICLDI2(K)
     313     &    +xfrb(nlayer+1-J,4) *XICLDI2(K)
     314
     315** OPTICAL CONSTANT : LIQUID DROP = THOLIN
     316
     317          IF(xfrb(nlayer+1-J,1).ge.0.1) THEN
     318            XNI=XIMGI(K)
     319            XNR=REALI(K)
     320          ENDIF
     321
     322          IF (XNI.gt.1.e-10  .and. XNR.gt.1.00) THEN
     323            CALL CMIE(1.E-2/WNOI(K),XNR,XNI,
     324     &      rcdb(nlayer+1-J),
     325     &      QEXTC,QSCTC,QABSC,CBARC)
     326          ELSE
     327            PRINT*,' WARNING XNR/XNI in optci: ',XNR,XNI
     328            QEXTC=0.
     329            QSCTC=0.
     330            QABSC=0.
     331            CBARC=0.
     332          ENDIF
     333        ELSE
     334          QEXTC=0.
     335          QSCTC=0.
     336          QABSC=0.
     337          CBARC=0.
    288338        ENDIF
     339
     340        DO inq=1,nrad         !BOUCLE SUR LES NQMX TAILLE D"AEROSOLS
     341          QC1(inq,K)=QEXTC/xnuf
     342          QC2(inq,K)=QSCTC/xnuf
     343          QC3(inq,K)=QABSC/xnuf
     344          QC4(inq,K)=CBARC
     345          TNUAGE=QC1(inq,K)*zqaer_1pt(nlayer+1-J,inq+nrad)*1.e-4
     346     &          +TNUAGE
     347          TNUAGESCAT=QC2(inq,K)*zqaer_1pt(nlayer+1-J,inq+nrad)*1.e-4
     348     &              +TNUAGESCAT
     349          CNBAR=QC4(inq,K)*QC2(inq,K)
     350     &         *zqaer_1pt(nlayer+1-J,inq+nrad)*1.e-4+CNBAR
     351        ENDDO
     352
     353        IF(TNUAGESCAT.EQ.0.) THEN
     354             CNBAR=0.
     355        ELSE
     356             CNBAR=CNBAR/TNUAGESCAT
     357        ENDIF
     358      ENDIF    ! Cond CLD
    289359
    290360
    … …  
    329399C
    330400
    331       IF (TAEROS + TAUCLD .GT. 0.) THEN
    332          COSBI_1pt(J,K)=(CBAR*TAEROSCAT + CBARC*TAUCLD )
    333      &                     /(TAEROSCAT + TAUCLD)
     401      DTAUI_1pt(J,K)=TAUGAS+TAEROS+TNUAGE
     402      DTAUIP_1pt(J,K)=TAUGAS+TAEROS
     403
     404      TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
     405      TAUHID_1pt(J,K)=TAUHI_1pt(K)
     406
     407      TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
     408      TAUGID_1pt(J,K)=TAUGI_1pt(K)
     409
     410      TAUCI_1pt(K)=TAUCI_1pt(K) + TNUAGE
     411      TAUCID_1pt(J,K)=TAUCI_1pt(K)
     412 
     413C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
     414
     415      TLIMIT=1.E-16
     416
     417
     418      IF (TAEROSCAT + TNUAGESCAT .GT. 0.) THEN
     419         COSBI_1pt(J,K)=(CBAR*TAEROSCAT + CNBAR*TNUAGESCAT )
     420     &                     /(TAEROSCAT + TNUAGESCAT)
    334421      ELSE
    335422         COSBI_1pt(J,K)=0.0
    336423      ENDIF
    337424
    338       DTAUI_1pt(J,K)=TAUGAS+TAEROS+TAUCLD
    339 
    340       TAUHI_1pt(K)=TAUHI_1pt(K) + TAEROS
    341       TAUHID_1pt(J,K)=TAUHI_1pt(K)
    342 
    343       TAUGI_1pt(K)=TAUGI_1pt(K) + TAUGAS
    344       TAUGID_1pt(J,K)=TAUGI_1pt(K)
    345 
    346       TAUCI_1pt(K)=TAUCI_1pt(K) + TAUCLD
    347  
    348 c     if (ig.eq.1) then
    349 c     if (k.eq.NSPECI/2.or.k.eq.1.or.k.eq.NSPECI) then
    350 c      print*,'@IR',K,J,DTAUI_1pt(J,K),TAUGAS,TAEROS,TAUCLD         
    351 c     endif
    352 c     endif
    353 
    354 
    355 C ??FLAG? SERIOUS PROBLEM WITH THE CODE HERE!
    356 
    357       TLIMIT=1.E-16
    358 
    359       IF (DTAUI_1pt(J,K) .GT. TLIMIT) THEN
    360 
    361 c***************** ECHANGE
    362 c        WBARI(J,K)=(QSCT*XNUMB(J) + QSCTC*XNCLD(J))
    363 c****************
    364 CFC        WBARI_1pt(J,K)=(TAEROSCAT + QSCTC*XNCLD(J))
    365 c****************
    366            WBARI_1pt(J,K)=TAEROSCAT
    367      &   /DTAUI_1pt(J,K)
    368 c        if(iwarning.eq.0)
    369 c     s   print*,'WARNING!!! dans optci xncld non initialise'
    370 c        iwarning=1
    371 
     425      IF (TAEROSCAT  .GT. 0.) THEN
     426         COSBIP_1pt(J,K)=(CBAR*TAEROSCAT)
     427     &                     /(TAEROSCAT)
    372428      ELSE
    373 
     429         COSBIP_1pt(J,K)=0.0
     430      ENDIF
     431
     432*---------
     433
     434      IF (DTAUI_1pt(J,K) .GT.  TLIMIT) THEN
     435          WBARI_1pt(J,K)=(TAEROSCAT+TNUAGESCAT) /DTAUI_1pt(J,K)
     436      ELSE
    374437         WBARI_1pt(J,K)=0.0
    375 c        PRINT*,'gaz ',taugas,'aerosols',taeros,'nuages',taucld
    376 c        WRITE (6,999) J,K,DTAUI_1pt(J,K)
    377  999     FORMAT ('  WARNING TAU MINIMUM AT J,K,DTAU:',2I3,1PE10.3)
    378438         DTAUI_1pt(J,K)=TLIMIT
    379 
    380439      ENDIF
     440
     441      IF (DTAUIP_1pt(J,K) .GT.  TLIMIT) THEN
     442          WBARIP_1pt(J,K)=(TAEROSCAT) /DTAUIP_1pt(J,K)
     443      ELSE
     444         WBARIP_1pt(J,K)=0.0
     445         DTAUIP_1pt(J,K)=TLIMIT
     446      ENDIF
     447
    381448
    382449c     IF (IPRINT .GT. 9)
    … …  
    397464        DO 119 K=1,NSPECI
    398465           TAUI_1pt(1,K)=0.0
     466           TAUIP_1pt(1,K)=0.0
    399467        DO 119 J=1,NLAYER
    400468           TAUI_1pt(J+1,K)=TAUI_1pt(J,K)+DTAUI_1pt(J,K)
     469           TAUIP_1pt(J+1,K)=TAUIP_1pt(J,K)+DTAUIP_1pt(J,K)
    401470 119    CONTINUE
    402471

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F

    r119 r175  
    3131     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
    3232
    33       COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
    34      &             , RCLDI(NSPECI), XICLDI(NSPECI)
     33      COMMON /CLOUD/
     34     &               RCLDI(NSPECI), XICLDI(NSPECI)
    3535     &             , RCLDV(NSPECV), XICLDV(NSPECV)
     36     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
     37     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
    3638
    3739      COMMON /TAUS/   TAUHI(ngrid,NSPECI), TAUCI(ngrid,NSPECI)
    … …  
    4143
    4244      COMMON /TAUD/   TAUHID(ngrid,NLAYER,NSPECI)
     45     &               ,TAUCID(ngrid,NLAYER,NSPECI)
    4346     &               ,TAUGID(ngrid,NLAYER,NSPECI)
    4447     &               ,TAUHVD(ngrid,NLAYER,NSPECV)
     48     &               ,TAUCVD(ngrid,NLAYER,NSPECV)
    4549     &               ,TAUGVD(ngrid,NLAYER,NSPECV)
    4650
    … …  
    4953     &               ,WBARV(ngrid,NLAYER,NSPECV,4)
    5054     &               ,COSBV(ngrid,NLAYER,NSPECV,4)
     55     &               ,DTAUVP(ngrid,NLAYER,NSPECV,4)
     56     &               ,TAUVP(ngrid,NLEVEL,NSPECV,4)
     57     &               ,WBARVP(ngrid,NLAYER,NSPECV,4)
     58     &               ,COSBVP(ngrid,NLAYER,NSPECV,4)
    5159
    5260      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
    … …  
    5765      COMMON /CONST/ RGAS,RHOP,PI,SIGMA
    5866      COMMON /part/ v(nrad),rayon(nrad),vrat,dr(nrad),dv(nrad)
     67
     68c-----Rayons nuages et "composition" de la goutte
     69c     sur la grille ...
     70      integer ncount(ngrid,NLAYER)
     71      real    rmcbar(ngrid,NLAYER)
     72      real    xfbar(ngrid,NLAYER,4)
     73      COMMON/rnuabar/ncount,rmcbar,xfbar
    5974
    6075      REAL xv1(klev,NSPECV)
    … …  
    7489      data ioptv,iwarning,seulmtunpt/0,0,0/
    7590
    76       real   zqaer_1pt(NLAYER,nrad)
    77       real   TAUHVD_1pt(NLAYER,NSPECV)
    78       real   TAUGVD_1pt(NLAYER,NSPECV)
    79       real   TAUHV_1pt(NSPECV),TAUCV_1pt(NSPECV)
    80       real   TAURV_1pt(NSPECV),TAUGV_1pt(NSPECV)
    81       real   DTAUV_1pt(NLAYER,NSPECV,4),TAUV_1pt(NLEVEL,NSPECV,4)
    82       real   WBARV_1pt(NLAYER,NSPECV,4)
    83       real   COSBV_1pt(NLAYER,NSPECV,4)
     91      real   zqaer_1pt(NLAYER,2*nrad)
     92#include "optcv_1pt.h"
     93
    8494      character*100 dummy
    8595      real   dummy2,dummy3
    … …  
    8797C*
    8898C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
    89 C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
     99C IT CALCULATES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VIS
    90100C LAYER: WBAR, DTAU, COSBAR
    91101C LEVEL: TAU
    … …  
    95105       print*,'ATTENTION, TAU UNIFORME DANS OPTCV'
    96106
    97 c      do nng=2,klon
    98 c       do i=1,klev           
    99 c        do j=1,nqtot
    100 c          sum=sum+qaer(nng,i,j)*rayon(j)**3.*1.3333*3.1415*1000.
    101 c        enddo
    102 c       enddo
    103 c       enddo
    104 c       print*,sum/(klon-1),'SOMME COLONNE/OPTCV'
    105 
    106              
    107107C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    108108c INITIALISATIONS UNE SEULE FOIS
    … …  
    114114c il faut quand meme qu'on lise la look-up table de dim nrad=10
    115115c et si microfi=1, on doit avoir nmicro=nrad (dans microtab.h)
    116        if ((nmicro.ne.nrad).and.(microfi.eq.1)) then
    117           print*,"nmicro.ne.nrad",nmicro,nrad
    118           print*,"PROBLEME pour zqaer_1pt dans optcv !!"
    119           stop
     116c
     117c Nouvelle verif pour nuages :
     118c La condition ci-dessus n'est plus realisable !
     119c nmicro comprend maintenant aussi des glaces
     120c Donc on teste juste que nmicro soit > 2*nrad (ou nrad si on ne fait pas de nuages)
     121       if (microfi.ge.1) then
     122         if ((clouds.eq.1).and.(nmicro.lt.2*nrad)) then
     123           print*,"OPTCV :"
     124           print*,"clouds = 1 MAIS nmicro < 2*nrad"
     125           print*,"Probleme pour zqaer_1pt dans optcv."
     126           stop
     127         endif
     128         if ((clouds.eq.0).and.(nmicro.lt.nrad)) then
     129           print*,"OPTCV :"
     130           print*,"nmicro < nrad"
     131           print*,"Probleme pour zqaer_1pt dans optcv."
     132           stop
     133         endif
    120134       endif
    121 
    122135     
    123136      DO 130 K=1,NSPECV
    … …  
    130143C      XIMGV(K)=FITEDN(WLNV(K))
    131144C THE CLOUD IS CLEAR IN THE VISIBLE
    132       RCLDV(K)=1.27
    133       XICLDV(K)=1.E-7
     145      CALL LIQCH4(WLNV(K),TNR,TNI)
     146      RCLDV(K)=TNR
     147      XICLDV(K)=TNI
     148      CALL LIQC2H6(WLNV(K),TNR,TNI)
     149      RCLDV2(K)=TNR
     150      XICLDV2(K)=TNI
    134151 130  CONTINUE
    135152C
    … …  
    150167      DO 101 ig=1,klon       !c! BOUCLE SUR GRILLE HORIZONTALE
    151168
    152         if (microfi.eq.1) then
    153            do iq=1,nrad
    154               do j=1,NLAYER
    155                  zqaer_1pt(j,iq)=qaer(ig,j,iq)
    156               enddo
    157            enddo
     169        if (microfi.ge.1) then
     170           do iq=1,2*nrad
     171             if (clouds.eq.0.and.iq.gt.nrad) then
     172                zqaer_1pt(:,iq)=0.
     173             else
     174               do j=1,NLAYER
     175                  zqaer_1pt(j,iq)=qaer(ig,j,iq)
     176               enddo
     177             endif
     178           enddo
    158179        else
    159180         if (ig.eq.1)  then
    … …  
    181202c       if ((microfi.eq.0).or.(ig.eq.klon/2)) iout=1
    182203        if (seulmtunpt.eq.0) then
    183            call optcv_1pt(zqaer_1pt,ioptv,
    184      .            COSBV_1pt,DTAUV_1pt,TAUHV_1pt,TAUHVD_1pt,TAUCV_1pt,
    185      .       TAURV_1pt,TAUGV_1pt,TAUGVD_1pt,WBARV_1pt,TAUV_1pt,iout)
     204          call optcv_1pt2(zqaer_1pt,rmcbar(ig,:),xfbar(ig,:,:),
     205     &                   ioptv,IPRINT)
    186206           ioptv = 1
    187207        endif
    188208
    189209c Pas de microphysique, ni de composition variable: un seul passage
    190 c dans optci_1pt.
     210c dans optcv_1pt.
    191211        if ((microfi.eq.0).and.(ylellouch)) then
    192212           seulmtunpt = 1
    193213        endif
    194214       
    195         COSBV(ig,:,:,:)= COSBV_1pt(:,:,:)
    196         WBARV(ig,:,:,:)= WBARV_1pt(:,:,:)
    197         DTAUV(ig,:,:,:)= DTAUV_1pt(:,:,:)
    198         TAUHV(ig,:)    = TAUHV_1pt(:)
    199         TAUCV(ig,:)    = TAUCV_1pt(:)
    200         TAURV(ig,:)    = TAURV_1pt(:)
    201         TAUGV(ig,:)    = TAUGV_1pt(:)
    202         TAUV(ig,:,:,:) = TAUV_1pt(:,:,:)
    203         TAUHVD(ig,:,:) = TAUHVD_1pt(:,:)
    204         TAUGVD(ig,:,:) = TAUGVD_1pt(:,:)
     215        COSBV(ig,:,:,:)= COSBV_1pt(:,:,:)
     216        WBARV(ig,:,:,:)= WBARV_1pt(:,:,:)
     217        DTAUV(ig,:,:,:)= DTAUV_1pt(:,:,:)
     218        TAUV(ig,:,:,:) = TAUV_1pt(:,:,:)
     219
     220        COSBVP(ig,:,:,:)= COSBVP_1pt(:,:,:)
     221        WBARVP(ig,:,:,:)= WBARVP_1pt(:,:,:)
     222        DTAUVP(ig,:,:,:)= DTAUVP_1pt(:,:,:)
     223        TAUVP(ig,:,:,:) = TAUVP_1pt(:,:,:)
     224
     225        TAUHV(ig,:)    = TAUHV_1pt(:)
     226        TAUCV(ig,:)    = TAUCV_1pt(:)
     227        TAURV(ig,:)    = TAURV_1pt(:)
     228        TAUGV(ig,:)    = TAUGV_1pt(:)
     229
     230        TAUHVD(ig,:,:) = TAUHVD_1pt(:,:)
     231        TAUCVD(ig,:,:) = TAUCVD_1pt(:,:)
     232        TAUGVD(ig,:,:) = TAUGVD_1pt(:,:)
    205233
    206234 101  CONTINUE

  • trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt.F

    r102 r175  
    1       SUBROUTINE optcv_1pt(zqaer_1pt,ioptv,
    2      .            COSBV_1pt,DTAUV_1pt,TAUHV_1pt,TAUHVD_1pt,TAUCV_1pt,
    3      .       TAURV_1pt,TAUGV_1pt,TAUGVD_1pt,WBARV_1pt,TAUV_1pt,IPRINT)
     1      SUBROUTINE optcv_1pt(zqaer_1pt,rcdb,xfrb,ioptv,IPRINT)
    42
    53
    … …  
    1715C ioptv: premier appel, on ne calcule qu'une fois les QM et QF
    1816* nrad dans microtab.h
    19       real   zqaer_1pt(NLAYER,nrad)
    20       real   TAUHVD_1pt(NLAYER,NSPECV)
    21       real   TAUGVD_1pt(NLAYER,NSPECV)
    22       real   TAUHV_1pt(NSPECV),TAUCV_1pt(NSPECV)
    23       real   TAURV_1pt(NSPECV),TAUGV_1pt(NSPECV)
    24       real   DTAUV_1pt(NLAYER,NSPECV,4),TAUV_1pt(NLEVEL,NSPECV,4)
    25       real   WBARV_1pt(NLAYER,NSPECV,4)
    26       real   COSBV_1pt(NLAYER,NSPECV,4)
     17      real   zqaer_1pt(NLAYER,2*nrad)
     18#include "optcv_1pt.h"
    2719c   ---------
    2820
    … …  
    3830     & , REALI(NSPECI), XIMGI(NSPECI), REALV(NSPECV), XIMGV(NSPECV)
    3931
    40       COMMON /CLOUD/ RADCLD(NLAYER), XNCLD(NLAYER)
    41      &             , RCLDI(NSPECI), XICLDI(NSPECI)
     32      COMMON /CLOUD/
     33     &               RCLDI(NSPECI), XICLDI(NSPECI)
    4234     &             , RCLDV(NSPECV), XICLDV(NSPECV)
     35     &             , RCLDI2(NSPECI), XICLDI2(NSPECI)
     36     &             , RCLDV2(NSPECV), XICLDV2(NSPECV)
    4337
    4438      COMMON /SPECTV/ BWNV(NSPC1V),WNOV(NSPECV)
    … …  
    5549      REAL QM1(nrad,NSPECV),QM2(nrad,NSPECV)
    5650      REAL QM3(nrad,NSPECV),QM4(nrad,NSPECV)
     51      REAL QC1(nrad,NSPECV),QC2(nrad,NSPECV)
     52      REAL QC3(nrad,NSPECV),QC4(nrad,NSPECV)
     53
     54c---- NUAGES
     55      real TNUABS,TNUSCAT
     56      real   rcdb(NLAYER)
     57      real     xfrb(NLAYER,4)
    5758
    5859      save qf1,qf2,qf3,qf4,qm1,qm2,qm3,qm4
    5960 
     61      integer ilat,jalt
     62      common/toto/ilat,jalt
     63
    6064C*
    6165C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISIBLE
    … …  
    7781
    7882      DO 100 J=1,NLAYER         !a! BOUCLE SUR L"ALTITUDE
    79 
     83        jalt=j
    8084C #1:                   HAZE
    8185c---------------------------
    … …  
    149153            if(ioptv.eq.0.and.J.eq.1) then
    150154
    151         CALL OPTFRAC(XMONO,10000./WNOV(K)
    152      &                        ,QEXT,QSCT,QABS,QBAR)
    153 
    154 c       CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
    155 c    &   ,XMONO,QSCT,QEXT,QABS,QBAR)
     155c        CALL OPTFRAC(XMONO,10000./WNOV(K)
     156c     &                        ,QEXT,QSCT,QABS,QBAR)
     157
     158        CALL CFFFV11(1.e-2/WNOV(K),REALV(K),XIMGV(K),RF(inq),2.
     159     &   ,XMONO,QSCT,QEXT,QABS,QBAR)
    156160
    157161
    … …  
    215219       IF (TAEROSCAT.le.0.) CBAR=0.
    216220
    217       if (IPRINT.eq.1) then
    218       if (k.eq.NSPECV/2) then   
    219        print*,'@VI',K,J,TAEROS,TAEROSCAT,CBAR
    220        print*,'@  ',K,J,QF1(1,K),QF2(1,K),zqaer_1pt(NLAYER+1-J,1)
    221        print*,'@  ',K,J,QF1(3,K),QF2(3,K),zqaer_1pt(NLAYER+1-J,3)
    222        print*,'@  ',K,J,QF1(5,K),QF2(5,K),zqaer_1pt(NLAYER+1-J,5)
    223        print*,'@  ',K,J,QF1(7,K),QF2(7,K),zqaer_1pt(NLAYER+1-J,7)
    224        print*,'@  ',K,J,QF1(9,K),QF2(9,K),zqaer_1pt(NLAYER+1-J,9)
    225        print*
    226       endif
    227       endif
    228 
    229 
     221c      if (IPRINT.eq.1) then
     222c      if (k.eq.NSPECV/2) then   
     223c       write(*,1699) '@VI',K,J,TAEROS,TAEROSCAT,CBAR
     224c       write(*,1699) '@  ',K,J,QF1(1,K),QF2(1,K),zqaer_1pt(NLAYER+1-J,1)
     225c       write(*,1699) '@  ',K,J,QF1(3,K),QF2(3,K),zqaer_1pt(NLAYER+1-J,3)
     226c       write(*,1699) '@  ',K,J,QF1(5,K),QF2(5,K),zqaer_1pt(NLAYER+1-J,5)
     227c       write(*,1699) '@  ',K,J,QF1(7,K),QF2(7,K),zqaer_1pt(NLAYER+1-J,7)
     228c       write(*,1699) '@  ',K,J,QF1(9,K),QF2(9,K),zqaer_1pt(NLAYER+1-J,9)
     229c       print*
     230c      endif
     231c      endif
     232
     2331699  FORMAT(a3,2I3,3(ES15.7,1X))
    230234
    231235c*********** EN TRAVAUX ***************************
    … …  
    280284C NEXT COMPUTE TAU CLOUD
    281285
    282       TAUCLD=0.0
    283       CBARC =0.0
    284       QEXTC =0.0
    285       QSCTC =0.0
    286 c             XNCLD(J)=0.
    287       IF ( XNCLD(J) .GT. 0. .and .taufac.gt.0.) THEN
    288                 CALL XMIE(RADCLD(J),RCLDV(K),XICLDV(K),
    289      &                         QEXTC,QSCTC,QABSC,CBARC,WNOV(K))
    290                 TAUCLD=QEXTC*XNCLD(J)         
    291       ENDIF
    292 C
    293       TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
    294       TAUGVD_1pt(J,K)=TAURV_1pt(K)
    295 
    296       TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS          ! INTEGRATED Quant.
    297       TAUHVD_1pt(J,K)=TAUHV_1pt(K)
    298 
    299       TAUCV_1pt(K)=TAUCV_1pt(K)+TAUCLD
     286       IF (clouds.eq.0) THEN
     287         CNBAR=0.
     288         TNUSCAT=0.
     289         TNUABS=0.
     290         TBNUABS=0.
     291       ELSE
     292         CNBAR=0.
     293         TNUSCAT=0.
     294         TNUABS=0.
     295         TBNUABS=0.
     296         QEXTC=0.
     297         QSCTC=0.
     298         QABSC=0.
     299         CBARC=0.
     300
     301         do inq=1,nrad
     302           QC1(INQ,k)=0.
     303           QC2(INQ,k)=0.
     304           QC3(INQ,k)=0.
     305           QC4(INQ,k)=0.
     306         enddo
     307
     308         IF (rcdb(nlayer+1-J).gt.1.1e-10) THEN
     309
     310** OPTICAL CONSTANT : MIXING RULES
     311
     312           XNR=xfrb(nlayer+1-J,1)*REALV(K)                   !
     313     &        +xfrb(nlayer+1-J,2)*RCLDV(K)                   !
     314     &        +xfrb(nlayer+1-J,3)*RCLDV2(K)                  !
     315     &        +xfrb(nlayer+1-J,4)*RCLDV2(K)                  !
     316
     317           XNI=xfrb(nlayer+1-J,1)*XIMGV(K)
     318     &        +xfrb(nlayer+1-J,2)*XICLDV(K)
     319     &        +xfrb(nlayer+1-J,3)*XICLDV2(K)
     320     &        +xfrb(nlayer+1-J,4)*XICLDV2(K)
     321
     322** OPTICAL CONSTANT : LIQUID DROP = THOLIN
     323           IF(xfrb(nlayer+1-J,1).ge.0.01) THEN
     324             XNI=XIMGV(K)
     325             XNR=REALV(K)
     326           ENDIF
     327
     328           IF (XNI.gt.1.e-10  .and. XNR.gt.1.00) THEN
     329             CALL CMIE(1.E-2/WNOV(K),XNR,XNI,
     330     &       rcdb(nlayer+1-J),
     331     &       QEXTC,QSCTC,QABSC,CBARC)
     332           ELSE
     333             PRINT*,' WARNING XNR/XNI in optcv: ',XNR,XNI
     334             QEXTC=0.
     335             QSCTC=0.
     336             QABSC=0.
     337             CBARC=0.
     338             STOP
     339           ENDIF
     340         ELSE
     341           QEXTC=0.
     342           QSCTC=0.
     343           QABSC=0.
     344           CBARC=0.
     345         ENDIF
     346
     347         DO inq=1,nrad
     348
     349           QC1(INQ,k)=QEXTC/xnuf
     350           QC2(INQ,k)=QSCTC/xnuf
     351           QC3(INQ,k)=QABSC/xnuf
     352           QC4(INQ,k)=CBARC
     353
     354           TNUABS=QC1(inq,K)*zqaer_1pt(NLAYER+1-J,inq+nrad)*1.e-4
     355     &           +TNUABS
     356
     357           TNUSCAT=QC2(inq,K)*zqaer_1pt(NLAYER+1-J,inq+nrad)*1.e-4
     358     &            +TNUSCAT
     359
     360           CNBAR=QC4(inq,K)*QC2(inq,K)*
     361     &           zqaer_1pt(NLAYER+1-J,inq+nrad)*1.e-4 + CNBAR
     362
     363         ENDDO
     364
     365         IF(TNUSCAT.EQ.0) CNBAR=0.
     366         IF(TNUSCAT.NE.0.) CNBAR=CNBAR/TNUSCAT
     367
     368
     369       ENDIF  ! Cond. CLD
     370
     371       TAURV_1pt(K)=TAURV_1pt(K)+TAURAY
     372       TAUGVD_1pt(J,K)=TAURV_1pt(K)
     373
     374       TAUHV_1pt(K)=TAUHV_1pt(K)+TAEROS ! INTEGRATED Quant.
     375       TAUHVD_1pt(J,K)=TAUHV_1pt(K)
     376
     377       TAUCV_1pt(K)=TAUCV_1pt(K)+TNUABS
     378       TAUCVD_1pt(J,K)=TAUCV_1pt(K)
     379
    300380
    301381C #4:                  TAUGAS
    … …  
    311391
    312392
    313 C COMPUTE THE AVERAGE COSBAR AND WBAR
    314 C&&
    315 
    316 c     CBAR=MIN(1.0,1.05*CBAR)       ! THE HAZE FORWARD SCATTERING 5%(WHY?)
    317       COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CBARC*TAUCLD)
    318      &  /(TAEROSCAT+TAUCLD+TAURAY)     !CBAR_RAY=0.
    319 c        print*,'CBV',J,K,NT,CBAR,TAEROSCAT,CBARC,TAUCLD       
    320 
    321       DTAUV_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY+TAUCLD       !TOTAL TAU_EXT
    322       TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K)         !TAU_ABS_METH INTEG.
    323 
    324 C WE LET W RAYLEIGH BE .999 OR W=1 WHEN ONLY RAYLEIGH=PROBLEM FOR TRID
    325 c WE HAVE ASSUMED ABOVE THAT COSBAR FOR RAYLEIGH IS ZERO.
    326       if (IPRINT.eq.1) then
    327       if (k.eq.NSPECV/2) then   
    328        print*,'@VI',K,J,DTAUV_1pt(J,K,1),TAUGAS,TAEROS,TAUCLD
    329       endif
    330       endif
    331 
    332 
    333 c***************** ECHANGE
    334 c     WBARV(J,K,NT)=(QSCT*XNUMB(J)+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
    335 c****************
    336       WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + QSCTC*XNCLD(J) )
    337 c     WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
    338      &            /(TAUGAS+TAEROS+TAURAY+TAUCLD)
    339 c****************
    340       IF((TAEROS+TAUCLD+TAURAY+TAUCLD).le.0.) WBARV_1pt(J,K,NT)=0.
    341       IF((TAEROS+TAUCLD+TAURAY).le.0.) COSBV_1pt(J,K,NT)=0.
    342 
    343 c     print*,'WBV',J,K,NT,TAEROSCAT,TAURAY,QSCTC*XNCLD(J)
    344 c     print*,'WBV',J,K,NT,TAEROS,TAUGAS,TAURAY,TAUCLD
    345 c     print*,Z(j),J,K,NT,TAUV(1,j,K,NT),WBARV(1,j,K,NT),COSBV(1,j,K,NT)
     393*  COSBV ET COSBVP
     394*-----------------
     395
     396      IF(TAEROSCAT+TNUSCAT+TAURAY .ne. 0.) THEN
     397         COSBV_1pt(J,K,NT)=(CBAR*TAEROSCAT + CNBAR*TNUSCAT)
     398     &  /(TAEROSCAT+TNUSCAT+TAURAY) !CBAR_RAY=0.
     399      ELSE
     400         COSBV_1pt(J,K,NT)=0.
     401      ENDIF
     402
     403      IF(TAEROSCAT+TAURAY .ne. 0.) THEN
     404         COSBVP_1pt(J,K,NT)=(CBAR*TAEROSCAT)
     405     &  /(TAEROSCAT+TAURAY) !CBAR_RAY=0.
     406      ELSE
     407         COSBVP_1pt(J,K,NT)=0.
     408      ENDIF
     409
     410*  DTAUV ET DTAUVP
     411*-----------------
     412
     413      DTAUV_1pt(J,K,NT) =TAUGAS+TAEROS+TAURAY+TNUABS !TAU_ABS_METH
     414      DTAUVP_1pt(J,K,NT)=TAUGAS+TAEROS+TAURAY       !TAU_ABS_METH
     415
     416      TAUGV_1pt(K)=TAUGV_1pt(K)+TAUGAS*ATERM(NT,K) !INTEG.
     417
     418*  WBARV ET WBARVP
     419*-----------------
     420
     421      IF(TAUGAS+TAEROS+TAURAY+TNUABS .ne.  0.) THEN
     422         WBARV_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 + TNUSCAT)
     423     & /(TAUGAS+TAEROS+TAURAY+TNUABS)
     424      ELSE
     425         WBARV_1pt(J,K,NT)=0.
     426      ENDIF
     427
     428      IF(TAUGAS+TAEROS+TAURAY .ne.  0.) THEN
     429         WBARVP_1pt(J,K,NT)=(TAEROSCAT+TAURAY*0.9999999 )
     430     & /(TAUGAS+TAEROS+TAURAY)
     431      ELSE
     432         WBARVP_1pt(J,K,NT)=0.
     433      ENDIF
    346434
    347435 909  CONTINUE
    … …  
    358446           DO 119 NT=1,NTERM(K)
    359447           TAUV_1pt(1,K,NT)=0.0
     448           TAUVP_1pt(1,K,NT)=0.0
    360449             DO 119 J=1,NLAYER
    361450             TAUV_1pt(J+1,K,NT)=TAUV_1pt(J,K,NT)+DTAUV_1pt(J,K,NT)
     451             TAUVP_1pt(J+1,K,NT)=TAUVP_1pt(J,K,NT)+DTAUVP_1pt(J,K,NT)
    362452 119     CONTINUE
    363453

  • trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F

    r102 r175  
    77     .            rlat,rlon, tsol,tsoil,
    88     .           albe, solsw, sollw,
    9      .           fder,radsol,
     9     .           fder,radsol,resch4,
    1010     .           tabcntr0,
    1111     .           t_ancien,ancien_ok)
    … …  
    2525      REAL dtime
    2626      INTEGER radpas,chimpas
    27       REAL rlat(klon), rlon(klon)
     27      REAL rlat(klon), rlon(klon)   ! in degrees
    2828      REAL tsol(klon)
    2929      REAL tsoil(klon,nsoilmx)
    … …  
    4040      LOGICAL ancien_ok
    4141
     42      REAL resch4(klon)
     43      INTEGER ig0
     44
    4245      REAL xmin, xmax
    4346c
    … …  
    329332         ENDIF
    330333      ENDIF
     334c     Par defaut on cree 2 bandes de methane au pole Nord et au pole Sud
     335c     (entre 75 et 85 degres de latitude) de 2 metres.
     336c     Les poles sont sec !
     337      resch4(1) = 0.    ! pole nord = 1 point
     338      DO ig0=2,klon
     339          if ((rlat(ig0).ge.75..and.rlat(ig0).le.85.).or.
     340     &        (rlat(ig0).ge.-85.and.rlat(ig0).le.-75.)) then
     341            resch4(ig0) = 2.
     342          else
     343            resch4(ig0) = 0.
     344          endif
     345      ENDDO
     346      resch4(klon) = 0.   ! pole sud = 1 point
     347
     348      ierr = NF_INQ_VARID (nid, "RESCH4", nvarid)
     349      IF (ierr.NE.NF_NOERR) THEN
     350         PRINT*, "phyetat0: Le champ <RESCH4> est absent"
     351         PRINT*, "Pas de reservoir de methane mais je continue..."
     352         PRINT*, "Pour info, je met 2 metres de methane sur 2 bandes"
     353         PRINT*, "comprises entre 75 et 85 degres de latitude dans  "
     354         PRINT*, "chaque hemisphere."
     355      ELSE
     356#ifdef NC_DOUBLE
     357         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, resch4)
     358#else
     359         ierr = NF_GET_VAR_REAL(nid, nvarid,resch4)
     360#endif
     361         IF (ierr.NE.NF_NOERR) THEN
     362            PRINT*, "phyetat0: Lecture echouee pour <reservoir>"
     363            CALL abort
     364         ENDIF
     365      ENDIF
    331366c
    332367c Fermer le fichier:

  • trunk/LMDZ.TITAN/libf/phytitan/phyredem.F

    r102 r175  
    77     .           albedo,
    88     .           solsw, sollw,fder,
    9      .           radsol,
     9     .           radsol,resch4,
    1010     .           t_ancien)
    1111
    … …  
    3434      real fder(klon)
    3535      REAL radsol(klon)
     36      real resch4(klon)
    3637      REAL t_ancien(klon,klev)
    3738c
    … …  
    243244      ierr = NF_REDEF (nid)
    244245#ifdef NC_DOUBLE
     246      ierr = NF_DEF_VAR (nid, "RESCH4", NF_DOUBLE, 1, idim2,nvarid)
     247#else
     248      ierr = NF_DEF_VAR (nid, "RESCH4", NF_FLOAT, 1, idim2,nvarid)
     249#endif
     250      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
     251     .                        "Reservoir CH4 a la surface")
     252      ierr = NF_ENDDEF(nid)
     253#ifdef NC_DOUBLE
     254      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,resch4)
     255#else
     256      ierr = NF_PUT_VAR_REAL (nid,nvarid,resch4)
     257#endif
     258c
     259      ierr = NF_REDEF (nid)
     260#ifdef NC_DOUBLE
    245261      ierr = NF_DEF_VAR (nid, "TANCIEN", NF_DOUBLE, 1, idim3,nvarid)
    246262#else

  • trunk/LMDZ.TITAN/libf/phytitan/physiq.F

    r152 r175  
    7979#include "logic.h"
    8080#include "comorbit.h"
     81#include "microtab.h"
     82#include "diagmuphy.h"
     83#include "itemps.h"
    8184c======================================================================
    8285      LOGICAL ok_mensuel ! sortir le fichier mensuel
    … …  
    395398c Recuperation des TAU du TR
    396399      REAL t_tauhvd(klon,klev),t_khvd(klon,klev)
     400      REAL t_tcld(klon,klev),t_kcld(klon,klev)
     401      REAL t_kcvd(klon,klev)
    397402c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
    398403      INTEGER   ngrid
    … …  
    401406      PARAMETER (NSPECV=24,NSPECI=46,NLAYER=llm)
    402407      REAL TAUHID(ngrid,NLAYER,NSPECI)
     408     &               ,TAUCID(ngrid,NLAYER,NSPECI)
    403409     &               ,TAUGID(ngrid,NLAYER,NSPECI)
    404410     &               ,TAUHVD(ngrid,NLAYER,NSPECV)
     411     &               ,TAUCVD(ngrid,NLAYER,NSPECV)
    405412     &               ,TAUGVD(ngrid,NLAYER,NSPECV)
    406413
    407       COMMON /TAUD/   TAUHID,TAUGID,TAUHVD,TAUGVD
     414      COMMON /TAUD/   TAUHID,TAUCID,TAUGID,TAUHVD,TAUCVD,TAUGVD
     415* common relatifs au nuages
     416      real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4)
     417      integer ncount(ngrid,NLAYER)
     418      COMMON/rnuabar/ncount,rmcbar,xfbar
     419
     420       REAL ch4(klon,jjm+1),dch4(jjm+1)
     421       INTEGER ig0
     422       integer ich4
     423       common/ch4ind/ich4
     424
     425c     flux de chaleur latente d'evaporation CH4
     426      REAL fclat(klon)
     427c     reservoir de surface
     428      REAL,save,allocatable :: reservoir(:)
    408429
    409430c Declaration des constantes et des fonctions thermodynamiques
    … …  
    446467         allocate(solsw(klon),sollw(klon),sollwdown(klon))
    447468         allocate(source(klon,nqmax))
     469         allocate(reservoir(klon))
    448470
    449471         CALL suphec ! initialiser constantes et parametres phys.
    … …  
    458480c Initialiser les compteurs:
    459481c
    460          itap     = 0
    461          itaprad  = 0
    462          itapchim = 0
     482         itap        = 0
     483         itaprad     = 0
     484         itapchim    = 0
     485         ncount(:,:) = 0
    463486
    464487c         
    … …  
    469492     .       rlatd,rlond,ftsol,ftsoil,
    470493     .       falbe, solsw, sollw,
    471      .       dlw,radsol,
     494     .       dlw,radsol,reservoir,
    472495c     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
    473496     .       tabcntr0,
    … …  
    575598c Verifications:
    576599c
    577       if ((nmicro.eq.0).and.(microfi.eq.1)) then
    578         print*,"MICROPHYSIQUE ONLINE, MAIS NMICRO=0..."
    579         stop
    580       endif
     600         IF ((nmicro.eq.0).and.(microfi.eq.1)) THEN
     601           abort_message="MICROPHYSIQUE ONLINE, MAIS NMICRO=0..."
     602           call abort_gcm(modname,abort_message,1)
     603         ENDIF
     604         IF (microfi.lt.1.and.clouds.eq.1) THEN
     605          write(lunout,*)"microfi.lt.1.and.clouds.eq.1"
     606          abort_message =
     607     &    "Impossible de faire des nuages sans microphysique..."
     608          call abort_gcm(modname,abort_message,1)
     609         ENDIF
    581610         IF (nlon .NE. klon) THEN
    582611            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
    … …  
    650679      ENDDO
    651680
     681      rmcbar  = 0.
     682      xfbar   = 0.
    652683         
    653684      ENDIF ! debut
    654685c====================================================================
    655686c======================================================================
     687
     688c   Creer un reservoir de surface infini
     689c
     690      reservoir(:) = 2.
    656691
    657692c Mettre a zero des variables de sortie (pour securite)
    … …  
    10611096
    10621097      if (iflag_trac.eq.1) then
     1098c         call begintime(tt0)
    10631099         call phytrac (debut,lafin,
    10641100     .                   nqmax,nmicro,dtime,appel_chim,dtimechim,
    10651101     .                   paprs,pplay,delp,t,rmu0,fract,pdecli,zls,
    1066      .                   tr_seri,qaer,d_tr_mph,d_tr_kim)
    1067 
    1068         if (microfi.eq.1) then
     1102     .                   yu1,yv1,zzlev,zzlay,ftsol,
     1103     .                   tr_seri,qaer,d_tr_mph,d_tr_kim,
     1104     .                   fclat,reservoir)
     1105
     1106c         call endtime(tt0,tt1)
     1107c         ttphytra=ttphytra+tt1
     1108
     1109c ----- ICI on ajuste radsol en tenant compte du flux de chaleur latente
     1110c       d'evaporation du reservoir.
     1111c       NOTE : c'est pas tres elegant mais ca permet d'eviter d'aller
     1112c              toucher a clmain.
     1113        if (clouds.eq.1) then
     1114          radsol(:) = radsol(:)+fclat(:)    !test pas de flx de chaleur latente
     1115        endif
     1116
     1117        if (microfi.ge.1) then
    10691118         tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
    10701119     .                        + d_tr_mph(:,:,1:nmicro)*dtime
    10711120        endif
     1121c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
     1122c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
     1123c       retourne des traceurs nuls et il y a parfois des valeurs negatives qui trainent.
     1124c       Pour ne diffuser le probleme, on force les valeurs negatives a ZERO.
     1125        DO iq=1,nmicro
     1126          DO i=1,klon
     1127            DO l=1,klev
     1128              if (tr_seri(i,l,iq).lt.0.) then
     1129                 tr_seri(i,l,iq) = 0.
     1130              endif
     1131            ENDDO
     1132          ENDDO
     1133        ENDDO
     1134
     1135c condensation:
     1136c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
     1137        if ((clouds.eq.1.or.(chimi)).and.nqmax.gt.nmicro) then
     1138          tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
     1139     .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
     1140        endif
    10721141        if ((chimi).and.(nqmax.gt.nmicro)) then
    10731142c chimie:
    10741143         tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
    1075 c condensation:
    1076          tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
    1077      .                        + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
    10781144        endif
    10791145
    10801146      endif
     1147
     1148c       ch4=0.
     1149c       do l=1,llm
     1150c         ch4(1,l) = tr_seri(1,l,ich4)
     1151c         do j=2,jjm
     1152c           ig0=1+(j-2)*iim
     1153c           do i=1,iim
     1154c             ch4(j,l)= ch4(j,l)  + tr_seri(ig0+i,l,ich4)/iim
     1155c           enddo
     1156c         enddo
     1157c         ch4(jjm+1,l) = tr_seri(klon,l,ich4)
     1158c       enddo
     1159c       do j=1,jjm+1
     1160c         write(501,*) j,ch4(j,1)
     1161c       enddo
     1162c       do l=1,llm
     1163c         write(502,'(I3,49(ES24.17,1X))') l,
     1164c     &   (ch4(j,l),j=1,jjm+1)
     1165c       enddo
     1166c       write(501,*) ""
     1167c       write(502,*) ""
    10811168
    10821169c------------------
    … …  
    11071194c     print*,"sollw avant radlwsw=",sollw(klon/2)
    11081195c     print*,"avant radlwsw"
     1196
     1197c   ----------------
     1198c   Calcul du rayon moyen des gouttes et des fractions volumique pour le TR
     1199c  ----------------
     1200      IF (clouds.eq.1) THEN
     1201        DO i=1,klon
     1202          DO j=1,klev
     1203            rmcbar(i,j)=rmcbar(i,j)/MAX(FLOAT(ncount(i,j)),1.)
     1204            xfbar(i,j,:)=xfbar(i,j,:)/MAX(FLOAT(ncount(i,j)),1.)
     1205          ENDDO
     1206        ENDDO
     1207      ENDIF
     1208     
     1209c      call begintime(tt0)
    11091210      CALL radlwsw
    11101211     e            (dist, rmu0, fract, dtimerad, zzlev,
    … …  
    11151216     s             sollwdown,
    11161217     s             lwnet, swnet)
     1218c      call endtime(tt0,tt1)
     1219c      ttrad=ttrad+tt1
     1220
    11171221c     print*,"apres radlwsw"
     1222c     mise a zero du rayon moyen des gouttes et des fractions volumique
     1223      IF (clouds.eq.1) THEN
     1224        rmcbar(:,:)  = 0.
     1225        xfbar(:,:,:) = 0.
     1226        ncount(:,:)  = 0
     1227      ENDIF
    11181228
    11191229c     print*,"radsol apres radlwsw=",radsol(klon/2)
    … …  
    11221232      itaprad = 0
    11231233      DO k = 1, klev
    1124       DO i = 1, klon
     1234       DO i = 1, klon
    11251235         dtrad(i,k) = heat(i,k)-cool(i,k)     !K/s
    1126       ENDDO
     1236       ENDDO
    11271237      ENDDO
    11281238c       print*,"heat (K/s) =",heat(klon/2,:)
    … …  
    14261536      ENDDO
    14271537c
    1428 c incrementation de la tendance (repassee en extensif) sur qaer pour sorties
    1429       if (microfi.eq.1) then
    1430       do iq=1,nmicro
    1431          DO l=1,llm
    1432           DO i = 1, klon
    1433             qaer(i,l,iq) = qaer(i,l,iq) +
    1434      .                (d_tr_mph(i,l,iq)*delp(i,l)/RG)*dtime
    1435           ENDDO
    1436          ENDDO
    1437       enddo
    1438       endif   ! microfi
    14391538c=============================================================
    14401539c   Ecriture des sorties
    … …  
    15021601     .      falbe,
    15031602     .      solsw, sollw,dlw,
    1504      .      radsol,
     1603     .      radsol,reservoir,
    15051604c     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
    15061605     .      t_ancien)

  • trunk/LMDZ.TITAN/libf/phytitan/phytrac.F

    r104 r175  
    22     .                   nqmax,nmicro,ptimestep,appkim,dtkim,
    33     .                   pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
    4      .                   lonsol,tr_seri,qaer,d_tr_mph,d_tr_kim)
     4     .                   lonsol,
     5     .                   pu,pv,pzlev,pzlay,ftsol,
     6     .                   tr_seri,qaer,d_tr_mph,d_tr_kim,
     7     .                   fclat,reservoir)
    58
    69c======================================================================
    … …  
    2427c pdecli-------input-R-declinaison en radian
    2528c lonsol-------input-R-longitude solaire en radian
     29c pu-----------input-R-vitesse dans la direction X (de O a E) en m/s (1ere couche)
     30c pv-----------input-R-vitesse Y (de S a N) en m/s                   (1ere couche)
     31c pzlev--------input-R-altitude pour chaque inter-couche (en m)
     32c pzlay--------input-R-altitude pour chaque couche (en m)
     33c ftsol--------input-R-temperature au sol (en K)
    2634c tr_seri------input-R-mass mixing ratio traceurs (kg/kg)
    2735c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
    2836c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
     37c fclat--------output-R-flux de chaleur latente d'evaporation du reservoir CH4 (J/m2/s)
     38c reservoir----outpur-R-un reservoir de surface !!! (m)
    2939c======================================================================
    3040      USE infotrac
    … …  
    3444#include "clesphys.h"
    3545#include "YOMCST.h"
     46#include "microtab.h"
     47#include "varmuphy.h"
     48#include "diagmuphy.h"
     49#include "itemps.h"
    3650
    3751c======================================================================
    … …  
    4458      REAL ptemp(klon,klev)
    4559      REAL pmu0(klon), pfract(klon), pdecli, lonsol
     60      REAL pu(klon),pv(klon)
     61      REAL pzlev(klon,klev+1),pzlay(klon,klev)
     62      REAL ftsol(klon)
    4663      REAL tr_seri(klon,klev,nqmax)
    4764      REAL qaer(klon,klev,nqmax)
    4865      REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
     66      REAL fclat(klon)
     67      REAL reservoir(klon)
    4968
    5069c======================================================================
    5170c Local variables
     71      REAL qaer0(klon,klev,nqmax)
     72
     73c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
     74      INTEGER   ngrid,NLAYER
     75      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
     76      PARAMETER (NLAYER=llm)           ! = klev
     77* common relatifs au nuages
     78      real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4)
     79      integer ncount(ngrid,NLAYER)
     80      COMMON/rnuabar/ncount,rmcbar,xfbar
     81
     82      REAL rcloud(klon,klev,nrad),xfrac(klon,klev,4)
     83
     84      REAL vcl,nuc,xgsn,xmsn,xesn,xasn
     85
     86
     87      ReAL gaz1(klon,klev),gaz2(klon,klev),gaz3(klon,klev)
     88
     89      REAL socccld
    5290
    5391c grandeurs en moyennes zonales
    54       REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev)
     92      REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev),ztsol(jjm+1)
     93      REAL zzlev(jjm+1,klev+1),zzlay(jjm+1,klev)
    5594      REAL ztemp(jjm+1,klev),zmu0(jjm+1),zfract(jjm+1)
    5695      real temp_eq(klev),press_eq(klev)
    5796      REAL zqaer(jjm+1,klev,nqmax)    ! et non nmicro... Permet nmicro=0.
    5897      REAL zqaer0(jjm+1,klev,nqmax)
    59       REAL pdqmfi(jjm+1,klev,nqmax)
     98      REAL zdqmufi(jjm+1,klev,nqmax)
    6099      REAL ychim(jjm+1,klev,nqmax-nmicro)
     100      REAL zgaz1(jjm+1,klev),zgaz2(jjm+1,klev),zgaz3(jjm+1,klev)
     101      REAL zgaz10(jjm+1,klev),zgaz20(jjm+1,klev),zgaz30(jjm+1,klev)
    61102c La saturation n est calculee qu une seule fois: sauvegarde qysat
    62103c La chimie n est pas calculee tous les pas, il faut donc
    … …  
    65106     
    66107      character*10 nomqy(nqmax-nmicro+1)
    67       integer      i,j,l,iq,ig0
     108      integer      i,j,k,l,iq,ig0
    68109     
     110      REAL zprec(jjm+1,5),zsolesp(jjm+1,klev+1,3),
     111     &     zflxesp_i(jjm+1,klev,3)
     112      REAL ztau_drop(jjm+1,klev),ztau_aer(jjm+1,klev,nrad)
     113c
     114c    indice des esp chimiques utilisees dans la microfi 
     115      integer icldch4,icldc2h6,icldc2h2
     116      save icldch4,icldc2h6,icldc2h2
     117     
     118      real fte,ftm,Lvch4
     119
     120      REAL tmp,ex,kmin,kmax,dqsq
     121      REAL dqch4
     122c      REAL ch4(jjm+1),ch4b(jjm+1),dch4(jjm+1),ch4c(jjm+1,llm)
     123c       integer ich4
     124c       common/ch4ind/ich4
     125
    69126c======================================================================
    70127c======================================================================
    … …  
    72129      if (firstcall) then
    73130       allocate(qysat(klev,nqmax-nmicro),pdyfi(jjm+1,klev,nqmax-nmicro))
     131
     132c  -------- Quelques verifications au demarrage sur les tailles des tableaux.
     133         IF (microfi.ge.1) then
     134c        Faire de la microphysique sans traceurs... bon courage !
     135           if (nmicro.le.0) then
     136             print*,"aLeRtE cRiTiQuE !!!"
     137             print*,"Vous faites de la microphysique sans traceurs"
     138             print*,"microphysique..."
     139             print*,"Je m'arrete et vous laisse reflechir !"
     140             stop
     141           endif
     142c        Nombre de traceurs incompatibles avec la microphysique.
     143           if (nmicro.ne.ntype*nrad) then
     144             print*,"aLeRtE cRiTiQuE !!!"
     145             print*,"Nb trac imcompatible avec la microphysique."
     146             print*,nmicro,ntype*nrad
     147             stop
     148           endif
     149         ENDIF
     150
    74151      endif
     152
     153c RAZ des sorties : les moyennes se font directement dans IOIPSL :
     154c
     155          flxesp_i(:,:,:) = 0.
     156          tau_drop(:,:)   = 0.
     157          tau_aer(:,:,:)  = 0.
     158          solesp(:,:,:)   = 0.
     159          prec(:,:)       = 0.
    75160
    76161c-----------------------------------------------------------------------
    … …  
    80165c     print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
    81166
    82       zplev = 0.0
    83       zplay = 0.0
    84       ztemp = 0.0
    85       zqaer = 0.0
    86       ychim = 0.0
    87       zmu0  = 0.0
    88       zfract= 0.0
    89      
     167c   -------------------
     168c   Gestion de la temperature et de la pression :
     169c   soit la chimie est active, soit la microphysique se fait en 2D.
     170      IF (chimi.or.microfi.eq.1) THEN
     171 
     172        zplev = 0.0
     173        zplay = 0.0
     174        zzlev = 0.0
     175        zzlay = 0.0
     176        ztemp = 0.0
     177        zqaer = 0.0
     178        ychim = 0.0
     179        zmu0  = 0.0
     180        zfract= 0.0
     181        zgaz1 = 0.0
     182        zgaz2 = 0.0
     183        zgaz3 = 0.0
     184        zprec = 0.0
     185        zflxesp_i = 0.0
     186        ztau_drop = 0.0
     187        ztau_aer  = 0.0
     188        zsolesp   = 0.0
     189   
    90190       do l=1,llm+1
    91191         zplev(1,l) = pplev(1,l)
     192         zzlev(1,l) = pzlev(1,l)
    92193         do j=2,jjm
    93194            ig0=1+(j-2)*iim
    94195            do i=1,iim
    95196               zplev(j,l) = zplev(j,l) + pplev(ig0+i,l)/iim
     197               zzlev(j,l) = zzlev(j,l) + pzlev(ig0+i,l)/iim
    96198            enddo
    97199         enddo
    98200         zplev(jjm+1,l) = pplev(klon,l)
     201         zzlev(jjm+1,l) = pzlev(klon,l)
    99202       enddo
    100203
    … …  
    102205         ztemp(1,l) = ptemp(1,l)
    103206         zplay(1,l) = pplay(1,l)
     207         zzlay(1,l) = pzlay(1,l)
    104208         do j=2,jjm
    105209            ig0=1+(j-2)*iim
    … …  
    107211               ztemp(j,l) = ztemp(j,l) + ptemp(ig0+i,l)/iim
    108212               zplay(j,l) = zplay(j,l) + pplay(ig0+i,l)/iim
     213               zzlay(j,l) = zzlay(j,l) + pzlay(ig0+i,l)/iim
    109214            enddo
    110215         enddo
    111216         ztemp(jjm+1,l) = ptemp(klon,l)
    112217         zplay(jjm+1,l) = pplay(klon,l)
     218         zzlay(jjm+1,l) = pzlay(klon,l)
    113219         temp_eq  = ztemp((jjm+1)/2,:)
    114220         press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
    115221       enddo
     222
     223      ENDIF   ! chimi or microfi=1
     224
     225c   -----------------------------
     226c   Gestion des variables de la microphysique :
     227c
     228c   -------------------
     229      if (microfi.ge.1) then
     230
     231c   Traceurs microphysiques: passage en extensif: n/kg --> n/m^2 (2D ou 3D passage obligatoire)
     232      DO iq=1,nmicro
     233c        print*,tname(iq)
     234        DO l=1,llm
     235          DO i = 1, klon
     236            qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
     237          ENDDO
     238        ENDDO
     239      ENDDO
     240c     copie du tableau de traceur :
     241      qaer0(:,:,:)=qaer(:,:,:)
     242c
     243c   -------------------
     244c    Extraction des gaz pour les nuages
     245c
     246c     recuperation des indices des gaz qui nous interesse       
     247      if (firstcall) then
     248        if (clouds.eq.1) then
     249          icldch4=-1
     250          icldc2h6=-1
     251          icldc2h2=-1
     252          do i=1,nqmax
     253            if (tname(i).eq."CH4") then
     254              icldch4=i
     255c              ich4=i
     256            elseif (tname(i).eq."C2H6") then
     257              icldc2h6=i
     258            elseif (tname(i).eq."C2H2") then
     259              icldc2h2=i
     260            endif
     261          enddo
     262          if (icldch4 .eq.-1 .or.
     263     &        icldc2h6.eq.-1 .or.
     264     &        icldc2h2.eq.-1 ) then
     265            print*, "Sacrebleu !!!"
     266            print*, "Vous voulez faire des nuages sans gaz."
     267            print*, "Mais vous etes inconscient. Je vais m'arreter la"
     268            print*, "pour vous laisser reflechir au probleme"
     269            STOP
     270          endif
     271        endif    ! clouds=1
     272      endif      ! firstcall
     273
     274c     Saturation et fraction molaire CLOUD
     275c     Calcul des saturations pour les esp chimique de la muphy des nuages.
     276c     On le fait ici pour les sortir dans physiq.F sans avoir a surcharger la routine.
     277c     Elles passent ensuite dans un common pour passer dans les I/O.
     278c
     279c-------------------------------------------
     280      IF (clouds.eq.1) THEN
     281        DO l=1,llm
     282          DO i = 1, klon
     283            call ch4sat(ptemp(i,l),pplay(i,l),tmp) !tmp en kg/kg !
     284            satch4(i,l) = tr_seri(i,l,icldch4)/(tmp*28./16.)
     285
     286            call c2h6sat(ptemp(i,l),pplay(i,l),tmp)
     287            satc2h6(i,l) =tr_seri(i,l,icldc2h6)/(tmp*28./30.)
     288
     289            call c2h2sat(ptemp(i,l),pplay(i,l),tmp)
     290            satc2h2(i,l) =tr_seri(i,l,icldc2h2)/(tmp*28./26.)
     291   
     292          ENDDO
     293        ENDDO
     294
     295c   Copie des gaz (en 3D)  <== UNIQUEMENT SI ON FAIT DES NUAGES
     296        gaz1(:,:) = tr_seri(:,:,icldch4)
     297        gaz2(:,:) = tr_seri(:,:,icldc2h6)
     298        gaz3(:,:) = tr_seri(:,:,icldc2h2)
     299
     300      ENDIF      ! clouds=1
     301       
     302      endif      ! microfi.ge.1
     303
     304c     -------------------
     305c     Si microfi = 1 on est en 2D :
     306c     conversion des inputs de muphys
     307      IF (microfi.eq.1) THEN
    116308
    117309         zmu0(1)   = pmu0(1)
    … …  
    126318         zmu0(jjm+1)   = pmu0(klon)
    127319         zfract(jjm+1) = pfract(klon)
    128 
    129 c TRACEURS MICROPHYSIQUES
    130 
    131       if (microfi.eq.1) then
    132 
    133       do iq=1,nmicro
    134 c        print*,tname(iq)
    135 
    136 c       Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
    137          DO l=1,llm
    138           DO i = 1, klon
    139             qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
    140           ENDDO
    141          ENDDO
    142 
     320c
     321c     traceurs 3D --> 2D
     322c
     323      do iq=1,nqmax
    143324       do l=1,llm
    144325         zqaer(1,l,iq) = qaer(1,l,iq)
    … …  
    152333       enddo
    153334      enddo
    154 
    155       endif   ! microfi
     335c       copie du tableau de traceur
     336        zqaer0(:,:,:) = zqaer(:,:,:)
     337c
     338c      gaz 3D --> 2D    <=== UNIQUEMENT SI ON FAIT DES NUAGES.
     339c
     340        if (clouds.eq.1) then
     341          do l=1,llm
     342            zgaz1(1,l) = gaz1(1,l)
     343            zgaz2(1,l) = gaz2(1,l)
     344            zgaz3(1,l) = gaz3(1,l)
     345            do j=2,jjm
     346              ig0=1+(j-2)*iim
     347              do i=1,iim
     348                zgaz1(j,l) = zgaz1(j,l) + gaz1(ig0+i,l)/iim
     349                zgaz2(j,l) = zgaz2(j,l) + gaz2(ig0+i,l)/iim
     350                zgaz3(j,l) = zgaz3(j,l) + gaz3(ig0+i,l)/iim
     351              enddo
     352            enddo
     353            zgaz1(jjm+1,l) = gaz1(klon,l)
     354            zgaz2(jjm+1,l) = gaz2(klon,l)
     355            zgaz3(jjm+1,l) = gaz3(klon,l)
     356          enddo
     357 
     358          zgaz10=zgaz1
     359          zgaz20=zgaz2
     360          zgaz30=zgaz3
     361        endif ! clouds=1
     362
     363      endif   ! microfi=1
    156364
    157365c AUTRES TRACEURS
    … …  
    170378         ychim(jjm+1,l,iq-nmicro) = tr_seri(klon,l,iq)
    171379       enddo
    172          nomqy(iq-nmicro) = tname(iq)
    173 
     380       nomqy(iq-nmicro) = tname(iq)
    174381c       print*,iq-nmicro,nomqy(iq-nmicro)
    175 
    176382      enddo
    177383      nomqy(nqmax-nmicro+1) = "HV"
    … …  
    190396
    191397c-----------------------------------------------------------------------
    192 c     Appel de la microphysique
     398c     Appel de la microphysique   en 2D/3D !!!!!!
    193399c    --------------------------
    194400
    195        pdqmfi = 0.
    196        zqaer0 = zqaer
    197        
    198401       IF(firstcall) THEN
    199402        print*,'MICROPHYSIQUE ',MICROFI
    200403       ENDIF
    201404
    202        IF(MICROFI.eq.1) THEN
    203 
    204          IF(firstcall) THEN
    205           print*,'OH! On passe dans la microphysique'
     405c       call begintime(tt0)
     406       IF (MICROFI.eq.0) THEN
     407c        PAS DE MICROPHYSIQUE :
     408c        On appelle juste rdf pour creer la grille de rayons.
     409         IF (firstcall) THEN
     410          print*,'MICROPHYSIQUE OFF-LINE',MICROFI
     411           call rdf()
    206412         ENDIF
    207 
    208          CALL pg2(zplev,ztemp,zqaer,pdqmfi    ! tendances 2D, /s
    209      &   ,nmicro,ptimestep,zmu0,zfract,lastcall)
    210 
     413c        NOTES :
     414c        L'appel de rdf ne sert a rien ici mis a part pour le TR. Si cet
     415c        appel a deja lieu dans le TR inutile de le refaire ici.
     416c        Je ne sais pas exactement comment marche les modules en F90
     417c        Mais je recopie les valeurs du common/part/ de rdf pour
     418c        les mettre dans un common interne a la microphysique (voir varmuphy.h)
     419c        DONC J'AI BESOIN D'AVOIR ACCES A L'ANCIEN COMMON !!!
     420c
     421       ELSEIF (MICROFI.eq.1) THEN
     422c      MICROPHYSIQUE 2D :
     423c      Les input/output comportent le prefixe z pour 2D :)
     424         zdqmufi = 0.  ! ne sert que pour chimi pour condensation
     425         call muphys(jjm+1,
     426     &        zplev,zplay,zzlev,zzlay,
     427     &        ztemp,zqaer,zgaz1,zgaz2,zgaz3,
     428     &        nmicro,ptimestep,
     429     &        zmu0,zfract,
     430c -------- sorties diagnostiques
     431     &        zflxesp_i,
     432     &        ztau_drop,ztau_aer,
     433     &        zsolesp,zprec)
    211434       ELSE
    212 
    213         IF(firstcall) THEN
    214           print*,'MICROPHYSIQUE OFF-LINE',MICROFI
    215           if (nmicro.gt.0) then
    216             CALL pg2(zplev,ztemp,zqaer,pdqmfi
    217      &    ,nmicro,ptimestep,zmu0,zfract,lastcall)
    218           endif
    219         ENDIF
    220 
     435c      MICROPHYSIQUE 3D :
     436c      Les input sont des champs 3D directement !
     437         call muphys(klon,
     438     &        pplev,pplay,pzlev,pzlay,
     439     &        ptemp,qaer,gaz1,gaz2,gaz3,
     440     &        nmicro,ptimestep,
     441     &        pmu0,pfract,     
     442c ------ sorties diagnostiques
     443     &        flxesp_i,
     444     &        tau_drop,tau_aer,
     445     &        solesp,prec)
     446c
     447c    NOTES :
     448c    Ici toutes nos sorties sont des champs 3D...(meme les diagnostiques)
     449c    On a rien a faire mis a part copier les dq dans les d_tr
     450c
    221451       ENDIF
    222 
    223        zqaer  = zqaer+pdqmfi*ptimestep
    224        pdqmfi = (zqaer-zqaer0)/ptimestep
     452c       call endtime(tt0,tt1)
     453c       ttmuphys=ttmuphys+tt1
     454
     455c-----------------------------------------------------------------------
     456c     Mise a jour des sorties de muphys
     457c    -------------
     458c       En 2D on copie les sorties de muphys de la grille LATxALT
     459c       sur la grille complete.
     460       IF (microfi.eq.1) THEN
     461c        precipitations
     462         DO l=1,5
     463           prec(1,l) = zprec(1,l)
     464           ig0 = 2
     465           DO j=2,jjm
     466             DO i = 1, iim
     467               prec(ig0,l) = zprec(j,l)
     468               ig0 = ig0 + 1
     469             ENDDO
     470           ENDDO
     471           prec(ig0,l) = zprec(jjm+1,l)
     472         ENDDO
     473c        taux sedimentation
     474         DO l=1,llm
     475c          taux sed goutte
     476           IF (clouds.eq.1) THEN
     477             tau_drop(1,l) = ztau_drop(1,l)
     478             ig0 = 2
     479             DO j=2,jjm
     480               DO i = 1, iim
     481                 tau_drop(ig0,l) = ztau_drop(j,l)
     482                 ig0 = ig0 + 1
     483               ENDDO
     484             ENDDO
     485             tau_drop(ig0,l) = ztau_drop(jjm+1,l)
     486           ENDIF
     487c          taux sed aer
     488           DO iq=1,nrad
     489             tau_aer(1,l,iq)  = ztau_aer(1,l,iq)
     490             ig0 = 2
     491             DO j=2,jjm
     492               DO i = 1, iim
     493                 tau_aer(ig0,l,iq)  = ztau_aer(j,l,iq)
     494                 ig0 = ig0 + 1
     495               ENDDO
     496             ENDDO
     497             tau_aer(ig0,l,iq) = ztau_aer(jjm+1,l,iq)
     498           ENDDO
     499         ENDDO
     500c        flux glace / production glace
     501         IF (clouds.eq.1) THEN
     502           DO iq=1,3
     503             DO l=1,llm+1
     504               if (l.le.llm) flxesp_i(1,l,iq) = zflxesp_i(1,l,iq)
     505               solesp(1,l,iq) = zsolesp(1,l,iq)
     506               ig0 = 2
     507               DO j=2,jjm
     508                 DO i = 1, iim
     509                   if (l.le.llm) flxesp_i(ig0,l,iq)=zflxesp_i(1,l,iq)
     510                   solesp(ig0,l,iq) = zsolesp(1,l,iq)   
     511                   ig0 = ig0 + 1
     512                 ENDDO
     513               ENDDO
     514               if (l.le.llm) flxesp_i(ig0,l,iq)=zflxesp_i(jjm+1,l,iq)
     515               solesp(ig0,l,iq) = zsolesp(jjm+1,l,iq)
     516             ENDDO
     517           ENDDO
     518         ENDIF
     519       ENDIF
    225520       
     521c-----------------------------------------------------------------------
     522c     Gestion des sources
     523c    -------------
     524c
     525       IF (clouds.eq.1) THEN
     526         IF (microfi.eq.1) THEN
     527c          On repasse les gaz en 3D si on a fait de la microphysique en 2D
     528           DO l=1,llm
     529             gaz1(1,l) = zgaz1(1,l)
     530             gaz2(1,l) = zgaz2(1,l)
     531             gaz3(1,l) = zgaz3(1,l)
     532             ig0 = 2
     533             DO j=2,jjm
     534               DO i = 1, iim
     535                 gaz1(ig0,l) = zgaz1(j,l)* gaz1(ig0,l) /zgaz10(j,l)
     536                 gaz2(ig0,l) = zgaz2(j,l)* gaz2(ig0,l) /zgaz20(j,l)
     537                 gaz3(ig0,l) = zgaz3(j,l)* gaz3(ig0,l) /zgaz30(j,l)
     538                 ig0 = ig0 + 1
     539               ENDDO
     540             ENDDO
     541             gaz1(ig0,l) = zgaz1(jjm+1,l)
     542             gaz2(ig0,l) = zgaz2(jjm+1,l)
     543             gaz3(ig0,l) = zgaz3(jjm+1,l)
     544           ENDDO
     545         ENDIF
     546c        Mise a jour du reservoir de CH4 (ie : seul le CH4 remplit le reservoir)
     547         DO i=1,klon
     548            reservoir(i) = reservoir(i)+prec(i,1)
     549         ENDDO
     550c        Calcul des sources :
     551c        ch4=0.
     552c        ch4(1) = gaz1(1,1)
     553c         do j=2,jjm
     554c           ig0=1+(j-2)*iim
     555c           do i=1,iim
     556c             ch4(j)= ch4(j) + gaz1(ig0+i,1)/iim
     557c           enddo
     558c         enddo
     559c         ch4(jjm+1) = gaz1(ig0,1)
     560
     561         CALL sources(klon,klev,ptimestep,z0,
     562     &                pu,pv,pplev,pzlay,pzlev,
     563     &                gaz1,gaz2,gaz3,
     564     &                ftsol,solesp,reservoir)
     565 
     566c        ch4b=0.
     567c        ch4b(1) = gaz1(1,1)
     568c         do j=2,jjm
     569c           ig0=1+(j-2)*iim
     570c           do i=1,iim
     571c             ch4b(j)= ch4b(j) + gaz1(ig0+i,1)/iim
     572c           enddo
     573c         enddo
     574c         ch4b(jjm+1) = gaz1(ig0,1)
     575c         do j=1,jjm+1
     576c           write(499,*) j,ch4(j),ch4b(j)
     577c         enddo
     578c         write(499,*) ""
     579       ENDIF
    226580c-----------------------------------------------------------------------
    227581c     Condensation
    228582c    -------------
    229583
    230       if ((chimi).and.(nqmax.gt.nmicro)) then
    231 
    232 c   tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
     584      IF ((chimi).and.(nqmax.gt.nmicro)) then
     585
     586c   tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
    233587c        print*,'Condensation'
    234588
    … …  
    237591              do j=1,jjm+1
    238592                 if (ychim(j,l,iq).gt.qysat(l,iq)) then
    239            pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
    240      .                            / ptimestep                  ! / dt
     593           zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
     594     .                             / ptimestep                  ! / dt
    241595                 endif
    242596              enddo
    … …  
    249603c eventuellement, modif initiale de la compo
    250604c
    251 c   tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
     605c   tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
    252606c
    253607c     if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
    … …  
    260614c             do j=1,jjm+1
    261615c                if (ychim(j,l,iq).le.0.015) then
    262 c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
     616c          zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
    263617c    .                            / ptimestep                  ! / dt
    264618c                endif
    … …  
    278632c             do j=1,jjm+1
    279633c                if (ychim(j,l,iq).gt.1.e-5) then
    280 c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
     634c          zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
    281635c    .                            / ptimestep                  ! / dt
    282636c                endif
    … …  
    288642c             do j=1,jjm+1
    289643c                if (ychim(j,l,iq).gt.3.e-5) then
    290 c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
     644c          zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
    291645c    .                            / ptimestep                  ! / dt
    292646c                endif
    … …  
    298652c     endif
    299653c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     654         
     655c ----- commentaire de fin (mise a jour des profil de fraction molaire)
    300656         
    301657c-----------------------------------------------------------------------
    … …  
    326682
    327683c TRACEURS MICROPHYSIQUES
    328 
    329       if (microfi.eq.1) then
    330 
    331       DO iq=1,nmicro
    332          DO l=1,llm
    333             d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
    334             ig0          = 2
    335             DO j=2,jjm
     684c                                       
     685c ---> pas de microphysique
     686       IF (microfi.eq.0) THEN
     687         DO iq=1,nmicro
     688           d_tr_mph(:,:,iq)=0.
     689         ENDDO
     690       ENDIF
     691c ---> microphysique 2D
     692      IF (microfi.eq.1) THEN
     693c  on repasse le champ de traceurs en 3D (pas les tendances)
     694         DO iq=1,nmicro
     695           DO l=1,llm
     696             qaer(1,l,iq) = zqaer(1,l,iq)
     697             ig0          = 2
     698             DO j=2,jjm
    336699               DO i = 1, iim
    337                   d_tr_mph(ig0,l,iq)  = pdqmfi(j,l,iq)
    338      &                 *qaer(ig0,l,iq)/zqaer0(j,l,iq)
    339                   ig0             = ig0 + 1
     700c    un petit patch :
     701c    Si la moyenne zonale au depart est "nulle" :
     702c    On a quand meme le droit de produire des traceurs dans la cellule.
     703c    On considere donc que la valeur de sortie 3D correspond a la valeur de sortie 2D.
     704c    Cela permet aussi entre autre d'eviter les NaN pour les traceurs des nuages !
     705c    (au dessus de la tropo pas de nuages donc qaer(nrad+1:ntype*nrad) = 0 !!!)
     706                 IF (zqaer0(j,l,iq).lt.1e-100) THEN
     707                   qaer(ig0,l,iq) = zqaer(j,l,iq)
     708                 ELSE
     709                   qaer(ig0,l,iq) = zqaer(j,l,iq) *
     710     &             qaer0(ig0,l,iq)/zqaer0(j,l,iq)
     711                 ENDIF
     712                 ig0 = ig0 + 1
    340713               ENDDO
    341             ENDDO
    342             d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
    343          ENDDO
    344       ENDDO
    345 
    346       do iq=1,nmicro
    347          DO l=1,llm
    348           DO i = 1, klon
     714             ENDDO
     715             qaer(ig0,l,iq) = zqaer(jjm+1,l,iq)
     716           ENDDO
     717         ENDDO
     718c        La tendances correspond a (qaer-qaer0)/ptimestep
     719         DO iq=1,nmicro
     720           DO i=1,klon
     721             DO l=1,llm
     722               d_tr_mph(i,l,iq) = (qaer(i,l,iq)-qaer0(i,l,iq))/
     723     &                            ptimestep
     724             ENDDO
     725           ENDDO
     726         ENDDO
     727c ---> microphysique 3D
     728       ELSEIF(microfi.gt.1) THEN
     729         DO iq=1,nmicro
     730           DO l=1,llm
     731             DO i = 1, klon
     732               d_tr_mph(i,l,iq)=(qaer(i,l,iq)-qaer0(i,l,iq))/ptimestep
     733             ENDDO
     734           ENDDO
     735         ENDDO
     736
     737         do iq=1,nmicro
     738          DO l=1,llm
     739           DO i = 1, klon
    349740c  Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
    350             d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
     741             d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
     742           ENDDO
    351743          ENDDO
    352          ENDDO
    353       enddo
    354 
    355       endif   ! microfi
     744         enddo
     745
     746       ENDIF   ! microfi
    356747
    357748c AUTRES TRACEURS
    … …  
    360751c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
    361752c          a  d_tr_kim (tendance chimique 3D en /s, passee a physiq)
    362 c  et de pdqmfi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
     753c  et de zdqmufi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
    363754
    364755      DO iq=nmicro+1,nqmax
    365756         DO l=1,llm
    366757            d_tr_kim(1,l,iq) = pdyfi(1,l,iq-nmicro)
    367             d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
     758            d_tr_mph(1,l,iq) = zdqmufi(1,l,iq)
    368759            ig0          = 2
    369760            DO j=2,jjm
    … …  
    371762                  d_tr_kim(ig0,l,iq)  = pdyfi(j,l,iq-nmicro)
    372763     &                 *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
    373                   d_tr_mph(ig0,l,iq)  = pdqmfi(j,l,iq)
     764                  d_tr_mph(ig0,l,iq)  = zdqmufi(j,l,iq)
    374765     &                 *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
    375766                  ig0             = ig0 + 1
    … …  
    377768            ENDDO
    378769            d_tr_kim(ig0,l,iq) = pdyfi(jjm+1,l,iq-nmicro)
    379             d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
     770            d_tr_mph(ig0,l,iq) = zdqmufi(jjm+1,l,iq)
    380771         ENDDO
    381772      ENDDO
    382773
    383774      endif   ! chimi
     775
     776c--------------------------------------------------
     777c  CONDENSATION VIA MICROFI
     778c----------------------
     779c La microphysique avec nuages doit se faire obligatoirement en 3D.
     780c Rien n'empeche de faire la chimie en 2D. Cependant pour prendre en compte la
     781c condensation due a la microfi (en 3D) on recalcule la tendance finale pour
     782c les especes concernees (CH4, C2H6 pour le moment).
     783       IF (microfi.ge.1.and.clouds.eq.1) THEN
     784         DO i=1,klon
     785           DO l=1,klev
     786c     condensation CH4
     787             d_tr_mph(i,l,icldch4)=(gaz1(i,l)-tr_seri(i,l,icldch4))
     788     &                            /ptimestep
     789c     condensation C2H6
     790             d_tr_mph(i,l,icldc2h6)=(gaz2(i,l)-tr_seri(i,l,icldc2h6))
     791     &                             /ptimestep
     792c     condensation C2H2
     793             d_tr_mph(i,l,icldc2h2)=(gaz3(i,l)-tr_seri(i,l,icldc2h2))
     794     &                             /ptimestep
     795           ENDDO
     796         ENDDO
     797       ENDIF
     798c        ch4c=0.
     799c       do l=1,llm
     800c       ch4c(1,l) = tr_seri(1,l,icldch4)
     801c       do j=2,jjm
     802c         ig0=1+(j-2)*iim
     803c         do i=1,iim
     804c            ch4c(j,l)= ch4c(j,l)+tr_seri(ig0+i,l,icldch4)/iim
     805c         enddo
     806c       enddo
     807c        ch4c(jjm+1,l) = tr_seri(klon,l,icldch4)
     808c      enddo
     809c       do l=1,llm
     810c         write(500,*) pplay(25,l),ch4c(25,l)
     811c       enddo
     812c       write(500,*) ""
     813
     814
     815c--------------------------------------------------
     816c  MISE A JOUR CH4 : (pour refixer la fraction
     817c                     molaire)
     818c--------------------------------------------------
     819c       IF (firstcall) THEN
     820c         do i=1,klon
     821c           do j=1,llm
     822c             call ch4sat(ptemp(i,j),pplay(i,j),tmp) !tmp en kg/kg !
     823c             tmp=0.95*0.85*tmp*28./16.
     824c             if (pplay(i,j).lt.20000.) then
     825c               dqch4 = 1.4e-2         
     826c             else
     827c               dqch4 = tmp
     828c             endif
     829c             d_tr_mph(i,j,icldch4)=(-tr_seri(i,j,icldch4)+dqch4)/
     830c     &       ptimestep
     831c           enddo
     832c         enddo
     833c         
     834c       ENDIF
     835
     836c--------------------------------------------------
     837c CALCUL DU FLUX DE CHALEUR LATENTE D'EVAPORATION
     838c DU METHANE
     839c--------------------------------------------------
     840       IF (clouds.eq.1) THEN
     841         DO i=1,klon
     842           fte= (1.-ftsol(i)/305.5)
     843           ftm= (1.-ftsol(i)/190.5)
     844           if(ftm.le.1.e-3) ftm=1.e-3
     845           if(fte.le.1.e-3) fte=1.e-3
     846           Lvch4 =8.314*190.4*
     847     &     (7.08*ftm**0.354+10.95*1.1e-2*ftm**0.456)
     848     &     /mch4
     849           ! solcxhy en m3/m2 {ok}
     850           ! 425 en kg/m3
     851           ! Lv en J/kg       {ok}
     852           ! ptimestep en s   {ok}
     853           fclat(i)=(solesp(i,klev+1,1)*Lvch4*rhoi_ch4)   ! en J/m2/s
     854         ENDDO
     855       ENDIF
     856
     857c--------------------------------------------------
     858c  GESTION DES RAYONS DE GOUTTES POUR TR
     859c--------------------------------------------------
     860       IF (clouds.eq.1) THEN
     861
     862c Calcul du rayon des gouttes par bin ...
     863c----------------------------------------
     864         DO i=1,klon
     865           DO j=1,klev
     866             DO iq=1,nrad
     867*      Rayon minimum selon la quantité de noyaux
     868               IF (qaer(i,j,iq+nrad) .le.   1.e-5) THEN
     869                  rcloud(i,j,iq) = 1.e-10
     870               ELSE
     871                 rcloud(i,j,iq)=
     872     &           ((qaer(i,j,iq+2*nrad)/qaer(i,j,iq+nrad)+
     873     &           qaer(i,j,iq+3*nrad)/qaer(i,j,iq+nrad) +
     874     &           v_e(iq))*0.75/3.1415926353)**(1./3.)
     875               ENDIF
     876             ENDDO
     877           ENDDO
     878         ENDDO
     879
     880c .... et de leur rayon moyen total (tt bins confondu)
     881c------------------------------------------------------
     882         DO i=1,klon
     883           socccld=0.
     884           DO j=klev,1,-1    !de haut en bas pour le calcul des opacites
     885             vcl=0.
     886             nuc=0.
     887             xgsn=0.
     888             xmsn=0.
     889             xesn=0.
     890             xasn=0.
     891             DO iq=1,nrad
     892               vcl=vcl+qaer(i,j,iq+2*nrad)+
     893     &         qaer(i,j,iq+3*nrad)+
     894     &         qaer(i,j,iq+4*nrad)+
     895     &         v_e(iq)*qaer(i,j,iq+nrad)            ! volume des gouttes
     896               nuc=nuc+qaer(i,j,iq+nrad)            ! nombre de noyaux
     897               xgsn=xgsn+qaer(i,j,iq+nrad)*v_e(iq)  ! volume de noyaux
     898               xmsn=xmsn+qaer(i,j,iq+2*nrad)        ! volume de methane
     899               xesn=xesn+qaer(i,j,iq+3*nrad)        ! volume d' ethane
     900               xasn=xasn+qaer(i,j,iq+4*nrad)        ! volume d' acethylene
     901             ENDDO
     902             IF (nuc .le.  1.e-5) THEN
     903               rmcloud(i,j)=1.0e-10
     904               xfrac(i,j,:)=0.
     905             ELSE
     906               IF(xgsn/vcl.lt.0.  .or. xgsn/vcl.gt.1.001)
     907     &         print*, 'PB AVEC XFRAC:', i,j,xgsn,vcl
     908               rmcloud(i,j)=          ! rayon moyen des gouttes
     909     &         (vcl/nuc*0.75/3.1415926353)**(1./3.)
     910               xfrac(i,j,1)=xgsn/vcl         ! fraction volumique noyau/goutte
     911               xfrac(i,j,2)=xmsn/vcl         ! fraction volumique CH4/goutte
     912               xfrac(i,j,3)=xesn/vcl         ! fraction volumique C2H6/goutte
     913               xfrac(i,j,4)=xasn/vcl         ! fraction volumique C2H2/goutte
     914c              calcul du rayon moyen (moyenne temporelle)
     915               rmcbar(i,j)=rmcbar(i,j)+rmcloud(i,j)
     916               xfbar(i,j,:)=xfbar(i,j,:)+xfrac(i,j,:)
     917               ncount(i,j) = ncount(i,j)+1
     918             ENDIF
     919             socccld=socccld+3.1415926*(rmcloud(i,j)**2.)*nuc
     920             occcld(i,j)=socccld
     921           ENDDO
     922         ENDDO
     923c
     924c      OCCCLD
     925c      Calcul le nombre d'occurence d'un nuage
     926c      d'opacité comprise en kmin et kmax
     927c          k        kmin            kmax
     928c          1   0.0000000      0.10000000   
     929c          2  0.10000000      0.17782794   
     930c          3  0.17782794      0.31622776   
     931c          4  0.31622776      0.56234139   
     932c          5  0.56234139       1.0000000   
     933c          6   1.0000000       1.7782795   
     934c          7   1.7782795       3.1622777   
     935c          8   3.1622777       5.6234136   
     936c          9   5.6234136       10.000000   
     937c         10   10.000000       17.782795   
     938c         11   17.782795       31.622778   
     939c         12   31.622778       100000.00
     940c
     941         DO i=1,klon
     942           DO j=1,klev
     943             DO k=1,12
     944               ex=10.**(0.25)
     945               kmin=0.
     946               kmax=1.e5
     947               if(k.ne.1)  kmin=0.1*ex**(k-2)
     948               if(k.ne.12) kmax=0.1*ex**(k-1)
     949               if(occcld(i,j).ge.kmin .and. occcld(i,j).lt.kmax)
     950     &         occcld_m(i,j,k)=1.
     951             ENDDO
     952           ENDDO
     953         ENDDO
     954       ENDIF  ! fin condition clouds => pas besoin de calculer des rayons
    384955
    385956      RETURN

  • trunk/LMDZ.TITAN/libf/phytitan/radlwsw.F

    r121 r175  
    5959      real solsw(klon), sollw(klon)
    6060      real sollwdown(klon)
    61       REAL swnet(klon,kflev+1),lwnet(klon,kflev+1)
     61      REAL swnet(klon,klev+1),lwnet(klon,klev+1)
    6262c
    6363c LOCAL VARIABLES
    … …  
    6565      real zp(klon,klev+1),zt(klon,klev+1),zz(klon,klev+1)
    6666      real zq(klon,klev,nq)
    67       real zheat(klon,klev), zcool(klon,klev)
    68       REAL zswnet(klon,kflev+1),zlwnet(klon,kflev+1)
    69      
     67      real zheatc(klon,klev), zcoolc(klon,klev)
     68      real zheatp(klon,klev), zcoolp(klon,klev)
     69      REAL zswnetc(klon,klev+1),zlwnetp(klon,klev+1)
     70      REAL zswnetp(klon,klev+1),zlwnetc(klon,klev+1)
     71      REAL zsollwdownc(klon),zsollwdownp(klon)
     72      INTEGER icld
     73
    7074
    7175c =======================================
    … …  
    123127      print*,'On calcule le rayonnement SW'
    124128
    125       CALL heating(dist,rmu0,fract,zheat,zswnet)
     129       IF (clouds.eq.1) THEN
     130         ICLD = 1   ! colonne avec nuages
     131         CALL heating(dist,rmu0,fract,zheatc,zswnetc,icld)
     132       ELSE
     133         zheatc  = 0.
     134         zswnetc = 0.
     135       ENDIF
     136       ICLD = 0   ! colonne sans nuages
     137       CALL heating(dist,rmu0,fract,zheatp,zswnetp,icld)
    126138
    127139c inversion de l'axe vertical
    128              do l=1,klev
    129                 do i=1,klon
    130                    heat(i,l)=zheat(i,klev+1-l)
    131                 enddo
    132              enddo
    133              do l=1,klev+1
    134                 do i=1,klon
    135                    swnet(i,l)=zswnet(i,klev+2-l)
    136                 enddo
    137              enddo
     140       do l=1,klev
     141         do i=1,klon
     142           heat(i,l)=zheatc(i,klev+1-l)*xnuf +
     143     &               zheatp(i,klev+1-l)*(1.-xnuf)
     144         enddo
     145       enddo
     146       do l=1,klev+1
     147         do i=1,klon
     148           swnet(i,l)=zswnetc(i,klev+2-l)*xnuf +
     149     &                zswnetp(i,klev+2-l)*(1.-xnuf)
     150         enddo
     151       enddo
    138152
    139153      solsw = swnet(:,1)
    … …  
    146160      print*,'On calcule le rayonnement LW'
    147161
    148       CALL cooling(klev+1,zp,zt,zz,zcool,zlwnet,sollwdown)
     162       IF (clouds.eq.1) THEN
     163         ICLD = 1
     164         CALL cooling(klon,klev+1,zp,zt,zz,zcoolc,zlwnetc,zsollwdownc,
     165     &   icld)
     166       ELSE
     167         zcoolc      = 0.
     168         zlwnetc     = 0.
     169         zsollwdownc = 0.
     170       ENDIF
     171       ICLD = 0
     172       CALL cooling(klon,klev+1,zp,zt,zz,zcoolp,zlwnetp,zsollwdownp,
     173     & icld)
    149174
    150175c inversion de l'axe vertical
    151              do l=1,klev
    152                 do i=1,klon
    153                    cool(i,l)=zcool(i,klev+1-l)
    154                 enddo
    155              enddo
    156              do l=1,klev+1
    157                 do i=1,klon
    158                    lwnet(i,l)=zlwnet(i,klev+2-l)
    159                 enddo
    160              enddo
     176       do l=1,klev
     177         do i=1,klon
     178           cool(i,l)=zcoolc(i,klev+1-l)*xnuf +
     179     &               zcoolp(i,klev+1-l)*(1.-xnuf)
     180         enddo
     181       enddo
     182       do l=1,klev+1
     183         do i=1,klon
     184           lwnet(i,l)=zlwnetc(i,klev+2-l)*xnuf +
     185     &                zlwnetp(i,klev+2-l)*(1.-xnuf)
     186         enddo
     187       enddo
     188   
     189       do i=1,klon
     190         sollwdown(i)=zsollwdownc(i)*xnuf +
     191     &                zsollwdownp(i)*(1.-xnuf)
     192       enddo
    161193
    162194      sollw  = -lwnet(:,1)

  • trunk/LMDZ.TITAN/libf/phytitan/radtitan.F

    r106 r175  
    8181c   ---------------------------------------------
    8282
    83       REAL DTAUP(ngrid,NLAYER,NSPECI)
    84       REAL UBARI,UBARV,UBAR0
    8583      REAL DZED(NLAYER)
    8684      REAL Z(NLEVEL),PRESS(NLEVEL),DEN(NLEVEL),TEMP(NLEVEL)
    87       REAL  DTDP(NLAYER),CONVEQ
    8885      REAL  CH4(NLEVEL),XN2(NLEVEL),H2(NLEVEL),AR(NLEVEL)
    8986      REAL  XMU(NLEVEL),GAS1(NLAYER),COLDEN(NLAYER)
    9087      REAL  C2H2(NLAYER),C2H6(NLAYER),HCN(NLAYER)
    91       REAL  SOLARF(NSPECV),PEXPON(NSPECV),ATERM(4,NSPECV)
    92       INTEGER NTERM(NSPECV)
    93       REAL  BTERM(4,NSPECV)
    94       REAL  RADIUS(NLAYER), XNUMB(NLAYER),REALI(NSPECI), XIMGI(NSPECI)
    95       REAL  REALV(NSPECV), XIMGV(NSPECV)
    96       REAL  RADCLD(NLAYER), XNCLD(NLAYER),RCLDI(NSPECI),  XICLDI(NSPECI)
    97       REAL  RCLDV(NSPECV),  XICLDV(NSPECV)
    98       REAL  TAUHI(ngrid,NSPECI),TAUCI(ngrid,NSPECI)
    99       REAL  TAUGI(ngrid,NSPECI), TAUGV(ngrid,NSPECV)
    100       REAL  TAURV(ngrid,NSPECV),TAUHV(ngrid,NSPECV)
    101       REAL  TAUCV(ngrid,NSPECV)
    102 c
    103       REAL  DTAUI(ngrid,NLAYER,NSPECI)
    104       REAL  TAUI(ngrid,NLEVEL,NSPECI)
    105       REAL  WBARI(ngrid,NLAYER,NSPECI)
    106       REAL  COSBI(ngrid,NLAYER,NSPECI)
    107       REAL  BWNI(NSPC1I),WNOI(NSPECI)
    108       REAL  WLNI(NSPECI),DWNI(NSPECI)
    109 c
    110       REAL DTAUV(ngrid,NLAYER,NSPECV,4)
    111       REAL TAUV(ngrid,NLEVEL,NSPECV,4)
    112       REAL WBARV(ngrid,NLAYER,NSPECV,4)
    113       REAL COSBV(ngrid,NLAYER,NSPECV,4)
    114       REAL BWNV(NSPC1V), WNOV(NSPECV),DWNV(NSPECV), WLNV(NSPECV)
    115       REAL FNETV(ngrid,NLEVEL),FUPV(ngrid,NLEVEL,NSPECV) 
    116       REAL FDV(ngrid,NLEVEL,NSPECV),FMNETV(ngrid,NLEVEL)
    117       REAL FMUPV(NLEVEL),FMDV(NLEVEL)
    118       REAL FNET(ngrid,NLEVEL),FMNET(ngrid,NLEVEL)
    119       REAL THEAT(ngrid,NLAYER)
    120       REAL CSUBP,RSFI,RSFV,F0PI
    12188      REAL  RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
    122       REAL RGAS,RHOP,PII,SIGMA
    123       REAL TIDAL
    124 
    125       COMMON /IRTAUS/ DTAUP
     89      REAL RGAS,RHOP,PI,SIGMA
    12690
    12791      COMMON /VERTICAL/ DZED
    … …  
    12993      COMMON /ATM/ Z,PRESS
    13094     &            ,DEN,TEMP
    131 
    132       COMMON /LAPSE/ DTDP,CONVEQ
    133       COMMON /UBARED/ UBARI,UBARV,UBAR0
    134 
    13595
    13696
    … …  
    143103      COMMON /STRAT2/ HCN
    144104
    145       COMMON /VISGAS/ SOLARF,NTERM
    146      &               ,PEXPON,ATERM
    147      &               ,BTERM
    148 
    149       COMMON /AERSOL/ RADIUS, XNUMB
    150      &               ,REALI, XIMGI
    151      &               ,REALV, XIMGV
    152 
    153       COMMON /CLOUD/ RADCLD, XNCLD
    154      &             , RCLDI,  XICLDI
    155      &             , RCLDV,  XICLDV
    156 
    157       COMMON /TAUS/   TAUHI,TAUCI,
    158      &                TAUGI, TAUGV,
    159      &                TAURV,TAUHV,
    160      &                TAUCV
    161 
    162 *     INFRARED  CHARACTERISTICS
    163 *------------------------------
    164 
    165 
    166       COMMON /SPECTI/ BWNI,WNOI
    167      &               ,DWNI,WLNI
    168 
    169       COMMON /OPTICI/ DTAUI,
    170      &                TAUI,
    171      &                WBARI,
    172      &                COSBI
    173 
    174 
    175 
    176 *     VISIBLE  CHARACTERISTICS
    177 *------------------------------
    178  
    179 
    180 
    181       COMMON /OPTICV/ DTAUV
    182      &               ,TAUV
    183      &               ,WBARV
    184      &               ,COSBV
    185 
    186       COMMON /SPECTV/ BWNV, WNOV
    187      &               ,DWNV, WLNV
    188 
    189       COMMON /FLUXvV/ FNETV,     
    190      &               FUPV,
    191      &               FDV,
    192      &               FMNETV
    193 
    194       COMMON /FLUX/  FNET,      FMNET
    195      &              ,THEAT
    196 
    197 
    198       COMMON /PLANT/ CSUBP,RSFI,RSFV,F0PI
    199105      COMMON /ADJUST/ RHCH4,FH2,FHAZE,FHVIS,FHIR,TAUFAC,RCLOUD,FARGON
    200       COMMON /CONST/RGAS,RHOP,PII,SIGMA
    201       COMMON /IO/ TIDAL
    202 
    203 * nrad dans microtab.h
    204       REAL volume(nrad),rayon(nrad),vrat,
    205      &      drayon(nrad),dvolume(nrad)
    206 
    207       common/part/volume,rayon,vrat,
    208      &      drayon,dvolume
     106      COMMON /CONST/RGAS,RHOP,PI,SIGMA
     107
    209108c-----------------------------------------------------------------------
    210109c   1. Initialisations:
    … …  
    368267      print*,'On sort de optci'
    369268
    370 C THIS ROUTINE HAS ALREADY SET THE DTAUI(J,K) VALUES BUT MUST BE PASSED
    371         DO 225 IG=1,klon
    372          DO 220 J=1,NLAYER
    373           DO 230 K=1,NSPECI
    374               DTAUP(IG,J,K)=DTAUI(IG,J,K)
    375 230       CONTINUE
    376 220      CONTINUE
    377 225     CONTINUE
    378 
    379269C  NOW, THIS COMPUTATION IS DONE FOR EACH VALUE OF klon
    380270C  INFRARED. AND THEN IN THE VISIBLE.
    … …  
    397287c on ne recalcule pas optci si microfi=0 et compo lellouch
    398288c -----------------------------
    399       IF ((MICROFI.eq.1).or.(.not.ylellouch)) THEN 
     289      IF ((MICROFI.ge.1).or.(.not.ylellouch)) THEN 
    400290      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
    401291       print*,'aerosol/gas/cloud properties'
    402292       CALL OPTCI(ycomp,qaer,nmicro,IPRINT)        ! #1
    403         DO  IG=1,klon
    404          DO  J=1,NLAYER
    405           DO  K=1,NSPECI
    406               DTAUP(IG,J,K)=DTAUI(IG,J,K)
    407           ENDDO
    408          ENDDO
    409         ENDDO
    410293      ENDIF
    411294      ENDIF
    … …  
    413296c ni optcv si microfi=0
    414297
    415       IF (MICROFI.eq.1) THEN 
     298      IF (MICROFI.ge.1) THEN 
    416299      IF(notfirstcall)    THEN  !F au 1er appel   T aux autres appels!!
    417300       print*,'aerosol/gas/cloud properties'

  • trunk/LMDZ.TITAN/libf/phytitan/sfluxv.F

    r104 r175  
    1       SUBROUTINE SFLUXV(IPRINT,IG,dist_sol)
     1      SUBROUTINE SFLUXV(IPRINT,IG,dist_sol,icld)
    22
    33      use dimphy
    … …  
    1010      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
    1111c
    12       INTEGER NLAYER,NLEVEL,NSPECV,NSPC1V
     12      INTEGER NLAYER,NLEVEL,NSPECV,NSPC1V,icld
    1313      PARAMETER (NLAYER=llm,NLEVEL=NLAYER+1)
    1414      PARAMETER (NSPECV=24,NSPC1V=25)
    … …  
    2424     &     ,WBARV(ngrid,NLAYER,NSPECV,4)
    2525     &     ,COSBV(ngrid,NLAYER,NSPECV,4)
     26     &     ,DTAUVP(ngrid,NLAYER,NSPECV,4)
     27     &     ,TAUVP(ngrid,NLEVEL,NSPECV,4)
     28     &     ,WBARVP(ngrid,NLAYER,NSPECV,4)
     29     &     ,COSBVP(ngrid,NLAYER,NSPECV,4)
    2630      REAL BWNV(NSPC1V),WNOV(NSPECV)
    2731     &     ,DWNV(NSPECV),WLNV(NSPECV)
    … …  
    4347     &                ,WBARV
    4448     &                ,COSBV
     49     &                ,DTAUVP
     50     &                ,TAUVP
     51     &                ,WBARVP
     52     &                ,COSBVP
    4553
    4654      COMMON /SPECTV/ BWNV,WNOV
    … …  
    8694C LOOP OVER THE NTERMS BEGINING HERE
    8795      DO 912 NT=1,NTERM(K)
    88       BSURF=0.+ RSFV*UBAR0*F0PI*EXP(-TAUV(ig,NLEVEL,K,NT)/UBAR0)
     96      IF (ICLD.eq.1) THEN
     97        BSURF=0.+ RSFV*UBAR0*F0PI*EXP(-TAUV(ig,NLEVEL,K,NT)/UBAR0)
     98      ELSE
     99        BSURF=0.+ RSFV*UBAR0*F0PI*EXP(-TAUVP(ig,NLEVEL,K,NT)/UBAR0)
     100      ENDIF
    89101C
    90102C* WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
    … …  
    99111C  USE DT0,T0,WB0,CO0 INSTEAD OF DTAUV(ig,1,K,NT)..etc...
    100112
    101        DO  J=1,NLAYER       
    102              DT0(J)=DTAUV(ig,J,K,NT)
    103              T0(J) =TAUV(ig,J,K,NT)
    104              WB0(J)=WBARV(ig,J,K,NT)
    105              CO0(J)=COSBV(ig,J,K,NT)
    106        ENDDO
     113       IF (ICLD.EQ.1) THEN
     114         DO  J=1,NLAYER       
     115           DT0(J)=DTAUV(ig,J,K,NT)
     116           T0(J) =TAUV(ig,J,K,NT)
     117           WB0(J)=WBARV(ig,J,K,NT)
     118           CO0(J)=COSBV(ig,J,K,NT)
     119         ENDDO
     120         T0(NLEVEL)=TAUV(ig,NLEVEL,K,NT)
     121       ELSE
     122         DO  J=1,NLAYER       
     123           DT0(J)=DTAUVP(ig,J,K,NT)
     124           T0(J) =TAUVP(ig,J,K,NT)
     125           WB0(J)=WBARVP(ig,J,K,NT)
     126           CO0(J)=COSBVP(ig,J,K,NT)
     127         ENDDO
     128         T0(NLEVEL)=TAUVP(ig,NLEVEL,K,NT)
    107129
    108              T0(NLEVEL)=TAUV(ig,NLEVEL,K,NT)
     130       ENDIF
    109131
    110132c       PRINT*,'entree gfluxv #: ',ig,K

  • trunk/LMDZ.TITAN/libf/phytitan/write_histday.h

    r110 r175  
    9494     .                                   iim*jjmp1*klev,ndex3d)
    9595c
    96       ENDIF !lev_histday.GE.2
    97 c
    98 c-------------------------------------------------------
    99       IF(lev_histday.GE.3) THEN
    100 c
    10196cccccccccccccccccc  Tracers
    10297c
    10398         if (iflag_trac.eq.1) THEN
    104           if (microfi.eq.1) then
    105            DO iq=1,nmicro
    106        CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qaer(1,1,iq), zx_tmp_3d)
    107        CALL histwrite(nid_day,tname(iq),itau_w,zx_tmp_3d,
    108      .                                   iim*jjmp1*klev,ndex3d)
    109            ENDDO
     99          if (microfi.ge.1) then
     100c           DO iq=1,nmicro
     101c      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qaer(1,1,iq), zx_tmp_3d)
     102c      CALL histwrite(nid_day,tname(iq),itau_w,zx_tmp_3d,
     103c     .                                   iim*jjmp1*klev,ndex3d)
     104c           ENDDO
     105c    -------   NB AER TOT
     106               do i=1,klon
     107                 do j=1,klev
     108                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,1:nrad))
     109                 enddo
     110               enddo
     111       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     112       CALL histwrite(nid_day,"qaer",itau_w,zx_tmp_3d,
     113     .                                   iim*jjmp1*klev,ndex3d)
     114c    -------   NB NOY TOT
     115               do i=1,klon
     116                 do j=1,klev
     117                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,nrad+1:2*nrad))
     118                 enddo
     119               enddo
     120       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     121       CALL histwrite(nid_day,"qnoy",itau_w,zx_tmp_3d,
     122     .                                   iim*jjmp1*klev,ndex3d)
     123c    -------   V GLA1 TOT
     124               do i=1,klon
     125                 do j=1,klev
     126                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,2*nrad+1:3*nrad))
     127                 enddo
     128               enddo
     129       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     130       CALL histwrite(nid_day,"qgl1",itau_w,zx_tmp_3d,
     131     .                                   iim*jjmp1*klev,ndex3d)
     132c    -------   V GLA2 TOT
     133               do i=1,klon
     134                 do j=1,klev
     135                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,3*nrad+1:4*nrad))
     136                 enddo
     137               enddo
     138       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     139       CALL histwrite(nid_day,"qgl2",itau_w,zx_tmp_3d,
     140     .                                   iim*jjmp1*klev,ndex3d)
     141c    -------   V GLA3 TOT
     142               do i=1,klon
     143                 do j=1,klev
     144                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,4*nrad+1:5*nrad))
     145                 enddo
     146               enddo
     147       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     148       CALL histwrite(nid_day,"qgl3",itau_w,zx_tmp_3d,
     149     .                                   iim*jjmp1*klev,ndex3d)
     150c --------------
     151c ----- SATURATION ESP NUAGES
     152             if (clouds.eq.1) then
     153
     154       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, satch4,zx_tmp_3d)
     155       CALL histwrite(nid_day,"ch4sat", itau_w, zx_tmp_3d,
     156     .                                   iim*jjmp1*klev,ndex3d)
     157
     158       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, satc2h6,zx_tmp_3d)
     159       CALL histwrite(nid_day,"c2h6sat", itau_w, zx_tmp_3d,
     160     .                                   iim*jjmp1*klev,ndex3d)
     161
     162       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, satc2h2,zx_tmp_3d)
     163       CALL histwrite(nid_day,"c2h2sat", itau_w, zx_tmp_3d,
     164     .                                   iim*jjmp1*klev,ndex3d)
     165c --------------
     166c ----- RESERVOIR DE SURFACE
     167       CALL gr_fi_ecrit(1, klon,iim,jjmp1,reservoir,zx_tmp_2d)
     168       CALL histwrite(nid_day,"reserv",itau_w,zx_tmp_2d,
     169     .                        iim*jjmp1,ndex2d)
     170c --------------
     171c ----- PRECIPITATIONS
     172c       -----  CH4
     173       CALL gr_fi_ecrit(1, klon,iim,jjmp1,prec(:,1),zx_tmp_2d)
     174       CALL histwrite(nid_day,"prech4",itau_w,zx_tmp_2d,
     175     .                        iim*jjmp1,ndex2d)
     176c       -----  C2H6
     177       CALL gr_fi_ecrit(1, klon,iim,jjmp1,prec(:,2),zx_tmp_2d)
     178       CALL histwrite(nid_day,"prec2h6",itau_w,zx_tmp_2d,
     179     .                        iim*jjmp1,ndex2d)
     180c       -----  C2H2
     181       CALL gr_fi_ecrit(1, klon,iim,jjmp1,prec(:,3),zx_tmp_2d)
     182       CALL histwrite(nid_day,"prec2h2",itau_w,zx_tmp_2d,
     183     .                        iim*jjmp1,ndex2d)
     184c
     185c --------------
     186c ----- FLUX GLACE
     187       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,flxesp_i(1,1,1),zx_tmp_3d)
     188       CALL histwrite(nid_day,"flxgl1", itau_w, zx_tmp_3d,
     189     .                                   iim*jjmp1*klev,ndex3d)
     190       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,flxesp_i(1,1,2),zx_tmp_3d)
     191       CALL histwrite(nid_day,"flxgl2", itau_w, zx_tmp_3d,
     192     .                                   iim*jjmp1*klev,ndex3d)
     193       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,flxesp_i(1,1,3),zx_tmp_3d)
     194       CALL histwrite(nid_day,"flxgl3", itau_w, zx_tmp_3d,
     195     .                                   iim*jjmp1*klev,ndex3d)
     196c
     197c --------------
     198c ----- RAYON MOYEN GOUTTE
     199       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, rmcloud,zx_tmp_3d)
     200       CALL histwrite(nid_day,"rcldbar", itau_w, zx_tmp_3d,
     201     .                                   iim*jjmp1*klev,ndex3d)
     202c
     203             endif
    110204          endif
     205c
     206c --------------
     207c ----- TRACEURS CHIMIQUES
    111208          if (nmicro.lt.nqmax) then
    112209           DO iq=nmicro+1,nqmax
    … …  
    118215         endif
    119216c
     217      ENDIF !lev_histday.GE.2
     218c
     219c-------------------------------------------------------
     220      IF(lev_histday.GE.3) THEN
     221c
    120222cccccccccccccccccc  Radiative transfer
    121223c
    … …  
    143245     .                                   iim*jjmp1*klev,ndex3d)
    144246c
    145 c 3D adding Tau and k  (31/08/10)
    146 c
     247c --------------
     248c ----- OPACITE BRUME
    147249       do k=7,NSPECV,10
    148250         do i=1,klon
    … …  
    158260       enddo      ! fin boucle NSPECV
    159261
     262       do k=8,NSPECI,10
     263         do i=1,klon
     264         do l=1,klev
     265           t_tauhvd(i,l)=TAUHID(i,klev-l+1,k)
     266         enddo
     267         enddo
     268         write(str1,'(i2.2)') k
     269      zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
     270      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     271      CALL histwrite(nid_day,"thi"//str1,itau_w,zx_tmp_3d,
     272     .                                   iim*jjmp1*klev,ndex3d)
     273       enddo      ! fin boucle NSPECI
     274c
     275c --------------
     276c ----- EXTINCTION BRUME
    160277       do k=7,NSPECV,10
    161278         do i=1,klon
    162279         do l=1,klev
    163          if(l.ne.klev)
    164      s    t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
    165      s    -TAUHVD(i,klev-l+1-1,k)
    166 
     280          if(l.ne.klev)
     281     s     t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
     282     s                -TAUHVD(i,klev-l+1-1,k)
    167283          if(l.eq.klev)
    168      s    t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
     284     s     t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
    169285
    170286         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    … …  
    178294       enddo      ! fin boucle NSPECV
    179295
     296       do k=8,NSPECI,10
     297         do i=1,klon
     298         do l=1,klev
     299          if(l.ne.klev)
     300     s     t_khvd(i,l)=TAUHID(i,klev-l+1,k)
     301     s                -TAUHID(i,klev-l+1-1,k)
     302          if(l.eq.klev)
     303     s     t_khvd(i,l)=TAUHID(i,klev-l+1,k)
     304
     305         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
     306         enddo
     307         enddo
     308         write(str1,'(i2.2)') k
     309      zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
     310      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     311      CALL histwrite(nid_day,"khi"//str1,itau_w,zx_tmp_3d,
     312     .                                   iim*jjmp1*klev,ndex3d)
     313       enddo      ! fin boucle NSPECI
     314c
     315c --------------
     316c ----- OPACITE GAZ
    180317       do k=7,NSPECV,10
    181318         do i=1,klon
    … …  
    191328       enddo      ! fin boucle NSPECV
    192329
     330       do k=8,NSPECI,10
     331         do i=1,klon
     332         do l=1,klev
     333           t_tauhvd(i,l)=TAUGID(i,klev-l+1,k)
     334         enddo
     335         enddo
     336         write(str1,'(i2.2)') k
     337      zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
     338      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     339      CALL histwrite(nid_day,"tgi"//str1,itau_w,zx_tmp_3d,
     340     .                                   iim*jjmp1*klev,ndex3d)
     341       enddo      ! fin boucle NSPECI
     342c
     343c --------------
     344c ----- EXTINCTION GAZ
    193345       do k=7,NSPECV,10
    194346         do i=1,klon
    195347         do l=1,klev
    196          if(l.ne.klev)
    197      s    t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
    198      s    -TAUGVD(i,klev-l+1-1,k)
    199 
     348          if(l.ne.klev)
     349     s     t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
     350     s                -TAUGVD(i,klev-l+1-1,k)
    200351          if(l.eq.klev)
    201      s    t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
     352     s     t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
    202353
    203354         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    … …  
    214365         do i=1,klon
    215366         do l=1,klev
    216            t_tauhvd(i,l)=TAUHID(i,klev-l+1,k)
    217          enddo
    218          enddo
    219          write(str1,'(i2.2)') k
    220       zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
    221       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    222       CALL histwrite(nid_day,"thi"//str1,itau_w,zx_tmp_3d,
     367          if(l.ne.klev)
     368     s     t_khvd(i,l)=TAUGID(i,klev-l+1,k)
     369     s                -TAUGID(i,klev-l+1-1,k)
     370
     371          if(l.eq.klev)
     372     s     t_khvd(i,l)=TAUGID(i,klev-l+1,k)
     373
     374         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
     375         enddo
     376         enddo
     377         write(str1,'(i2.2)') k
     378      zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
     379      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     380      CALL histwrite(nid_day,"kgi"//str1,itau_w,zx_tmp_3d,
    223381     .                                   iim*jjmp1*klev,ndex3d)
    224382       enddo      ! fin boucle NSPECI
    225383
    226        do k=8,NSPECI,10
    227          do i=1,klon
    228          do l=1,klev
    229          if(l.ne.klev)
    230      s    t_khvd(i,l)=TAUHID(i,klev-l+1,k)
    231      s    -TAUHID(i,klev-l+1-1,k)
    232 
    233           if(l.eq.klev)
    234      s    t_khvd(i,l)=TAUHID(i,klev-l+1,k)
    235 
    236          t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    237          enddo
    238          enddo
    239          write(str1,'(i2.2)') k
    240       zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
    241       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    242       CALL histwrite(nid_day,"khi"//str1,itau_w,zx_tmp_3d,
    243      .                                   iim*jjmp1*klev,ndex3d)
    244        enddo      ! fin boucle NSPECI
    245 
    246        do k=8,NSPECI,10
    247          do i=1,klon
    248          do l=1,klev
    249            t_tauhvd(i,l)=TAUGID(i,klev-l+1,k)
    250          enddo
    251          enddo
    252          write(str1,'(i2.2)') k
    253       zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
    254       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    255       CALL histwrite(nid_day,"tgi"//str1,itau_w,zx_tmp_3d,
    256      .                                   iim*jjmp1*klev,ndex3d)
    257        enddo      ! fin boucle NSPECI
    258 
    259        do k=8,NSPECI,10
    260          do i=1,klon
    261          do l=1,klev
    262          if(l.ne.klev)
    263      s    t_khvd(i,l)=TAUGID(i,klev-l+1,k)
    264      s    -TAUGID(i,klev-l+1-1,k)
    265 
    266           if(l.eq.klev)
    267      s    t_khvd(i,l)=TAUGID(i,klev-l+1,k)
    268 
    269          t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    270          enddo
    271          enddo
    272          write(str1,'(i2.2)') k
    273       zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
    274       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    275       CALL histwrite(nid_day,"kgi"//str1,itau_w,zx_tmp_3d,
    276      .                                   iim*jjmp1*klev,ndex3d)
    277        enddo      ! fin boucle NSPECI
    278 
     384c --------------
     385c ----- OPACITE NUAGES (ATTENTION PROXY)
     386         if (clouds.eq.1) then
     387           zx_tmp_fi3d(1:klon,1:klev)=occcld(1:klon,1:klev)
     388           CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     389           CALL histwrite(nid_day,"tcld",itau_w,zx_tmp_3d,
     390     .                                   iim*jjmp1*klev,ndex3d)
     391c --------------
     392c ----- EXTINCTION NUAGES (ATTENTION PROXY)
     393           do i=1,klon
     394             t_kcld(i,klev)=occcld(i,klev)
     395     .       /(zzlev(i,klev+1)-zzlev(i,klev))
     396             do j=klev-1,1,-1
     397               t_kcld(i,j)=(occcld(i,j)-occcld(i,j+1))
     398     .         /(zzlev(i,j+1)-zzlev(i,j))
     399             enddo
     400           enddo
     401           zx_tmp_fi3d(1:klon,1:klev)=t_kcld(1:klon,1:klev)
     402           CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     403           CALL histwrite(nid_day,"kcld",itau_w,zx_tmp_3d,
     404     .                                   iim*jjmp1*klev,ndex3d)
     405        endif 
     406c
    279407      ENDIF !lev_histday.GE.3
    280408c

  • trunk/LMDZ.TITAN/libf/phytitan/write_histins.h

    r110 r175  
    143143     .                                   iim*jjmp1*klev,ndex3d)
    144144c
    145 c 3D adding Tau and k  (31/08/10)
    146 c
     145c --------------
     146c ----- OPACITE BRUME
    147147       do k=7,NSPECV,10
    148148         do i=1,klon
    … …  
    158158       enddo      ! fin boucle NSPECV
    159159
     160       do k=8,NSPECI,10
     161         do i=1,klon
     162         do l=1,klev
     163           t_tauhvd(i,l)=TAUHID(i,klev-l+1,k)
     164         enddo
     165         enddo
     166         write(str1,'(i2.2)') k
     167      zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
     168      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     169      CALL histwrite(nid_ins,"thi"//str1,itau_w,zx_tmp_3d,
     170     .                                   iim*jjmp1*klev,ndex3d)
     171       enddo      ! fin boucle NSPECI
     172c --------------
     173c ----- EXTINCTION BRUME
    160174       do k=7,NSPECV,10
    161175         do i=1,klon
    162176         do l=1,klev
    163          if(l.ne.klev)
    164      s    t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
    165      s    -TAUHVD(i,klev-l+1-1,k)
     177          if(l.ne.klev)
     178     s     t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
     179     s                -TAUHVD(i,klev-l+1-1,k)
    166180
    167181          if(l.eq.klev)
    168      s    t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
     182     s     t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
    169183
    170184         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    … …  
    178192       enddo      ! fin boucle NSPECV
    179193
     194       do k=8,NSPECI,10
     195         do i=1,klon
     196         do l=1,klev
     197          if(l.ne.klev)
     198     s     t_khvd(i,l)=TAUHID(i,klev-l+1,k)
     199     s                -TAUHID(i,klev-l+1-1,k)
     200
     201          if(l.eq.klev)
     202     s     t_khvd(i,l)=TAUHID(i,klev-l+1,k)
     203
     204         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
     205         enddo
     206         enddo
     207         write(str1,'(i2.2)') k
     208      zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
     209      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     210      CALL histwrite(nid_ins,"khi"//str1,itau_w,zx_tmp_3d,
     211     .                                   iim*jjmp1*klev,ndex3d)
     212       enddo      ! fin boucle NSPECI
     213c --------------
     214c ----- OPACITE GAZ
    180215       do k=7,NSPECV,10
    181216         do i=1,klon
    … …  
    191226       enddo      ! fin boucle NSPECV
    192227
     228       do k=8,NSPECI,10
     229         do i=1,klon
     230         do l=1,klev
     231           t_tauhvd(i,l)=TAUGID(i,klev-l+1,k)
     232         enddo
     233         enddo
     234         write(str1,'(i2.2)') k
     235      zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
     236      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     237      CALL histwrite(nid_ins,"tgi"//str1,itau_w,zx_tmp_3d,
     238     .                                   iim*jjmp1*klev,ndex3d)
     239       enddo      ! fin boucle NSPECI
     240c --------------
     241c ----- EXTINCTION GAZ
    193242       do k=7,NSPECV,10
    194243         do i=1,klon
    195244         do l=1,klev
    196          if(l.ne.klev)
    197      s    t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
    198      s    -TAUGVD(i,klev-l+1-1,k)
     245          if(l.ne.klev)
     246     s     t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
     247     s                -TAUGVD(i,klev-l+1-1,k)
    199248
    200249          if(l.eq.klev)
    201      s    t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
     250     s     t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
    202251
    203252         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    … …  
    214263         do i=1,klon
    215264         do l=1,klev
    216            t_tauhvd(i,l)=TAUHID(i,klev-l+1,k)
    217          enddo
    218          enddo
    219          write(str1,'(i2.2)') k
    220       zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
    221       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    222       CALL histwrite(nid_ins,"thi"//str1,itau_w,zx_tmp_3d,
    223      .                                   iim*jjmp1*klev,ndex3d)
    224        enddo      ! fin boucle NSPECI
    225 
    226        do k=8,NSPECI,10
    227          do i=1,klon
    228          do l=1,klev
    229          if(l.ne.klev)
    230      s    t_khvd(i,l)=TAUHID(i,klev-l+1,k)
    231      s    -TAUHID(i,klev-l+1-1,k)
     265          if(l.ne.klev)
     266     s     t_khvd(i,l)=TAUGID(i,klev-l+1,k)
     267     s                -TAUGID(i,klev-l+1-1,k)
    232268
    233269          if(l.eq.klev)
    234      s    t_khvd(i,l)=TAUHID(i,klev-l+1,k)
    235 
    236          t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    237          enddo
    238          enddo
    239          write(str1,'(i2.2)') k
    240       zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
    241       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    242       CALL histwrite(nid_ins,"khi"//str1,itau_w,zx_tmp_3d,
    243      .                                   iim*jjmp1*klev,ndex3d)
    244        enddo      ! fin boucle NSPECI
    245 
    246        do k=8,NSPECI,10
    247          do i=1,klon
    248          do l=1,klev
    249            t_tauhvd(i,l)=TAUGID(i,klev-l+1,k)
    250          enddo
    251          enddo
    252          write(str1,'(i2.2)') k
    253       zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
    254       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    255       CALL histwrite(nid_ins,"tgi"//str1,itau_w,zx_tmp_3d,
    256      .                                   iim*jjmp1*klev,ndex3d)
    257        enddo      ! fin boucle NSPECI
    258 
    259        do k=8,NSPECI,10
    260          do i=1,klon
    261          do l=1,klev
    262          if(l.ne.klev)
    263      s    t_khvd(i,l)=TAUGID(i,klev-l+1,k)
    264      s    -TAUGID(i,klev-l+1-1,k)
    265 
    266           if(l.eq.klev)
    267      s    t_khvd(i,l)=TAUGID(i,klev-l+1,k)
     270     s     t_khvd(i,l)=TAUGID(i,klev-l+1,k)
    268271
    269272         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))

  • trunk/LMDZ.TITAN/libf/phytitan/write_histmth.h

    r119 r175  
    9292      CALL histwrite(nid_mth,"tops",itau_w,zx_tmp_2d,iim*jjmp1,ndex2d)
    9393c
     94cccccccccccccccccc  Tracers
     95c
    9496         if (iflag_trac.eq.1) THEN
    95 c         if (microfi.eq.1) then
     97          if (microfi.ge.1) then
    9698c          DO iq=1,nmicro
    9799c      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qaer(1,1,iq), zx_tmp_3d)
    … …  
    99101c    .                                   iim*jjmp1*klev,ndex3d)
    100102c          ENDDO
    101 c         endif
     103c    -------   NB AER TOT
     104               do i=1,klon
     105                 do j=1,klev
     106                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,1:nrad))
     107                 enddo
     108               enddo
     109       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     110       CALL histwrite(nid_mth,"qaer",itau_w,zx_tmp_3d,
     111     .                                   iim*jjmp1*klev,ndex3d)
     112c    -------   NB NOY TOT
     113               do i=1,klon
     114                 do j=1,klev
     115                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,nrad+1:2*nrad))
     116                 enddo
     117               enddo
     118       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     119       CALL histwrite(nid_mth,"qnoy",itau_w,zx_tmp_3d,
     120     .                                   iim*jjmp1*klev,ndex3d)
     121c    -------   V GLA1 TOT
     122               do i=1,klon
     123                 do j=1,klev
     124                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,2*nrad+1:3*nrad))
     125                 enddo
     126               enddo
     127       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     128       CALL histwrite(nid_mth,"qgl1",itau_w,zx_tmp_3d,
     129     .                                   iim*jjmp1*klev,ndex3d)
     130c    -------   V GLA2 TOT
     131               do i=1,klon
     132                 do j=1,klev
     133                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,3*nrad+1:4*nrad))
     134                 enddo
     135               enddo
     136       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     137       CALL histwrite(nid_mth,"qgl2",itau_w,zx_tmp_3d,
     138     .                                   iim*jjmp1*klev,ndex3d)
     139c    -------   V GLA3 TOT
     140               do i=1,klon
     141                 do j=1,klev
     142                   zx_tmp_fi3d(i,j)= SUM(qaer(i,j,4*nrad+1:5*nrad))
     143                 enddo
     144               enddo
     145       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     146       CALL histwrite(nid_mth,"qgl3",itau_w,zx_tmp_3d,
     147     .                                   iim*jjmp1*klev,ndex3d)
     148c --------------
     149c ----- SATURATION ESP NUAGES
     150             if (clouds.eq.1) then
     151
     152       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, satch4,zx_tmp_3d)
     153       CALL histwrite(nid_mth,"ch4sat", itau_w, zx_tmp_3d,
     154     .                                   iim*jjmp1*klev,ndex3d)
     155
     156       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, satc2h6,zx_tmp_3d)
     157       CALL histwrite(nid_mth,"c2h6sat", itau_w, zx_tmp_3d,
     158     .                                   iim*jjmp1*klev,ndex3d)
     159
     160       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, satc2h2,zx_tmp_3d)
     161       CALL histwrite(nid_mth,"c2h2sat", itau_w, zx_tmp_3d,
     162     .                                   iim*jjmp1*klev,ndex3d)
     163c --------------
     164c ----- RESERVOIR DE SURFACE
     165       CALL gr_fi_ecrit(1, klon,iim,jjmp1,reservoir,zx_tmp_2d)
     166       CALL histwrite(nid_mth,"reserv",itau_w,zx_tmp_2d,
     167     .                        iim*jjmp1,ndex2d)
     168c --------------
     169c ----- PRECIPITATIONS
     170c       -----  CH4
     171       CALL gr_fi_ecrit(1, klon,iim,jjmp1,prec(:,1),zx_tmp_2d)
     172       CALL histwrite(nid_mth,"prech4",itau_w,zx_tmp_2d,
     173     .                        iim*jjmp1,ndex2d)
     174c       -----  C2H6
     175       CALL gr_fi_ecrit(1, klon,iim,jjmp1,prec(:,2),zx_tmp_2d)
     176       CALL histwrite(nid_mth,"prec2h6",itau_w,zx_tmp_2d,
     177     .                        iim*jjmp1,ndex2d)
     178c       -----  C2H2
     179       CALL gr_fi_ecrit(1, klon,iim,jjmp1,prec(:,3),zx_tmp_2d)
     180       CALL histwrite(nid_mth,"prec2h2",itau_w,zx_tmp_2d,
     181     .                        iim*jjmp1,ndex2d)
     182c
     183c --------------
     184c ----- FLUX GLACE
     185       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,flxesp_i(1,1,1),zx_tmp_3d)
     186       CALL histwrite(nid_mth,"flxgl1", itau_w, zx_tmp_3d,
     187     .                                   iim*jjmp1*klev,ndex3d)
     188       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,flxesp_i(1,1,2),zx_tmp_3d)
     189       CALL histwrite(nid_mth,"flxgl2", itau_w, zx_tmp_3d,
     190     .                                   iim*jjmp1*klev,ndex3d)
     191       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,flxesp_i(1,1,3),zx_tmp_3d)
     192       CALL histwrite(nid_mth,"flxgl3", itau_w, zx_tmp_3d,
     193     .                                   iim*jjmp1*klev,ndex3d)
     194c
     195c --------------
     196c ----- RAYON MOYEN GOUTTE
     197       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, rmcloud,zx_tmp_3d)
     198       CALL histwrite(nid_mth,"rcldbar", itau_w, zx_tmp_3d,
     199     .                                   iim*jjmp1*klev,ndex3d)
     200c
     201             endif
     202          endif
     203c
     204c --------------
     205c ----- TRACEURS CHIMIQUES
    102206          if (nmicro.lt.nqmax) then
    103207           DO iq=nmicro+1,nqmax
    … …  
    145249     .                                   iim*jjmp1*klev,ndex3d)
    146250c
    147 c 3D adding Tau and k  (31/08/10)
    148 c
     251c --------------
     252c ----- OPACITE BRUME
    149253       do k=7,NSPECV,10
    150254         do i=1,klon
    … …  
    160264       enddo      ! fin boucle NSPECV
    161265
     266       do k=8,NSPECI,10
     267         do i=1,klon
     268         do l=1,klev
     269           t_tauhvd(i,l)=TAUHID(i,klev-l+1,k)
     270         enddo
     271         enddo
     272         write(str1,'(i2.2)') k
     273      zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
     274      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     275      CALL histwrite(nid_mth,"thi"//str1,itau_w,zx_tmp_3d,
     276     .                                   iim*jjmp1*klev,ndex3d)
     277       enddo      ! fin boucle NSPECI
     278c
     279c --------------
     280c ----- EXTINCTION BRUME
    162281       do k=7,NSPECV,10
    163282         do i=1,klon
    164283         do l=1,klev
    165          if(l.ne.klev)
    166      s    t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
    167      s    -TAUHVD(i,klev-l+1-1,k)
    168 
     284          if(l.ne.klev)
     285     s     t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
     286     s                -TAUHVD(i,klev-l+1-1,k)
    169287          if(l.eq.klev)
    170      s    t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
     288     s     t_khvd(i,l)=TAUHVD(i,klev-l+1,k)
    171289
    172290         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    … …  
    180298       enddo      ! fin boucle NSPECV
    181299
     300       do k=8,NSPECI,10
     301         do i=1,klon
     302         do l=1,klev
     303          if(l.ne.klev)
     304     s     t_khvd(i,l)=TAUHID(i,klev-l+1,k)
     305     s                -TAUHID(i,klev-l+1-1,k)
     306          if(l.eq.klev)
     307     s     t_khvd(i,l)=TAUHID(i,klev-l+1,k)
     308
     309         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
     310         enddo
     311         enddo
     312         write(str1,'(i2.2)') k
     313      zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
     314      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     315      CALL histwrite(nid_mth,"khi"//str1,itau_w,zx_tmp_3d,
     316     .                                   iim*jjmp1*klev,ndex3d)
     317       enddo      ! fin boucle NSPECI
     318c
     319c --------------
     320c ----- OPACITE GAZ
    182321       do k=7,NSPECV,10
    183322         do i=1,klon
    … …  
    193332       enddo      ! fin boucle NSPECV
    194333
     334       do k=8,NSPECI,10
     335         do i=1,klon
     336         do l=1,klev
     337           t_tauhvd(i,l)=TAUGID(i,klev-l+1,k)
     338         enddo
     339         enddo
     340         write(str1,'(i2.2)') k
     341      zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
     342      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     343      CALL histwrite(nid_mth,"tgi"//str1,itau_w,zx_tmp_3d,
     344     .                                   iim*jjmp1*klev,ndex3d)
     345       enddo      ! fin boucle NSPECI
     346c
     347c --------------
     348c ----- EXTINCTION GAZ
    195349       do k=7,NSPECV,10
    196350         do i=1,klon
    197351         do l=1,klev
    198          if(l.ne.klev)
    199      s    t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
    200      s    -TAUGVD(i,klev-l+1-1,k)
    201 
     352          if(l.ne.klev)
     353     s     t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
     354     s                -TAUGVD(i,klev-l+1-1,k)
    202355          if(l.eq.klev)
    203      s    t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
     356     s     t_khvd(i,l)=TAUGVD(i,klev-l+1,k)
    204357
    205358         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    … …  
    216369         do i=1,klon
    217370         do l=1,klev
    218            t_tauhvd(i,l)=TAUHID(i,klev-l+1,k)
    219          enddo
    220          enddo
    221          write(str1,'(i2.2)') k
    222       zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
    223       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    224       CALL histwrite(nid_mth,"thi"//str1,itau_w,zx_tmp_3d,
     371          if(l.ne.klev)
     372     s     t_khvd(i,l)=TAUGID(i,klev-l+1,k)
     373     s                -TAUGID(i,klev-l+1-1,k)
     374
     375          if(l.eq.klev)
     376     s     t_khvd(i,l)=TAUGID(i,klev-l+1,k)
     377
     378         t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
     379         enddo
     380         enddo
     381         write(str1,'(i2.2)') k
     382      zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
     383      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     384      CALL histwrite(nid_mth,"kgi"//str1,itau_w,zx_tmp_3d,
    225385     .                                   iim*jjmp1*klev,ndex3d)
    226386       enddo      ! fin boucle NSPECI
    227387
    228        do k=8,NSPECI,10
    229          do i=1,klon
    230          do l=1,klev
    231          if(l.ne.klev)
    232      s    t_khvd(i,l)=TAUHID(i,klev-l+1,k)
    233      s    -TAUHID(i,klev-l+1-1,k)
    234 
    235           if(l.eq.klev)
    236      s    t_khvd(i,l)=TAUHID(i,klev-l+1,k)
    237 
    238          t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    239          enddo
    240          enddo
    241          write(str1,'(i2.2)') k
    242       zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
    243       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    244       CALL histwrite(nid_mth,"khi"//str1,itau_w,zx_tmp_3d,
    245      .                                   iim*jjmp1*klev,ndex3d)
    246        enddo      ! fin boucle NSPECI
    247 
    248        do k=8,NSPECI,10
    249          do i=1,klon
    250          do l=1,klev
    251            t_tauhvd(i,l)=TAUGID(i,klev-l+1,k)
    252          enddo
    253          enddo
    254          write(str1,'(i2.2)') k
    255       zx_tmp_fi3d(1:klon,1:klev)=t_tauhvd(1:klon,1:klev)
    256       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    257       CALL histwrite(nid_mth,"tgi"//str1,itau_w,zx_tmp_3d,
    258      .                                   iim*jjmp1*klev,ndex3d)
    259        enddo      ! fin boucle NSPECI
    260 
    261        do k=8,NSPECI,10
    262          do i=1,klon
    263          do l=1,klev
    264          if(l.ne.klev)
    265      s    t_khvd(i,l)=TAUGID(i,klev-l+1,k)
    266      s    -TAUGID(i,klev-l+1-1,k)
    267 
    268           if(l.eq.klev)
    269      s    t_khvd(i,l)=TAUGID(i,klev-l+1,k)
    270 
    271          t_khvd(i,l)=t_khvd(i,l)/(zzlev(i,l+1)-zzlev(i,l))
    272          enddo
    273          enddo
    274          write(str1,'(i2.2)') k
    275       zx_tmp_fi3d(1:klon,1:klev)=t_khvd(1:klon,1:klev)
    276       CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    277       CALL histwrite(nid_mth,"kgi"//str1,itau_w,zx_tmp_3d,
    278      .                                   iim*jjmp1*klev,ndex3d)
    279        enddo      ! fin boucle NSPECI
    280 
     388c --------------
     389c ----- OPACITE NUAGES (ATTENTION PROXY)
     390         if (clouds.eq.1) then
     391           zx_tmp_fi3d(1:klon,1:klev)=occcld(1:klon,1:klev)
     392           CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     393           CALL histwrite(nid_mth,"tcld",itau_w,zx_tmp_3d,
     394     .                                   iim*jjmp1*klev,ndex3d)
     395c --------------
     396c ----- EXTINCTION NUAGES (ATTENTION PROXY)
     397           do i=1,klon
     398             t_kcld(i,klev)=occcld(i,klev)
     399     .       /(zzlev(i,klev+1)-zzlev(i,klev))
     400             do j=klev-1,1,-1
     401               t_kcld(i,j)=(occcld(i,j)-occcld(i,j+1))
     402     .         /(zzlev(i,j+1)-zzlev(i,j))
     403             enddo
     404           enddo
     405           zx_tmp_fi3d(1:klon,1:klev)=t_kcld(1:klon,1:klev)
     406           CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     407           CALL histwrite(nid_mth,"kcld",itau_w,zx_tmp_3d,
     408     .                                   iim*jjmp1*klev,ndex3d)
     409c
     410c --------------
     411c ----- OCCURENCE NUAGES
     412           do k=1,12
     413             write(str1,'(i2.2)') k
     414             zx_tmp_fi3d(1:klon,1:klev)=occcld_m(1:klon,1:klev,k)
     415             CALL histwrite(nid_mth,"occcld"//str1,itau_w,zx_tmp_3d,
     416     .                                   iim*jjmp1*klev,ndex3d)
     417           enddo
     418c
     419        endif 
     420c
    281421      ENDIF !lev_histmth.GE.3
    282422c
Note: See TracChangeset for help on using the changeset viewer.