Changeset 1591 for trunk/LMDZ.VENUS/libf


Ignore:
Timestamp:
Aug 31, 2016, 4:31:29 PM (8 years ago)
Author:
slebonnois
Message:

SL: Mise a jour de la haute atmosphere, du transfert radiatif (solaire=Rainer Haus; IR: reglages de juin 2016), et implementation des variations de temperature potentielle dans la basse atmosphere (variation de la masse molaire moyenne)

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.VENUS/libf/phyvenus/clmain.F

    r1543 r1591  
    291291c        call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,ycoefm
    292292c    .   ,ycoefm0,ycoefh0)
    293 c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')
    294 c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')
     293c      call dump2d(nbp_lon,nbp_lat-2,ycoefm(2:klon-1,2), 'KZ         ')
     294c      call dump2d(nbp_lon,nbp_lat-2,ycoefm0(2:klon-1,2),'KZMIN      ')
    295295 
    296296         ycoefm0 = 1.e-3
    … …  
    349349        ENDIF
    350350
    351 c   iflag_pbl peut etre utilise comme longuer de melange
     351c   iflag_pbl peut etre utilise comme longueur de melange
    352352
    353353         if (iflag_pbl.ge.11) then
    … …  
    558558      PARAMETER (check=.false.)
    559559C
     560c----------------------
     561c ATMOSPHERE PROFONDE
     562      real p_lim,p_ref
     563      real rsmu_mod(klon,klev),zrsmu_mod(klon,klev)
     564      real ratio_mod(klon,klev)
     565! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
     566
     567      p_lim = 6e6
     568      p_ref = 8.9e6
     569
     570      DO k = 1, klev
     571      DO i = 1, klon
     572        ratio_mod(i,k) = 1.
     573        rsmu_mod(i,k) = RD
     574       if (pplay(i,k).gt.p_lim) then
     575          ratio_mod(i,k) = RMD /
     576     &  (RMD+0.56*(log(pplay(i,k))-log(p_lim))/(log(p_ref)-log(p_lim)))
     577           ratio_mod(i,k) = max(ratio_mod(i,k),RMD/(RMD+0.56))
     578           rsmu_mod(i,k) = RD*ratio_mod(i,k)
     579       endif
     580      ENDDO
     581      ENDDO
     582c----------------------
     583
    560584      if (iflag_pbl.eq.1) then
    561585        do k = 3, klev
    … …  
    584608      DO i = 1, knon
    585609            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
     610c----------------------
     611c ATMOSPHERE PROFONDE
     612            zrsmu_mod(i,k)    = (rsmu_mod(i,k)+rsmu_mod(i,k-1)) * 0.5
     613c----------------------
    586614      ENDDO
    587615      ENDDO
    … …  
    612640      DO i = 1, knon
    613641         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
    614      .                  *(paprs(i,k)/zt(i,k)/RD)**2
     642c----------------------
     643c ATMOSPHERE PROFONDE
     644     .                  *(paprs(i,k)/zt(i,k)/zrsmu_mod(i,k))**2
     645c----------------------
    615646         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
    616647      ENDDO
    617648      ENDDO
    618649c
    619 c Preparer les flux lies aux contre-gardients
     650c Preparer les flux lies aux contre-gardients  OBSOLETE
    620651c
    621652      DO k = 2, klev
    … …  
    9841015      DATA appel1er /.TRUE./
    9851016c
     1017c----------------------
     1018c ATMOSPHERE PROFONDE
     1019      real p_lim,p_ref
     1020      real modgam(klon,klev)
     1021
     1022      p_lim = 6e6
     1023      p_ref = 8.9e6
     1024
     1025      DO k = 1, klev
     1026      DO i = 1, klon
     1027        modgam(i,k) = 0.
     1028       if (pplay(i,k).gt.p_lim) then
     1029           modgam(i,k) = 0.56/(log(p_ref)-log(p_lim))/R
     1030       endif
     1031      ENDDO
     1032      ENDDO
     1033c----------------------
     1034
    9861035      isommet=klev
    9871036
    … …  
    10451094            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
    10461095            zdphi      = zmgeom(i,k)/2.
    1047             ztvd(i,k)  = (t(i,k)   + zdphi/cpdet(zt(i,k)))
    1048             ztvu(i,k)  = (t(i,k-1) - zdphi/cpdet(zt(i,k)))
     1096c----------------------
     1097c ATMOSPHERE PROFONDE
     1098            ztvd(i,k)  = t(i,k)   + zdphi*
     1099     &                      (1./cpdet(zt(i,k))+modgam(i,k))
     1100            ztvu(i,k)  = t(i,k-1) - zdphi*
     1101     &                      (1./cpdet(zt(i,k))+modgam(i,k))
     1102c----------------------
    10491103            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
    10501104      ENDDO
    … …  
    10541108        zt(i,1)   = ts(i)
    10551109        ztvu(i,1) = ts(i)
    1056         ztvd(i,1) = t(i,1)+zgeop(i,1)/cpdet(zt(i,1))
     1110c----------------------
     1111c ATMOSPHERE PROFONDE
     1112        ztvd(i,1) = t(i,1)+zgeop(i,1)*
     1113     &                      (1./cpdet(zt(i,1))+modgam(i,1))
     1114c----------------------
    10571115        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
    10581116      ENDDO

  • trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F

    r1530 r1591  
    174174c         endif
    175175
    176   !     n2
    177 c            akin2 = 5.6e-4
    178 c            cpin2 = 1.034e3
    179 
    180 
    181          
    182176         ! tell the world about it:
    183177         write(*,*) "concentrations: firstcall, nbq=",nbq
    … …  
    224218            end do
    225219            mmean(ig,l) = 1./mmean(ig,l)
    226             rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg K           
    227  
    228 c            write(*,*),'Mmean: ',ig, l, mmean(ig,l)
     220            rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg K
     221c            write(*,*),'Mmean in concentration2: ',ig, l, mmean(ig,l)
    229222         end do
    230223      end do
    … …  
    244237c            write(*,*),'Air density: ',ig, l, rho(0,l)           
    245238
    246 !!! --- INSERT N2 values ----
    247239!!  WARNING -> Cp here below doesn't depend on T (cpdet)
    248240

  • trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h

    r1442 r1591  
    99      real*8 dzres
    1010
    11       parameter (Pdiff=15.)      ! pressure below which diffusion is computed
    12       parameter (tdiffmin=5d0)
    13       parameter (dzres=2d0)      ! grid resolution (km) for diffusion
     11      parameter (Pdiff=1.e-1)      ! pressure below which diffusion is computed
     12      parameter (tdiffmin=5.d0)
     13      parameter (dzres=2.d0)     ! grid resolution (km) for diffusion
    1414     

  • trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90

    r1530 r1591  
    55    ! Parametrization of the momentum flux deposition due to a discrete
    66    ! number of gravity waves.
    7     ! F. Lott (version 9: 16 February, 2012), reproduce v3 but with only
    8     !         two waves present at each time step
     7    ! F. Lott
     8    ! Version 14, Gaussian distribution of the source
    99    ! LMDz model online version     
    10     ! ADAPTED FOR VENUS
     10    ! ADAPTED FOR VENUS /  F. LOTT + S. LEBONNOIS
     11    ! Version adapted on 03/04/2013:
     12    !      - input flux compensated in the deepest layers
    1113    !---------------------------------------------------------------------
    1214
    1315      use dimphy
     16      use moyzon_mod, only: plevmoy
    1417      implicit none
    1518
    … …  
    3942    ! O.3 INTERNAL ARRAYS
    4043
    41     INTEGER II, LL, IEQ
     44    INTEGER II, LL
    4245
    4346    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
    … …  
    6265    REAL ZOM(NW, KLON), ZOP(NW, KLON)
    6366
    64     ! Wave vertical velocities at the 2 1/2 lev surrounding the full level
     67    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
    6568    REAL WWM(NW, KLON), WWP(NW, KLON)
    66 
    67     REAL RUW0(NW, KLON) ! Fluxes at launching level
    68 
     69    ! Fluxes X and Y for each waves at 1/2 Levels
    6970    REAL RUWP(NW, KLON), RVWP(NW, KLON)
    70     ! Fluxes X and Y for each waves at 1/2 Levels
    71 
    72     INTEGER LAUNCH   ! Launching altitude
    73 
    74     REAL RUWMAX,SAT  ! saturation parameter
    75     REAL XLAUNCH     ! Controle the launching altitude
    7671    REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
    7772    REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
    7873
     74    REAL RUW0(NW, KLON) ! Fluxes at launching level
     75    REAL RUWMAX         ! Max value
     76    INTEGER LAUNCH      ! Launching altitude
     77    REAL XLAUNCH        ! Control the launching altitude
     78
    7979    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
    8080
     81    REAL SAT  ! saturation parameter
    8182    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
    8283
    … …  
    99100    logical output
    100101    data output/.false./
    101   ! CAUTION ! IF output is .true. THEN change NO to 10 at least !
     102!    data output/.true./
     103 ! CAUTION ! IF output is .true. THEN change NO to 10 at least !
    102104    character*14 outform
    103105    character*2  str2
     106    integer      ieq
    104107
    105108! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
    … …  
    127130    RUWMAX = 0.005      ! Max EP-Flux at Launch altitude
    128131    SAT    = 0.85       ! Saturation parameter: Sc in (12)
    129     RDISS  = 10.        ! Diffusion parameter
     132    RDISS  = 0.1        ! Diffusion parameter
    130133
    131134    DELTAT=24.*3600.    ! Time scale of the waves (first introduced in 9b)
    132135
    133     KMIN = 1.E-6        ! Min horizontal wavenumber
    134     KMAX = 2.E-5        ! Max horizontal wavenumber
     136!!!! TEST GG. Values corresponding to min/max horizontal wavel 50-500 km
     137!(similar to observations)
     138    KMIN = 1.E-5        ! Min horizontal wavenumber
     139    KMAX = 1.E-4        ! Max horizontal wavenumber
    135140    !Online output: one value only
    136141    if (output) then
    137       KMIN = 6.3E-6
    138       KMAX = 6.3E-6
     142      KMIN = 3E-5
     143      KMAX = 3E-5
    139144    endif
    140145    CMIN = 1.           ! Min phase velocity
    141146    CMAX = 61.          ! Max phase speed velocity
    142     XLAUNCH=0.6         ! Parameter that control launching altitude
    143 
    144     PR = 9.2e6          ! Reference pressure    ! VENUS!!
    145     TR = 240.           ! Reference Temperature ! VENUS: cloud layer
     147!    XLAUNCH=0.6         ! Parameter that control launching altitude
     148    XLAUNCH=5e-3        ! Value for top of cloud convective region
     149
     150!    PR = 9.2e6          ! Reference pressure    ! VENUS!!
     151    PR = 5e5            ! Reference pressure    ! VENUS: cloud layer
     152    TR = 300.           ! Reference Temperature ! VENUS: cloud layer
    146153    H0 = RD * TR / RG   ! Characteristic vertical scale height
    147154
    … …  
    162169    ENDIF
    163170
    164     ! 1.2 WAVES CHARACTERISTICS CHOSEN RANDOMLY
     171    !  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
     172    !-------------------------------------------------------------
     173
     174    !Online output
     175    if (output) OPEN(11,file="impact-gwno.dat")
     176
     177    ! Pressure and Inv of pressure, Temperature / at 1/2 level
     178    DO LL = 2, KLEV
     179       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
     180    end DO
     181
     182    PH(:, KLEV + 1) = 0.
     183    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
     184
     185    ! Launching altitude
     186
     187    DO LL = 1, NLEV
     188       IF (plevmoy(LL) / plevmoy(1) > XLAUNCH) LAUNCH = LL
     189    ENDDO
     190! test
     191!    print*,"launch=",LAUNCH
     192!    print*,"launch p,N2=",plevmoy(LAUNCH),pn2(nlon/2+1,LAUNCH)
     193
     194    ! Log pressure vert. coordinate
     195    DO LL = 1, KLEV + 1
     196       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
     197    end DO
     198
     199    if (output) then
     200    ! altitude above surface
     201       ZHbis(:,1) = 0.   
     202       DO LL = 2, KLEV + 1
     203          H0bis(:, LL-1) = RD * TT(:, LL-1) / RG
     204          ZHbis(:, LL) = ZHbis(:, LL-1) &
     205           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
     206       end DO
     207    endif
     208
     209    ! Winds and BV frequency
     210    DO LL = 2, KLEV
     211       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
     212       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
     213       ! BVSEC: BV Frequency
     214! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
     215       BV(:, LL) = MAX(BVSEC,SQRT(pn2(:,LL)))
     216    end DO
     217    BV(:, 1) = BV(:, 2)
     218    UH(:, 1) = 0.
     219    VH(:, 1) = 0.
     220    BV(:, KLEV + 1) = BV(:, KLEV)
     221    UH(:, KLEV + 1) = UU(:, KLEV)
     222    VH(:, KLEV + 1) = VV(:, KLEV)
     223
     224
     225    ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
    165226    !-------------------------------------------
    166227
    … …  
    168229    ! are used to produce the waves characteristics
    169230    ! in a stochastic way
     231
     232!! A REVOIR:
     233!! - utilisation de MOD ou bien de aleas ?
     234!! - distribution gaussienne des CPHA ? (avec signe ZP qui est ajuste apres)
    170235
    171236    JW = 0
    … …  
    187252                ! Intrinsic frequency
    188253                ZO(JW, II) = CPHA * ZK(JW, II)
     254                ! Intrinsic frequency  is imposed
     255                    ZO(JW, II) = ZO(JW, II)      &
     256                  + ZK(JW, II) * COS(ZP(JW)) * UH(II, LAUNCH) &
     257                  + ZK(JW, II) * SIN(ZP(JW)) * VH(II, LAUNCH)
    189258                ! Momentum flux at launch lev
    190259                ! RUW0(JW, II) = RUWMAX / REAL(NW) &
    … …  
    199268    end DO
    200269
    201     !  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
    202     !-------------------------------------------------------------
    203 
    204     IEQ = KLON / 2
    205     !Online output
    206     if (output) OPEN(11,file="impact-gwno.dat")
    207 
    208     ! Pressure and Inv of pressure, Temperature / at 1/2 level
    209     DO LL = 2, KLEV
    210        PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
    211     end DO
    212 
    213     PH(:, KLEV + 1) = 0.
    214     PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
    215 
    216     ! Launching altitude
    217 
    218     DO LL = 1, NLEV
    219        IF (PH(IEQ, LL) / PH(IEQ, 1) > XLAUNCH) LAUNCH = LL
    220     ENDDO
    221 
    222     ! Log pressure vert. coordinate (altitude above surface)
    223     ZHbis(:,1) = 0.   
    224     DO LL = 2, KLEV + 1
    225        H0bis(:, LL-1) = RD * TT(:, LL-1) / RG
    226        ZHbis(:, LL) = ZHbis(:, LL-1) &
    227            + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
    228     end DO
    229     ! Log pressure vert. coordinate
    230     DO LL = 1, KLEV + 1
    231        ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
    232     end DO
    233 
    234 
    235     ! Winds and BV frequency
    236     DO LL = 2, KLEV
    237        UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
    238        VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
    239        ! BVSEC: BV Frequency
    240 ! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
    241        BV(:, LL) = MAX(BVSEC,SQRT(pn2(:,LL)))
    242     end DO
    243     BV(:, 1) = BV(:, 2)
    244     UH(:, 1) = 0.
    245     VH(:, 1) = 0.
    246     BV(:, KLEV + 1) = BV(:, KLEV)
    247     UH(:, KLEV + 1) = UU(:, KLEV)
    248     VH(:, KLEV + 1) = VV(:, KLEV)
    249 
    250 
    251     ! 3. COMPUTE THE FLUXES
     270    ! 4. COMPUTE THE FLUXES
    252271    !--------------------------
    253272
    254     !  3.1  Vertical velocity at launching altitude to ensure
     273    !  4.1  Vertical velocity at launching altitude to ensure
    255274    !       the correct value to the imposed fluxes.
    256275    !
    … …  
    261280            - ZK(JW, :) * COS(ZP(JW)) * UH(:, LAUNCH) &
    262281            - ZK(JW, :) * SIN(ZP(JW)) * VH(:, LAUNCH)
    263        ! Vertical velocity at launch level, value to ensure the imposed
    264        ! mom flux:
    265        WWP(JW, :) = SQRT(ABS(ZOP(JW, :)) / MAX(BV(:, LAUNCH),BVSEC) &
    266             * RUW0(JW,:))
     282       ! WW is directly a flux, here, not vertical velocity anymore
     283       WWP(JW, :) = RUW0(JW,:)
    267284       RUWP(JW, :) = COS(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
    268285       RVWP(JW, :) = SIN(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
    … …  
    270287    end DO
    271288
    272     !  3.2 Uniform values below the launching altitude
    273 
    274     DO LL = 1, LAUNCH
    275        RUW(:, LL) = 0
    276        RVW(:, LL) = 0
    277        DO JW = 1, NW
    278           RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
    279           RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
    280        end DO
    281     end DO
    282 
    283     !  3.3 Loop over altitudes, with passage from one level to the
     289    !  4.2 Initial flux at launching altitude
     290
     291    RUW(:, LAUNCH) = 0
     292    RVW(:, LAUNCH) = 0
     293    DO JW = 1, NW
     294       RUW(:, LAUNCH) = RUW(:, LAUNCH) + RUWP(JW, :)
     295       RVW(:, LAUNCH) = RVW(:, LAUNCH) + RVWP(JW, :)
     296    end DO
     297
     298    !  4.3 Loop over altitudes, with passage from one level to the
    284299    !      next done by i) conserving the EP flux, ii) dissipating
    285300    !      a little, iii) testing critical levels, and vi) testing
    … …  
    287302
    288303    !Online output
    289     write(str2,'(i2)') NW+2
    290     outform="("//str2//"(E12.4,1X))"
    291     if (output) WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
     304    if (output) then
     305        ieq=nlon/2+1
     306        write(str2,'(i2)') NW+2
     307        outform="("//str2//"(E12.4,1X))"
     308        WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
    292309               (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW)), JW = 1, NW)
     310    endif
    293311
    294312    DO LL = LAUNCH, KLEV - 1
    … …  
    305323
    306324          WWP(JW, :) = MIN( &
    307                ! No breaking (Eq.6)
     325       ! No breaking (Eq.6)
    308326               WWM(JW, :) &
    309                * SQRT(BV(:, LL ) / BV(:, LL+1) &
    310                * ABS(ZOP(JW, :)) / MAX(ABS(ZOM(JW, :)), ZOISEC)) &
    311                ! Dissipation (Eq. 8):
     327       ! Dissipation (Eq. 8):
    312328               * EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) &
    313329               * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
    314330               / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
    315331               * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), &
    316                ! Critical levels (forced to zero if intrinsic
    317                ! frequency changes sign)
     332       ! Critical levels (forced to zero if intrinsic frequency changes sign)
    318333               MAX(0., SIGN(1., ZOP(JW, :) * ZOM(JW, :))) &
    319                ! Saturation (Eq. 12)
    320                * ZOP(JW, :)**2 / ZK(JW, :)/BV(:, LL+1) &
    321                * EXP(-ZH(:, LL + 1)/2./H0) * SAT)
     334       ! Saturation (Eq. 12)
     335               * ABS(ZOP(JW, :))**3 /BV(:, LL+1) &
     336               * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4) 
    322337       end DO
    323338
    … …  
    327342
    328343       DO JW = 1, NW
    329           RUWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
    330                *BV(:,LL+1)&
    331                * COS(ZP(JW)) *  WWP(JW, :)**2
    332           RVWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
    333                *BV(:,LL+1)&
    334                * SIN(ZP(JW)) *  WWP(JW, :)**2
     344          RUWP(JW, :) = SIGN(1.,ZOP(JW, :))*COS(ZP(JW))*WWP(JW, :)
     345          RVWP(JW, :) = SIGN(1.,ZOP(JW, :))*SIN(ZP(JW))*WWP(JW, :)
    335346       end DO
    336347       !
    … …  
    358369    end DO
    359370
     371    ! 5 CALCUL DES TENDANCES:
     372    !------------------------
     373
     374    ! 5.1 Rectification des flux au sommet et dans les basses couches:
     375! MODIF SL
     376 
     377! Attention, ici c'est le total sur toutes les ondes...
     378
     379    RUW(:, KLEV + 1) = 0.
     380    RVW(:, KLEV + 1) = 0.
     381
     382    ! Here, big change compared to FLott version:
     383    ! We compensate (RUW(:, LAUNCH), ie total emitted upward flux
     384    !  over the layers max(1,LAUNCH-3) to LAUNCH-1
     385    DO LL = 1, max(1,LAUNCH-3)
     386      RUW(:, LL) = 0.
     387      RVW(:, LL) = 0.
     388    end DO
     389    DO LL = max(2,LAUNCH-2), LAUNCH-1
     390       RUW(:, LL) = RUW(:, LL - 1) + RUW(:, LAUNCH) * &
     391            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     392       RVW(:, LL) = RVW(:, LL - 1) + RVW(:, LAUNCH) * &
     393            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     394    end DO
     395    ! This way, the total flux from GW is zero, but there is a net transport
     396    ! (upward) that should be compensated by circulation
     397    ! and induce additional friction at the surface
     398
    360399    !Online output
    361400    if (output) then
     401       DO LL = 1, KLEV - 1
     402           WRITE(11,*) ZHbis(IEQ, LL)/1000.,RUW(IEQ,LL)
     403       end DO
    362404       CLOSE(11)
    363405       stop
    364406    endif
    365 
    366     ! 4 CALCUL DES TENDANCES:
    367     !------------------------
    368 
    369     ! 4.1 Rectification des flux au sommet et dans les basses couches:
    370 
    371     RUW(:, KLEV + 1) = 0.
    372     RVW(:, KLEV + 1) = 0.
    373     RUW(:, 1) = RUW(:, LAUNCH)
    374     RVW(:, 1) = RVW(:, LAUNCH)
    375     DO LL = 2, LAUNCH
    376        RUW(:, LL) = RUW(:, LL - 1) + (RUW(:, LAUNCH + 1) - RUW(:, 1)) * &
    377             (PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
    378        RVW(:, LL) = RVW(:, LL - 1) + (RVW(:, LAUNCH + 1) - RVW(:, 1)) * &
    379             (PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
    380     end DO
    381407
    382408    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4

  • trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90

    r1530 r1591  
    226226!!!
    227227!! =======> Gab 29 oct 2014
    228 !!$$$$$$$   THIS UPDATE IS ALREADY DONE in physique (t_seri & tr_seri) $$$$$$
     228!!$$$   THIS UPDATE IS ALREADY DONE in physique (t_seri & tr_seri) $$$$$$
    229229
    230230! Update the temperature modified by other processes
    … …  
    449449!        Tdiff=Tdiff/2.
    450450
    451 !       if (ig .eq. ij0) then
    452 !       print*,'test',ig,Tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
    453 !       endif
    454         if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
    455         ntime=anint(dble(ptimestep)/Tdiff)
    456 !       print*,'ptime',ig,ptimestep,Tdiff,ntime,tdiffmin,Mraf(nlraf)
     451        if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
     452
     453!     if (tdiff .lt. tdiffmin*Mraf(nlraf)) then       
     454!        print*, 'Moldiff L454', tdiff, ptimestep, tdiffmin*Mraf(nlraf)
     455!        endif
     456
     457!      JYC + GG aout 2015 add this condition below 
     458        if (tdiff .ge. ptimestep) tdiff=ptimestep
     459        ! Number of time step
     460      ntime=anint(dble(ptimestep)/tdiff)
     461      !print*,'ptime',ig,tdiff,ntime,Mraf(nlraf)
     462
    457463! Adimensionned temperature
    458464
    … …  
    592598        CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
    593599     &   pp,mol_tr,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
     600
     601!!!  Call added by JY+GG Aout 2014
     602        CALL MMOY(massemoy,mol_tr,qnew,gcmind,nlayer,ncompdiff)
    594603
    595604        CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
    … …  
    9941003        SUBROUTINE ZVERT(P,T,M,Z,nl,ig)
    9951004        IMPLICIT NONE
     1005#include "YOMCST.h"
    9961006        INTEGER :: nl,l,ig
    9971007        REAL*8,dimension(nl) :: P,T,M,Z,H
    9981008        REAL*8 :: z0
    9991009        REAL*8 :: kbolt,masseU,Konst,g,Hpm
    1000         masseU=1.660538782d-27
    1001         kbolt=1.3806504d-23
     1010        masseU=1.e-3/RNAVO
     1011        kbolt=RKBOL
    10021012        Konst=kbolt/masseU
    1003         g=3.72D0
     1013        g=RG
    10041014
    10051015        z0=0d0
    … …  
    10251035        SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
    10261036        IMPLICIT NONE
    1027 
     1037#include "YOMCST.h"
    10281038        REAL*8 :: kbolt,masseU,Konst
    10291039        INTEGER :: nl,nq,l,iq
    10301040        REAL*8,DIMENSION(nl) :: P,T,M,RHON
    10311041        REAL*8,DIMENSION(nl,nq) :: RHOK,qq
    1032         masseU=1.660538782d-27
    1033         kbolt=1.3806504d-23
     1042        masseU=1.e-3/RNAVO
     1043        kbolt=RKBOL
    10341044        Konst=Kbolt/masseU
    10351045
    … …  
    10491059        use infotrac
    10501060        IMPLICIT NONE
    1051        
     1061#include "YOMCST.h"
    10521062        INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
    10531063        INTEGER :: gc(nq)
    … …  
    10601070        REAL*8 :: facZ,dZ,H
    10611071        REAL,DIMENSION(nq) :: mol_tr
    1062         masseU=1.660538782d-27
    1063         kbolt=1.3806504d-23
     1072        masseU=1.e-3/RNAVO
     1073        kbolt=RKBOL
    10641074        Konst=Kbolt/masseU
    1065         g=3.72d0
     1075        g=RG
    10661076
    10671077
    … …  
    12061216        SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
    12071217        IMPLICIT NONE
     1218#include "YOMCST.h"
    12081219        INTEGER :: l,nl,nn
    12091220        REAL*8,DIMENSION(nl) :: T,H,D,W,DT
    12101221        REAL*8 :: Dz,Hmol,masse
    12111222        REAL*8 :: kbolt,masseU,Konst,g
    1212         masseU=1.660538782d-27
    1213         kbolt=1.3806504d-23
     1223        masseU=1.e-3/RNAVO
     1224        kbolt=RKBOL
    12141225        Konst=Kbolt/masseU
    1215         g=3.72d0
     1226        g=RG
    12161227
    12171228        DT(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)
    … …  
    13661377        SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
    13671378        IMPLICIT NONE
     1379#include "YOMCST.h"
    13681380        INTEGER :: nl,nn,l,nq,il
    13691381        REAL,DIMENSION(nl+1) :: P
    13701382        REAL*8,DIMENSION(nl,nq) :: qold,qnew
    13711383        REAL*8 :: DM,Mold,Mnew,g
    1372         g=3.72d0
     1384        g=RG
    13731385        DM=0d0
    13741386        Mold=0d0
    … …  
    13811393        ENDDO
    13821394        IF (abs(DM/Mold) .gt. 1d-5) THEN
    1383         Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
     1395        Print*,'We dont conserve mass',nn,DM,Mold,Mnew,DM/Mold
    13841396        ENDIF
    13851397
    … …  
    13911403        use infotrac
    13921404        IMPLICIT NONE
     1405#include "YOMCST.h"
    13931406        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
    13941407        INTEGER,DIMENSION(1) :: indP
    … …  
    14031416        REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
    14041417        REAL*8 :: Znew,Znew2,Pnew,Pnew2
    1405         masseU=1.660538782d-27
    1406         kbolt=1.3806504d-23
     1418        masseU=1.e-3/RNAVO
     1419        kbolt=RKBOL
    14071420        Konst=Kbolt/masseU
    1408         g=3.72d0
     1421        g=RG
    14091422        Dz=Z(2)-Z(1)
    14101423        Znew=Z(nl)
    … …  
    14931506!        use infotrac
    14941507        IMPLICIT NONE
     1508#include "YOMCST.h"
    14951509        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
    14961510        INTEGER,DIMENSION(1) :: indP
    … …  
    15041518        REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
    15051519        REAL*8 :: Znew,Znew2,Pnew,Pnew2
    1506         masseU=1.660538782d-27
    1507         kbolt=1.3806504d-23
     1520        masseU=1.e-3/RNAVO
     1521        kbolt=RKBOL
    15081522        Konst=Kbolt/masseU
    1509         g=3.72d0
     1523        g=RG
    15101524        Dz=Z(2)-Z(1)
    15111525        Znew=Z(nl)

  • trunk/LMDZ.VENUS/libf/phyvenus/nlte_tcool.F

    r1530 r1591  
    8282            cpnew_ig(l)=cpnew(ig,l)
    8383
     84c            write (*,*) 'nlte_tcool z_ig', l, z_ig(l)   
    8485         enddo
    8586
    86                                 ! From GCM's grid to NLTE's grid
     87c     ! From GCM's grid to NLTE's grid
    8788         call NLTEdlvr11_ZGRID_02 (nlev,
    8889     $        p_ig, t_ig, z_ig,
    … …  
    302303      zmin = z_gcm(jlowerboundary)
    303304      zmax = z_gcm(jtopboundary)
    304 !      print*, zmin, zmax
     305
     306!      print*, 'nlte boundaries', zmin, zmax
     307!      print*, 'Top atmosphere', z_gcm(nlev-1)
    305308!      stop
    306309      deltaz = (zmax-zmin) / (nl-1)

  • trunk/LMDZ.VENUS/libf/phyvenus/nlthermeq.F

    r1530 r1591  
    6464      END IF
    6565
    66 c    write(*,*) 'LTE rad. calculations up to layer ',  nlaylte
    67 
    6866c
    6967 100    return
    … …  
    7371   20 format('               half-width (scale heights) ',f6.2)
    7472   30 format('          suggested LTE coverage at least ',f6.2,'Pa')
    75    40 format(' nlthermeq: purely NLTE contribution over (nlayer) ',f6.2) 
     73   40 format(' nlthermeq: purely NLTE contribution over (nlayer) ',f6.4)
    7674
    7775      end

  • trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F

    r1545 r1591  
    5454      REAL   PPB(klev+1)
    5555
    56       REAL   zfract, zrmu0
     56      REAL   zfract, zrmu0,latdeg
    5757
    5858      REAL   zheat(klev), zcool(klev)
    … …  
    273273c---------
    274274      znivs=zzlev(j,:)
     275      latdeg=abs(latitude_deg(j))
     276
    275277c      CALL SW_venus_ve_1Dglobave(zrmu0, zfract,   ! pour moy globale
    276278c      CALL SW_venus_ve(zrmu0, zfract,
    … …  
    279281c     S        ztopsw,zsolsw,ZFSNET)
    280282
    281 c      CALL SW_venus_cl_1Dglobave(zrmu0, zfract,   ! pour moy globale
    282 c      CALL SW_venus_cl(zrmu0, zfract,
    283 c      CALL SW_venus_dc_1Dglobave(zrmu0, zfract,   ! pour moy globale
    284       CALL SW_venus_dc(zrmu0, zfract,
     283c      CALL SW_venus_cl_1Dglobave(zrmu0,zfract,   ! pour moy globale
     284c      CALL SW_venus_cl(zrmu0,zfract,
     285c      CALL SW_venus_dc_1Dglobave(zrmu0,zfract,   ! pour moy globale
     286c      CALL SW_venus_dc(zrmu0,zfract,
     287      CALL SW_venus_rh(zrmu0,zfract,latdeg,
     288c      CALL SW_venus_rh_1Dglobave(zrmu0,zfract,   ! pour moy globale
    285289     S        PPB,temp,
    286290     S        zheat,
Note: See TracChangeset for help on using the changeset viewer.