Changeset 2622 for trunk/LMDZ.VENUS/libf


Ignore:
Timestamp:
Feb 16, 2022, 4:26:05 PM (3 years ago)
Author:
slebonnois
Message:

SL: VENUS update (i) bug correction (2 bugs, phytrac and physiq), affected meam molec mass computations... (ii) updates for VCD 2.0 (iii) aeropacity: for latitudinal variations of the cloud distribution

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r2560 r2622  
    22         aerosol,reffrad,nueffrad,QREFvis3d,tau_col)
    33
    4        use radinc_h, only : L_TAUMAX,naerkind
     4       use radinc_h, only : naerkind
    55       use aerosol_mod
    6                  
     6                   
    77       implicit none
    88#include "YOMCST.h"
    … …  
    4949      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
    5050
    51       real aerosol0, obs_tau_col_aurora, pm
    52       real pcloud_deck, cloud_slope
    53 
    54       real dp_strato(ngrid)
    55       real dp_tropo(ngrid)
    56       real dp_layer(ngrid)
    57 
    58       INTEGER l,ig,iq,iaer,ia,ll
     51      INTEGER l,ig,iaer,ll
    5952
    6053      LOGICAL,SAVE :: firstcall=.true.
    6154!$OMP THREADPRIVATE(firstcall)
    62       REAL CBRT
    63       EXTERNAL CBRT
    64 
    65       ! for fixed dust profiles
    66       real topdust, expfactor, zp
    67       REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
    68       REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
    69 
    70       real CLFtot
     55      character*80 abort_message
     56      character*20 modname
    7157
    7258      !  for venus clouds
    73       real :: p_bot2,p_bot,p_top,p_top2,h_bot2,h_bot,h_top,h_top2
    74       real :: mode_dens,h_lay,nmode1,nmode2,nmode2p,nmode3,nmodeuv
    75 
     59      real,ALLOCATABLE,SAVE :: pbot2(:,:),pbot(:,:),ptop(:,:),ptop2(:,:)
     60      real,ALLOCATABLE,SAVE :: hbot2(:,:),hbot(:,:),htop(:,:),htop2(:,:)
     61      real,ALLOCATABLE,SAVE :: nmode(:,:)
     62!$OMP THREADPRIVATE(pbot2,pbot,ptop,ptop2)
     63!$OMP THREADPRIVATE(hbot2,hbot,htop,htop2)
     64!$OMP THREADPRIVATE(nmode)
     65      real :: mode_dens,hlay
     66
     67! ---------------------------------------------------------
     68! ---------------------------------------------------------
    7669! First call only
    7770      IF (firstcall) THEN
    78 
    79       ! identify tracers
    80         write(*,*) "Active aerosols found in aeropacity:"
    81 
    82         if (iaero_venus1.ne.0) then
     71     
     72        modname = 'aeropacity'
     73
     74! verify tracers
     75        write(*,*) "Active aerosols found in aeropacity, designed for Venus"
     76
     77        if (iaero_venus1.eq.1) then
    8378          print*,'iaero_venus1= ',iaero_venus1
    84         endif
    85         if (iaero_venus2.ne.0) then
     79        else
     80            abort_message='iaero_venus1 is not 1...'
     81            call abort_physic(modname,abort_message,1)
     82        endif
     83        if (iaero_venus2.eq.2) then
    8684          print*,'iaero_venus2= ',iaero_venus2
    87         endif
    88         if (iaero_venus2p.ne.0) then
     85        else
     86            abort_message='iaero_venus2 is not 2...'
     87            call abort_physic(modname,abort_message,1)
     88        endif
     89        if (iaero_venus2p.eq.3) then
    8990          print*,'iaero_venus2p= ',iaero_venus2p
    90         endif
    91         if (iaero_venus3.ne.0) then
     91        else
     92            abort_message='iaero_venus2p is not 3...'
     93            call abort_physic(modname,abort_message,1)
     94        endif
     95        if (iaero_venus3.eq.4) then
    9296          print*,'iaero_venus3= ',iaero_venus3
    93         endif
    94         if (iaero_venusUV.ne.0) then
     97        else
     98            abort_message='iaero_venus3 is not 4...'
     99            call abort_physic(modname,abort_message,1)
     100        endif
     101        if (iaero_venusUV.eq.5) then
    95102          print*,'iaero_venusUV= ',iaero_venusUV
    96         endif
     103        else
     104            abort_message='iaero_venus1UV is not 5...'
     105            call abort_physic(modname,abort_message,1)
     106        endif
     107
     108! cloud model
     109        allocate(pbot2(ngrid,naerkind),pbot(ngrid,naerkind))
     110        allocate(ptop(ngrid,naerkind),ptop2(ngrid,naerkind))
     111        allocate(hbot2(ngrid,naerkind),hbot(ngrid,naerkind))
     112        allocate(htop(ngrid,naerkind),htop2(ngrid,naerkind))
     113        allocate(nmode(ngrid,naerkind))
     114        call cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode)
    97115
    98116        firstcall=.false.
    99117      ENDIF ! of IF (firstcall)
    100 
    101 
    102 !     ---------------------------------------------------------
    103 
    104 !==================================================================
    105 !         Venus clouds (4 modes)
    106 !   S. Lebonnois (jan 2016)
    107 !==================================================================
    108 ! distributions from Haus et al, 2016
    109 ! mode             1      2      2p     3
    110 ! r (microns)     0.30   1.05   1.40   3.65
    111 ! sigma           1.56   1.29   1.23   1.28
    112 ! reff (microns)  0.49   1.23   1.56   4.25
    113 ! nueff           0.21   0.067  0.044  0.063
    114 ! (nueff=exp(ln^2 sigma)-1)
    115 !
    116 ! p_bot <=> zb ; p_top <=> zb+zc ; h_bot <=> Hlo ; h_top <=> Hup
    117 ! p<p_top: N(i+1)=N(i)*(p(i+1)/p(i))**(h_lay(i)/h_top)     
    118 !    h_lay=R(T(i)+T(i+1))/2/g  (in m)
    119 !    R=8.31/mu (mu in kg/mol)
    120 ! p>p_bot: N(i-1)=N(i)*(p(i)/p(i-1))**(h_lay(i)/h_bot)     
    121 ! N is in m-3
    122 !
    123 ! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*h_lay*(-dp)/p
    124 
    125 ! Mode 1
    126       if (iaero_venus1 .ne.0) then
    127           iaer=iaero_venus1
    128 
    129 !       1. Initialization
    130           aerosol(1:ngrid,1:nlayer,iaer)=0.0
    131 ! two scales heights for below and above max density (cf Haus et al 16)
    132           p_bot2 = 2.0e5 ! Pa   2.0e5
    133             h_bot2 = 5.0e3 ! m
    134           p_bot  = 1.2e5   !    1.2e5
    135             h_bot  = 1.0e3
    136           nmode1 = 1.935e8 ! m-3
    137           p_top  = 1.0e4
    138             h_top  = 3.5e3
    139           p_top2 = 4.5e2
    140             h_top2 = 2.0e3
    141 
    142 !       2. Opacity calculation
    143 
    144           DO ig=1,ngrid
     118! ---------------------------------------------------------
     119! ---------------------------------------------------------
     120
     121    aerosol(:,:,:) = 0.
     122        tau_col(:) = 0.
     123
     124      do ig=1,ngrid
     125        do iaer=1,naerkind
     126
     127!       Opacity calculation
     128
    145129! determine the approximate middle of the mode layer
    146130           ll=1
    147131           DO l=1,nlayer-1    ! high p to low p
    148              IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
     132             IF (pplay(ig,l) .gt. (ptop(ig,iaer)+pbot(ig,iaer))/2) ll=l
    149133           ENDDO
    150 ! from there go down and up for profile N
    151            mode_dens = nmode1  ! m-3
     134! from there go DOWN for profile N
     135           mode_dens = nmode(ig,iaer)  ! m-3
    152136           DO l=ll,1,-1
    153              h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
    154              IF (pplay(ig,l) .le. p_bot) THEN
    155                mode_dens = nmode1  ! m-3
    156              ELSEIF (pplay(ig,l) .gt. p_bot .and. pplay(ig,l) .le. p_bot2) THEN
    157                mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
    158              ELSEIF (pplay(ig,l) .gt. p_bot2) THEN
    159                mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot2)
     137             hlay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
     138             IF (pplay(ig,l) .le. pbot(ig,iaer)) THEN
     139               mode_dens = nmode(ig,iaer)  ! m-3
     140             ELSEIF (pplay(ig,l) .gt. pbot(ig,iaer) .and. pplay(ig,l) .le. pbot2(ig,iaer)) THEN
     141               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot(ig,iaer))
     142             ELSEIF (pplay(ig,l) .gt. pbot2(ig,iaer)) THEN
     143               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot2(ig,iaer))
    160144             ENDIF
    161145             mode_dens=max(1.e-30,mode_dens)
    162              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    163               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    164               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    165 !             aerosol(ig,l,iaer) = mode_dens
     146             if (iaer.lt.5) then
     147              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
     148               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
     149               mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
     150             else  ! for UV absorber
     151! normalized to 0.35 microns (peak of absorption)
     152              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
     153             endif
    166154           ENDDO
    167            mode_dens = nmode1  ! m-3
     155! ... then UP for profile N
     156           mode_dens = nmode(ig,iaer)  ! m-3
    168157           DO l=ll+1,nlayer-1
    169              h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
    170              IF (pplay(ig,l) .gt. p_top) THEN
    171                mode_dens = nmode1  ! m-3
    172              ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. p_top2) THEN
    173                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
    174              ELSEIF (pplay(ig,l) .le. p_top2) THEN
    175                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top2)
     158             hlay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
     159             IF (pplay(ig,l) .gt. ptop(ig,iaer)) THEN
     160               mode_dens = nmode(ig,iaer)  ! m-3
     161             ELSEIF (pplay(ig,l) .le. ptop(ig,iaer) .and. pplay(ig,l) .gt. ptop2(ig,iaer)) THEN
     162               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop(ig,iaer))
     163             ELSEIF (pplay(ig,l) .le. ptop2(ig,iaer)) THEN
     164               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop2(ig,iaer))
    176165             ENDIF
    177166             mode_dens=max(1.e-30,mode_dens)
    178              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    179               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    180               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    181 !             aerosol(ig,l,iaer) = mode_dens
     167             if (iaer.lt.5) then
     168              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
     169               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
     170               mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
     171             else  ! for UV absorber
     172! normalized to 0.35 microns (peak of absorption)
     173              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
     174             endif
    182175           ENDDO
    183           ENDDO
    184 
    185       end if ! mode 1
    186 
    187 ! Mode 2
    188       if (iaero_venus2 .ne.0) then
    189           iaer=iaero_venus2
    190 
    191 !       1. Initialization
    192           aerosol(1:ngrid,1:nlayer,iaer)=0.0
    193           p_bot = 1.0e4
    194             h_bot = 3.0e3
    195           nmode2 = 1.00e8 ! m-3
    196           p_top = 8.0e3
    197             h_top = 3.5e3
    198           p_top2 = 4.5e2
    199             h_top2 = 2.0e3
    200 
    201 !       2. Opacity calculation
    202 
    203           DO ig=1,ngrid
    204 ! determine the approximate middle of the mode layer
    205            ll=1
    206            DO l=1,nlayer-1    ! high p to low p
    207              IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
    208            ENDDO
    209 ! from there go down and up for profile N
    210            mode_dens = nmode2  ! m-3
    211            DO l=ll,1,-1
    212              h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
    213              IF (pplay(ig,l) .le. p_bot) THEN
    214                mode_dens = nmode2  ! m-3
    215              ELSEIF (pplay(ig,l) .gt. p_bot) THEN
    216                mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
    217              ENDIF
    218              mode_dens=max(1.e-30,mode_dens)
    219              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    220               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    221               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    222 !             aerosol(ig,l,iaer) = mode_dens
    223            ENDDO
    224            mode_dens = nmode2  ! m-3
    225            DO l=ll+1,nlayer-1
    226              h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
    227              IF (pplay(ig,l) .gt. p_top) THEN
    228                mode_dens = nmode2  ! m-3
    229              ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. p_top2) THEN
    230                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
    231              ELSEIF (pplay(ig,l) .le. p_top2) THEN
    232                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top2)
    233              ENDIF
    234              mode_dens=max(1.e-30,mode_dens)
    235              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    236               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    237               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    238 !             aerosol(ig,l,iaer) = mode_dens
    239            ENDDO
    240           ENDDO
    241 
    242       end if ! mode 2
    243 
    244 ! Mode 2p
    245       if (iaero_venus2p .ne.0) then
    246           iaer=iaero_venus2p
    247 
    248 !       1. Initialization
    249           aerosol(1:ngrid,1:nlayer,iaer)=0.0
    250           p_bot = 1.2e5  ! 1.2e5
    251             h_bot = 0.1e3
    252           nmode2p = 5.0e7 ! m-3
    253           p_top = 2.3e4
    254             h_top = 1.0e3
    255 
    256 !       2. Opacity calculation
    257 
    258           DO ig=1,ngrid
    259 ! determine the approximate middle of the mode layer
    260            ll=1
    261            DO l=1,nlayer-1    ! high p to low p
    262              IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
    263            ENDDO
    264 ! from there go down and up for profile N
    265            mode_dens = nmode2p  ! m-3
    266            DO l=ll,1,-1
    267              h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
    268              IF (pplay(ig,l) .le. p_bot) THEN
    269                mode_dens = nmode2p  ! m-3
    270              ELSEIF (pplay(ig,l) .gt. p_bot) THEN
    271                mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
    272              ENDIF
    273              mode_dens=max(1.e-30,mode_dens)
    274              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    275               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    276               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    277 !             aerosol(ig,l,iaer) = mode_dens
    278            ENDDO
    279            mode_dens = nmode2p  ! m-3
    280            DO l=ll+1,nlayer-1
    281              h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
    282              IF (pplay(ig,l) .gt. p_top) THEN
    283                mode_dens = nmode2p  ! m-3
    284              ELSEIF (pplay(ig,l) .le. p_top) THEN
    285                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
    286              ENDIF
    287              mode_dens=max(1.e-30,mode_dens)
    288              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    289               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    290               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    291 !             aerosol(ig,l,iaer) = mode_dens
    292            ENDDO
    293           ENDDO
    294 
    295       end if ! mode 2p
    296 
    297 ! Mode 3
    298       if (iaero_venus3 .ne.0) then
    299           iaer=iaero_venus3
    300 
    301 !       1. Initialization
    302           aerosol(1:ngrid,1:nlayer,iaer)=0.0
    303           p_bot = 1.2e5  ! 1.2e5
    304             h_bot = 0.5e3
    305           nmode3 = 1.4e7 ! m-3
    306           p_top = 3.9e4
    307             h_top = 1.0e3
    308 
    309 !       2. Opacity calculation
    310 
    311           DO ig=1,ngrid
    312 ! determine the approximate middle of the mode layer
    313            ll=1
    314            DO l=1,nlayer-1    ! high p to low p
    315              IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
    316            ENDDO
    317 ! from there go down and up for profile N
    318            mode_dens = nmode3  ! m-3
    319            DO l=ll,1,-1
    320              h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
    321              IF (pplay(ig,l) .le. p_bot) THEN
    322                mode_dens = nmode3  ! m-3
    323              ELSEIF (pplay(ig,l) .gt. p_bot) THEN
    324                mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
    325              ENDIF
    326              mode_dens=max(1.e-30,mode_dens)
    327              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    328               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    329               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    330 !             aerosol(ig,l,iaer) = mode_dens
    331            ENDDO
    332            mode_dens = nmode3  ! m-3
    333            DO l=ll+1,nlayer-1
    334              h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
    335              IF (pplay(ig,l) .gt. p_top) THEN
    336                mode_dens = nmode3  ! m-3
    337              ELSEIF (pplay(ig,l) .le. p_top) THEN
    338                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
    339              ENDIF
    340              mode_dens=max(1.e-30,mode_dens)
    341              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
    342               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
    343               mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
    344 !             aerosol(ig,l,iaer) = mode_dens
    345            ENDDO
    346           ENDDO
    347 
    348       end if ! mode 3
    349 
    350 ! UV absorber
    351       if (iaero_venusUV .ne.0) then
    352           iaer=iaero_venusUV
    353 
    354 !       1. Initialization
    355           aerosol(1:ngrid,1:nlayer,iaer)=0.0
    356           p_bot = 3.3e4  ! 58 km
    357             h_bot = 1.0e3
    358           nmodeuv = 1.00e7 !m-3
    359           p_top = 3.7e3 ! 70 km
    360             h_top = 1.0e3
    361 
    362 !       2. Opacity calculation
    363 
    364           DO ig=1,ngrid
    365 ! determine the approximate middle of the mode layer
    366            ll=1
    367            DO l=1,nlayer-1    ! high p to low p
    368              IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
    369            ENDDO
    370 ! from there go down and up for profile N
    371            mode_dens = nmodeuv ! m-3
    372            DO l=ll,1,-1
    373              h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
    374              IF (pplay(ig,l) .le. p_bot) THEN
    375                mode_dens = nmodeuv ! m-3
    376              ELSEIF (pplay(ig,l) .gt. p_bot) THEN
    377                mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
    378              ENDIF
    379              mode_dens=max(1.e-30,mode_dens)
    380 ! normalized to 0.35 microns (peak of absorption)
    381              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
    382            ENDDO
    383            mode_dens = nmodeuv ! m-3
    384            DO l=ll+1,nlayer-1
    385              h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
    386              IF (pplay(ig,l) .gt. p_top) THEN
    387                mode_dens = nmodeuv ! m-3
    388              ELSEIF (pplay(ig,l) .le. p_top) THEN
    389                mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
    390              ENDIF
    391              mode_dens=max(1.e-30,mode_dens)
    392 ! normalized to 0.35 microns (peak of absorption)
    393              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
    394            ENDDO
    395           ENDDO
    396 
    397 !       3. Re-normalize to Haus et al 2015 total column optical depth
    398          tau_col(:)=0.0
    399          DO l=1,nlayer
    400           DO ig=1,ngrid
     176
     177           if (iaer.eq.5) then  ! for UV absorber
     178            DO l=1,nlayer
    401179               tau_col(ig) = tau_col(ig) &
    402180                     + aerosol(ig,l,iaer)
    403181            ENDDO
    404          ENDDO
    405          DO ig=1,ngrid
    406            DO l=1,nlayer-1
     182            DO l=1,nlayer
    407183                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig)
    408            ENDDO
    409          ENDDO
    410 
    411       end if ! UV absorber
     184            ENDDO
     185           endif
     186
     187        enddo
     188      enddo
    412189
    413190!==================================================================
    … …  
    429206
    430207      tau_col(:)=0.0
    431       do iaer = 1, naerkind
     208      do ig=1,ngrid
    432209         do l=1,nlayer
    433             do ig=1,ngrid
     210            do iaer = 1, naerkind
    434211               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
    435212            end do
    … …  
    459236         endif
    460237      end do
     238
    461239      return
    462240    end subroutine aeropacity
     241
     242! --------------------------------------------------------------------------
     243! --------------------------------------------------------------------------
     244subroutine cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode)
     245
     246!==================================================================
     247!         Venus clouds (4 modes)
     248!   S. Lebonnois (jan 2016)
     249!==================================================================
     250! distributions from Haus et al, 2016
     251! ------------------------------------
     252!
     253! mode             1      2      2p     3
     254! r (microns)     0.30   1.05   1.40   3.65
     255! sigma           1.56   1.29   1.23   1.28
     256! reff (microns)  0.49   1.23   1.56   4.25
     257! nueff           0.21   0.067  0.044  0.063
     258! (nueff=exp(ln^2 sigma)-1)
     259!
     260! pbot <=> zb ; ptop <=> zb+zc ; hbot <=> Hlo ; htop <=> Hup
     261! p<ptop: N(i+1)=N(i)*(p(i+1)/p(i))**(hlay(i)/htop)     
     262!    hlay=R(T(i)+T(i+1))/2/g  (in m)
     263!    R=8.31/mu (mu in kg/mol)
     264! p>pbot: N(i-1)=N(i)*(p(i)/p(i-1))**(hlay(i)/hbot)     
     265! N is in m-3
     266!
     267! values in each mode below from Table 1 and text (p.186)
     268!
     269! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*hlay*(-dp)/p
     270!
     271! Latitudinal dependence (values from tables 3 and 4):
     272! mode factor MF12(lat) for modes 1 and 2
     273! mode 2' unchanged
     274! mode factor MF3(lat) for mode 3
     275! + for mode 2 only : pbot(lat), ptop(lat), htop(lat), htop2(lat)
     276!==================================================================
     277
     278      use radinc_h, only : naerkind
     279      USE geometry_mod, ONLY: latitude_deg             
     280      implicit none
     281       
     282      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
     283      real,INTENT(out) :: pbot2(ngrid,naerkind),pbot(ngrid,naerkind)
     284      real,INTENT(out) :: ptop(ngrid,naerkind),ptop2(ngrid,naerkind)
     285      real,INTENT(out) :: hbot2(ngrid,naerkind),hbot(ngrid,naerkind)
     286      real,INTENT(out) :: htop(ngrid,naerkind),htop2(ngrid,naerkind)
     287      real,INTENT(out) :: nmode(ngrid,naerkind)
    463288     
     289! For latitudinal dependence
     290      integer :: i,j,ilat
     291      real :: latit,factlat
     292      integer,parameter :: nlat_H16 = 16
     293      real :: lat_H16(nlat_H16)
     294      real :: MF12(nlat_H16),MF3(nlat_H16)
     295      real :: pbotM2(nlat_H16),ptopM2(nlat_H16),htopM2(nlat_H16)
     296
     297!-------------------------
     298!!   Equatorial model 
     299!-------------------------
     300! initialization
     301      pbot2(:,:) = 1.e8
     302      pbot(:,:)  = 1.e8
     303      ptop2(:,:) = 0.
     304      ptop(:,:)  = 0.
     305      hbot2(:,:) = 1.
     306      hbot(:,:)  = 1.
     307      htop2(:,:) = 1.
     308      htop(:,:)  = 1.
     309      nmode(:,:) = 0.
     310       
     311! table of values(lat,mode)
     312
     313! Mode 1
     314      pbot2(:,1) = 2.0e5 ! Pa
     315      hbot2(:,1) = 5.0e3 ! m
     316      pbot(:,1)  = 1.2e5
     317      hbot(:,1)  = 1.0e3
     318      nmode(:,1) = 1.935e8 ! m-3
     319      ptop(:,1)  = 1.0e4
     320      htop(:,1)  = 3.5e3
     321      ptop2(:,1) = 4.5e2
     322      htop2(:,1) = 2.0e3
     323
     324! Mode 2
     325      pbot(:,2)  = 1.0e4
     326      hbot(:,2)  = 3.0e3
     327      nmode(:,2) = 1.00e8 ! m-3
     328      ptop(:,2)  = 8.0e3
     329      htop(:,2)  = 3.5e3
     330      ptop2(:,2) = 4.5e2
     331      htop2(:,2) = 2.0e3
     332
     333! Mode 2p
     334      pbot(:,3)  = 1.2e5
     335      hbot(:,3)  = 0.1e3
     336      nmode(:,3) = 5.0e7 ! m-3
     337      ptop(:,3)  = 2.3e4
     338      htop(:,3)  = 1.0e3
     339
     340! Mode 3
     341      pbot(:,4)  = 1.2e5
     342      hbot(:,4)  = 0.5e3
     343      nmode(:,4) = 1.4e7 ! m-3
     344      ptop(:,4)  = 3.9e4
     345      htop(:,4)  = 1.0e3
     346
     347! UV absorber
     348      pbot(:,5)  = 3.3e4  ! 58 km
     349      hbot(:,5)  = 1.0e3
     350      nmode(:,5) = 1.00e7 !m-3
     351      ptop(:,5)  = 3.7e3  ! 70 km
     352      htop(:,5)  = 1.0e3
     353
     354!----------------------------
     355!!   Latitudinal variations 
     356!----------------------------
     357      lat_H16(:) = (/0.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90./)
     358      MF12(:) = (/0.98,0.98,0.99,1.00,0.98,0.94,0.86,0.81,0.73,0.67,0.64,0.61,0.59,0.47,0.36,0.36/)
     359      MF3(:)  = (/1.30,1.30,1.26,1.23,1.17,1.13,1.06,1.03,1.04,1.09,1.22,1.51,1.82,2.02,2.09,2.09/)
     360      pbotM2(:) = (/1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,9.4e3, &
     361                    9.2e3,9.1e3,9.6e3,1.14e4,1.21e4,1.25e4,1.25e4/)
     362      ptopM2(:) = (/8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,7.7e3, &
     363                    7.5e3,7.5e3,8.0e3,9.2e3,1.0e4,1.06e4,1.06e4/)
     364      htopM2(:) = (/3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.4e3, &
     365                    3.2e3,2.6e3,2.0e3,2.0e3,0.6e3,0.5e3,0.5e3/)
     366     
     367      do i=1,ngrid
     368        latit = abs(latitude_deg(i))
     369        ilat=2
     370        do j=2,nlat_H16
     371           if (latit.gt.lat_H16(j-1)) ilat=j
     372        enddo
     373        factlat    = (lat_H16(ilat)-latit)/(lat_H16(ilat)-lat_H16(ilat-1))
     374        pbot(i,2)  = pbotM2(ilat)*(1.-factlat) + pbotM2(ilat-1)*factlat
     375        ptop(i,2)  = ptopM2(ilat)*(1.-factlat) + ptopM2(ilat-1)*factlat
     376        htop(i,2)  = htopM2(ilat)*(1.-factlat) + htopM2(ilat-1)*factlat
     377        htop2(i,2) = min(htop2(i,2),htop(i,2))
     378        nmode(i,1) = nmode(i,1)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat)
     379        nmode(i,2) = nmode(i,2)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat)
     380        nmode(i,3) = nmode(i,3)*(MF3(ilat) *(1.-factlat)+MF3(ilat-1) *factlat)
     381      enddo
     382     
     383      return
     384end subroutine cloud_haus16

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

    r2198 r2622  
    1111      parameter (Pdiff=1.e-1)      ! pressure below which diffusion is computed
    1212      parameter (tdiffmin=5.d0)
    13       parameter (dzres=0.1d0)    ! grid resolution (km) for diffusion
     13      parameter (dzres=0.2d0)    ! grid resolution (km) for diffusion
    1414     

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

    r2464 r2622  
    7676         xabsi(3,i)  = rm(i,i_o)
    7777         xabsi(8,i)  = rm(i,i_n2)
    78          xabsi(11,i)  = rm(i,i_co)
     78         xabsi(11,i) = rm(i,i_co)
    7979
    8080c         xabsi(6,i)  = rm(i,i_h2o2)

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

    r2580 r2622  
    9696! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid
    9797
    98 !     REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity  TestGW6
     98!----------------------------------------
     99! VCD 1.1 tuning
     100!     REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity
    99101!----------------------------------------
    100102! VCD 2.0 tuning
    101      REAL, parameter:: cmax = 301.   ! Max horizontal absolute phase velocity  TestGW6
    102 !----------------------------------------
    103 
     103      REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity
     104!----------------------------------------
     105   
    104106! generic : cmax=zonal wind at launch
    105107     REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity
    … …  
    126128
    127129     ! 0.3.2 Parameters of waves dissipations
    128      REAL, parameter:: sat   = 0.85   ! saturation parameter
    129      REAL, parameter:: rdiss = 0.1    ! coefficient of dissipation
     130!----------------------------------------
     131! VCD 1.1 tuning
     132!    REAL, parameter:: sat   = 0.85   ! saturation parameter
     133!    REAL, parameter:: rdiss = 0.1    ! coefficient of dissipation
     134!----------------------------------------
     135! VCD 2.0 tuning
     136     REAL, parameter:: sat   = 0.6    ! saturation parameter
     137     REAL, parameter:: rdiss = 8.e-4  ! coefficient of dissipation
     138!----------------------------------------
    130139     REAL, parameter:: zoisec = 1.e-8 ! security for intrinsic freguency
    131140
    132141     ! 0.3.3 Background flow at 1/2 levels and vertical coordinate
    133142     REAL H0bis(ngrid, nlayer)       ! characteristic height of the atmosphere
     143     REAL min_k(ngrid)               ! min(kstar,kmin)
    134144     REAL, save::  H0                ! characteristic height of the atmosphere
    135      REAL, parameter:: pr = 5e5      ! Reference pressure [Pa]   ! VENUS: cloud layer
     145     REAL, save::  MDR               ! characteristic mass density 
     146!----------------------------------------
     147! VCD 1.1 tuning
     148!    REAL, parameter:: pr = 5e5      ! Reference pressure [Pa]   ! VENUS: cloud layer
     149!----------------------------------------
     150! VCD 2.0 tuning
     151     REAL, parameter:: pr = 5e4      ! Reference pressure [Pa]   ! VENUS: cloud layer
     152!----------------------------------------
    136153     REAL, parameter:: tr = 300.     ! Reference temperature [K] ! VENUS: cloud layer
    137154     REAL zh(ngrid, nlayer + 1)      ! Log-pressure altitude (constant H0)
    … …  
    156173        ! Characteristic vertical scale height
    157174        H0 = RD * tr / RG
     175        ! Characteristic mass density at launch altitude
     176        MDR = pr / ( RD * tr )
    158177        ! Control
    159178        if (deltat .LT. dtime) THEN
    … …  
    214233     uh(:, nlayer + 1) = uu(:, nlayer)
    215234     vh(:, nlayer + 1) = vv(:, nlayer)
     235
     236     ! TN+GG April/June 2020
     237     ! "Individual waves are not supposed
     238     ! to occupy the entire domain, but
     239     ! only a faction of it" Lott+2012
     240     ! minimum value of k between kmin and kstar
     241     DO ii = 1, ngrid
     242        kstar = RPI / sqrt(cell_area(ii))
     243        min_k(ii) = max(kmin,kstar)
     244     end DO
    216245
    217246! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
    … …  
    235264           ! horizontal wavenumber amplitude
    236265!           zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
    237            ! TN+GG April/June 2020
    238            ! "Individual waves are not supposed
    239            ! to occupy the entire domain, but
    240            ! only a faction of it" Lott+2012
    241            kstar = RPI / sqrt(cell_area(ii))
    242            zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2
     266           zk(jw, ii) = min_k(ii) + (kmax - min_k(ii)) * ran_num_2
    243267
    244268           ! horizontal phase speed
    … …  
    320344           ! Intrinsic frequency
    321345           intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1)    &
    322                       - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
     346                                          - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
    323347           ! Vertical velocity in flux formulation
    324348           wwp(jw, :) = min(                                                              &
    … …  
    336360                      * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1)                        &
    337361                      * exp(-zh(:, ll + 1) / H0)                                          &
    338                       * sat**2 * kmin**2 / zk(jw, :)**4)
     362!! Correction for VCD 2.0
     363!                      * sat**2 * kmin**2 / zk(jw, :)**4)
     364                      * sat**2 * min_k(:)**2 * MDR / zk(jw, :)**4)                       
    339365        end DO
    340366

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

    r2580 r2622  
    15691569!-- TuneC
    15701570! VCD 1.1 tuning
    1571     v_phot(65:78,4) = v_phot(65:78,4)*10.
    1572     v_phot(52:59,3) = v_phot(52:59,3)*5.
     1571!   v_phot(65:78,4) = v_phot(65:78,4)*10.
     1572!   v_phot(52:59,3) = v_phot(52:59,3)*5.
    15731573!--
    15741574!-- TuneE
    15751575! VCD 2.0 tuning
    1576 !   v_phot(65:90,4) = v_phot(65:90,4)*10.
     1576    v_phot(65:90,4) = v_phot(65:90,4)*10.
    15771577!--
    15781578!   v_phot(:,4) = v_phot(:,4)*10.

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

    r2580 r2622  
    14491449            mmean(:,:) = 0.
    14501450            do iq = 1,nqmax - nmicro
    1451                mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)*m_tr(iq)
    1452                rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     1451               mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
    14531452            enddo
     1453            mmean(:,:) = 1./mmean(:,:)
     1454            rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
    14541455         endif
    14551456
    … …  
    18121813! Obsolete
    18131814! but used for VCD 1.1
    1814       call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
    1815      e               t_seri, u_seri, v_seri, paprs(klon/2+1,:),
     1815!     call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
     1816!    e               t_seri, u_seri, v_seri, paprs(klon/2+1,:),
     1817!    o               zustrhi,zvstrhi,
     1818!    o               d_t_hin, d_u_hin, d_v_hin)
     1819
     1820! New routine based on Generic
     1821! used after VCD 1.1, for VCD 2.0
     1822      call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,
     1823     e               t_seri, u_seri, v_seri,
    18161824     o               zustrhi,zvstrhi,
    18171825     o               d_t_hin, d_u_hin, d_v_hin)
    1818 
    1819 ! New routine based on Generic
    1820 ! used after VCD 1.1
    1821 !     call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,
    1822 !    e               t_seri, u_seri, v_seri,
    1823 !    o               zustrhi,zvstrhi,
    1824 !    o               d_t_hin, d_u_hin, d_v_hin)
    18251826
    18261827c  ajout des tendances

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

    r2580 r2622  
    7474      if (debutphy) then
    7575     
    76 !!--- Adjustment of Helium amount
    77 !        if (i_he/=0) then
    78 !           trac(:,:,i_he)=trac(:,:,i_he)*20.
    79 !        endif
    80 !!---
     76!--- Adjustment of Helium amount
     77!       if (i_he/=0) then
     78!          trac(:,:,i_he)=trac(:,:,i_he)/2.
     79!       endif
     80!---
    8181
    8282!-------------------------------------------------------------------
    … …  
    8686
    8787!!! in this reinitialisation, trac is VOLUME mixing ratio
    88             print*, "Tracers are re-initialised"
    89             trac(:,:,:) = 1.0e-30
    90 
    91             if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
    92      $           .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0)
    93      $           .and. (i_co2 /= 0)) then
    94 
    95                trac(:,1:22,i_ocs) = 3.e-6
    96                trac(:,1:22,i_co)  = 25.e-6
    97                trac(:,:,i_hcl)    = 0.4e-6
    98                trac(:,1:22,i_so2) = 7.e-6
    99                trac(:,1:22,i_h2o) = 30.e-6
    100                trac(:,:,i_n2)     = 0.35e-1
     88! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     89!           convert mass to volume mixing ratio
     90            do iq = 1,nqmax - nmicro
     91               trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
     92            end do
     93            print*, "SO2 is re-initialised"
     94            if (i_so2 /= 0) then
     95               trac(:,1:22,i_so2) = 10.e-6
     96
     97! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     98!           print*, "Tracers are re-initialised"
     99!           trac(:,:,:) = 1.0e-30
     100
     101!           if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
     102!    $           .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0)
     103!    $           .and. (i_co2 /= 0)) then
     104
     105!              trac(:,1:22,i_ocs) = 3.e-6
     106!              trac(:,1:22,i_co)  = 25.e-6
     107!              trac(:,:,i_hcl)    = 0.4e-6
     108!              trac(:,1:22,i_so2) = 7.e-6
     109!              trac(:,1:22,i_h2o) = 30.e-6
     110!              trac(:,:,i_n2)     = 0.35e-1
     111!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    101112   
    102113!          ensure that sum of mixing ratios = 1
    … …  
    124135            do iq = 1,nqmax - nmicro
    125136               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
    126                rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
    127137            enddo
     138            rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
    128139
    129140!           convert volume to mass mixing ratio
    … …  
    293304            mmean(:,:) = 0.
    294305            do iq = 1,nqmax - nmicro
    295                mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
    296                rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     306               mmean(:,:) = mmean(:,:)+ztrac(:,:,iq)*m_tr(iq)
    297307            enddo
    298 
     308            rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     309
     310!===================================================================
     311!     convert volume to mass mixing ratio / then tendencies in mmr
     312!===================================================================
    299313!     gas phase
    300314
Note: See TracChangeset for help on using the changeset viewer.