Ignore:
Timestamp:
Jan 27, 2016, 10:42:32 AM (8 years ago)
Author:
idelkadi
Message:

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cosp/dsd.F90

    r1907 r2428  
    1   subroutine dsd(Q,Re,D,N,nsizes,dtype,rho_a,tc, &
    2              dmin,dmax,apm,bpm,rho_c,p1,p2,p3,fc,scaled)
     1  subroutine dsd(Q,Re,Np,D,N,nsizes,dtype,rho_a,tk, &
     2             dmin,dmax,apm,bpm,rho_c,p1,p2,p3)
    33  use array_lib
    44  use math_lib
    … …  
    77! Purpose:
    88!   Create a discrete drop size distribution
    9 !   Part of QuickBeam v1.03 by John Haynes
     9!
     10!   Starting with Quickbeam V3, this routine now allows input of
     11!   both effective radius (Re) and total number concentration (Nt)
     12!   Roj Marchand July 2010
     13!
     14!   The version in Quickbeam v.104 was modified to allow Re but not Nt
     15!   This is a significantly modified form for the version     
     16!
     17!   Originally Part of QuickBeam v1.03 by John Haynes
    1018!   http://reef.atmos.colostate.edu/haynes/radarsim
    1119!
    1220! Inputs:
     21!
    1322!   [Q]        hydrometeor mixing ratio (g/kg)
    14 !   [Re]       Optional Effective Radius (microns).  0 = use default.
    15 !   [D]        discrete drop sizes (um)
     23!   [Re]       Optional Effective Radius (microns).  0 = use defaults (p1, p2, p3)
     24!
     25!   [D]        array of discrete drop sizes (um) where we desire to know the number concentraiton n(D).
    1626!   [nsizes]   number of elements of [D]
     27!
    1728!   [dtype]    distribution type
    1829!   [rho_a]    ambient air density (kg m^-3)
    19 !   [tc]       temperature (C)
     30!   [tk]       temperature (K)
    2031!   [dmin]     minimum size cutoff (um)
    2132!   [dmax]     maximum size cutoff (um)
    … …  
    2435!
    2536! Input/Output:
    26 !   [fc]       scaling factor for the distribution
    27 !   [scaled]   has this hydrometeor type been scaled?
    2837!   [apm]      a parameter for mass (kg m^[-bpm])
    2938!   [bmp]      b params for mass
    … …  
    4150!   01/31/06  Port from IDL to Fortran 90
    4251!   07/07/06  Rewritten for variable DSD's
    43 !   10/02/06  Rewritten using scaling factors (Roger Marchand and JMH)
     52!   10/02/06  Rewritten using scaling factors (Roger Marchand and JMH), Re added V1.04
     53!   July 2020 "N Scale factors" (variable fc) removed (Roj Marchand).
    4454 
    4555! ----- INPUTS ----- 
    4656 
    47   integer*4, intent(in) :: nsizes
     57  integer, intent(in) :: nsizes
    4858  integer, intent(in) :: dtype
    49   real*8, intent(in) :: Q,D(nsizes),rho_a,tc,dmin,dmax, &
    50     rho_c,p1,p2,p3
    51    
    52 ! ----- INPUT/OUTPUT -----
    53 
    54   real*8, intent(inout) :: fc(nsizes),apm,bpm,Re
    55   logical, intent(inout) :: scaled 
    56    
     59  real*8, intent(in)  :: Q,Re,Np,D(nsizes)
     60  real*8, intent(in)  :: rho_a,tk,dmin,dmax,rho_c,p1,p2,p3
     61   
     62  real*8, intent(inout) :: apm,bpm 
     63 
    5764! ----- OUTPUTS -----
    5865
    … …  
    6067 
    6168! ----- INTERNAL -----
    62  
     69 
     70   real*8 :: fc(nsizes)
     71
    6372  real*8 :: &
    64   N0,D0,vu,np,dm,ld, &                  ! gamma, exponential variables
    65   dmin_mm,dmax_mm,ahp,bhp, &            ! power law variables
    66   rg,log_sigma_g, &                     ! lognormal variables
    67   rho_e                                 ! particle density (kg m^-3)
     73  N0,D0,vu,local_np,dm,ld, &            ! gamma, exponential variables
     74  dmin_mm,dmax_mm,ahp,bhp, &        ! power law variables
     75  rg,log_sigma_g, &         ! lognormal variables
     76  rho_e                 ! particle density (kg m^-3)
    6877 
    6978  real*8 :: tmp1, tmp2
    70   real*8 :: pi,rc
     79  real*8 :: pi,rc,tc
    7180
    7281  integer k,lidx,uidx
    7382
     83  tc = tk - 273.15
    7484  pi = acos(-1.0)
    7585 
    76 ! // if density is constant, store equivalent values for apm and bpm
     86  ! // if density is constant, store equivalent values for apm and bpm
    7787  if ((rho_c > 0) .and. (apm < 0)) then
    7888    apm = (pi/6)*rho_c
    … …  
    8090  endif
    8191 
     92  ! will preferentially use Re input over Np.
     93  ! if only Np given then calculate Re
     94  ! if neigher than use other defaults (p1,p2,p3) following quickbeam documentation
     95  if(Re==0 .and. Np>0) then
     96   
     97        call calc_Re(Q,Np,rho_a, &
     98             dtype,dmin,dmax,apm,bpm,rho_c,p1,p2,p3, &
     99             Re)
     100  endif
     101 
     102 
    82103  select case(dtype)
    83104 
    … …  
    85106! // modified gamma                                        !
    86107! ---------------------------------------------------------!
    87 ! :: N0 = total number concentration (m^-3)
    88 ! :: np = fixed number concentration (kg^-1)
     108! :: np = total number concentration
    89109! :: D0 = characteristic diameter (um)
    90 ! :: dm = mean diameter (um)
     110! :: dm = mean diameter (um) - first moment over zeroth moment
    91111! :: vu = distribution width parameter
    92112
    93113  case(1) 
    94     if (abs(p1+1) < 1E-8) then
    95 
    96 !     // D0, vu are given 
     114 
     115    if( abs(p3+2) < 1E-8) then
     116 
     117    if( Np>1E-30) then
     118   
     119        ! Morrison scheme with Martin 1994 shape parameter (NOTE: vu = pc +1)
     120        ! fixed Roj. Dec. 2010 -- after comment by S. Mcfarlane
     121        vu = (1/(0.2714 + 0.00057145*Np*rho_a*1E-6))**2.0 ! units of Nt = Np*rhoa = #/cm^3
     122    else
     123        print *, 'Error: Must specify a value for Np in each volume', &
     124             ' with Morrison/Martin Scheme.'
     125            stop   
     126    endif
     127   
     128    elseif (abs(p3+1) > 1E-8) then
     129
     130      ! vu is fixed in hp structure 
    97131      vu = p3
    98      
    99       if(Re.le.0) then
    100         dm = p2
    101         D0 = gamma(vu)/gamma(vu+1)*dm
    102       else
    103         D0 = 2.0*Re*gamma(vu+2)/gamma(vu+3)
    104       endif
    105      
    106       if (scaled .eqv. .false.) then
    107      
     132
     133    else
     134
     135      ! vu isn't specified
     136     
     137      print *, 'Error: Must specify a value for vu for Modified Gamma distribution'
     138      stop   
     139     
     140    endif
     141     
     142      if(Re>0) then
     143
     144    D0 = 2.0*Re*gamma(vu+2)/gamma(vu+3)
     145   
    108146        fc = ( &
    109147             ((D*1E-6)**(vu-1)*exp(-1*D/D0)) / &
    110148             (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm)) &
    111              ) * 1E-12
    112         scaled = .true.
    113 
    114       endif       
    115 
    116       N = fc*rho_a*(Q*1E-3)
    117    
    118     elseif (abs(p2+1) < 1E-8) then
    119 
    120 !     // N0, vu are given   
    121       np = p1
    122       vu = p3
    123       tmp1 = (Q*1E-3)**(1./bpm)
    124      
    125       if (scaled .eqv. .false.) then
    126 
    127         fc = (D*1E-6 / (gamma(vu)/(apm*np*gamma(vu+bpm)))** &
     149         ) * 1E-12
     150
     151        N = fc*rho_a*(Q*1E-3)
     152       
     153      elseif( p2+1 > 1E-8) then     ! use default value for MEAN diameter
     154     
     155        dm = p2
     156    D0 = gamma(vu)/gamma(vu+1)*dm
     157
     158        fc = ( &
     159             ((D*1E-6)**(vu-1)*exp(-1*D/D0)) / &
     160             (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm)) &
     161         ) * 1E-12
     162
     163        N = fc*rho_a*(Q*1E-3)
     164       
     165      elseif(abs(p3+1) > 1E-8)  then! use default number concentration
     166   
     167        local_np = p1 ! total number concentration / pa check
     168   
     169    tmp1 = (Q*1E-3)**(1./bpm)
     170     
     171        fc = (D*1E-6 / (gamma(vu)/(apm*local_np*gamma(vu+bpm)))** &
    128172             (1./bpm))**vu
    129              
    130         scaled = .true.
    131 
     173         
     174        N = ( &
     175          (rho_a*local_np*fc*(D*1E-6)**(-1.))/(gamma(vu)*tmp1**vu) * &
     176          exp(-1.*fc**(1./vu)/tmp1) &
     177      ) * 1E-12
     178
     179      else
     180     
     181        print *, 'Error:  No default value for Dm or Np provided!  '
     182        stop
     183       
    132184      endif
    133 
    134       N = ( &
    135           (rho_a*np*fc*(D*1E-6)**(-1.))/(gamma(vu)*tmp1**vu) * &
    136           exp(-1.*fc**(1./vu)/tmp1) &
    137           ) * 1E-12
    138 
    139     else
    140 
    141 !     // vu isn't given
    142       print *, 'Error: Must specify a value for vu'
    143       stop
    144    
    145     endif
     185   
    146186   
    147187! ---------------------------------------------------------!
    … …  
    152192
    153193  case(2)
    154     if (abs(p1+1) > 1E-8) then
    155 
    156 !     // N0 has been specified, determine ld
    157       N0 = p1
    158 
    159       if(Re>0) then
    160 
    161         ! if Re is set and No is set than the distribution is fully defined.
    162         ! so we assume Re and No have already been chosen consistant with 
    163         ! the water content, Q.
    164 
    165         ! print *,'using Re pass ...'
    166 
    167         ld = 1.5/Re   ! units 1/um
    168 
    169         N = ( &
    170                 N0*exp(-1*ld*D) &
    171         ) * 1E-12
    172    
    173       else
    174 
    175         tmp1 = 1./(1.+bpm)
    176      
    177         if (scaled .eqv. .false.) then
    178                 fc = ((apm*gamma(1.+bpm)*N0)**tmp1)*(D*1E-6)
    179                 scaled = .true.
    180 
    181         endif
     194 
     195    if(Re>0) then
     196 
     197        ld = 1.5/Re   ! units 1/um
     198       
     199    fc = (ld*1E6)**(1.+bpm)/(apm*gamma(1+bpm))* &
     200               exp(-1.*(ld*1E6)*(D*1E-6))*1E-12
    182201     
    183         N = ( &
    184                 N0*exp(-1.*fc*(1./(rho_a*Q*1E-3))**tmp1) &
    185         ) * 1E-12
    186 
    187       endif     
     202        N = fc*rho_a*(Q*1E-3)
     203       
     204    elseif (abs(p1+1) > 1E-8) then
     205
     206    ! use N0 default value
     207   
     208        N0 = p1
     209
     210        tmp1 = 1./(1.+bpm)
     211 
     212        fc = ((apm*gamma(1.+bpm)*N0)**tmp1)*(D*1E-6)
     213 
     214        N = ( &
     215            N0*exp(-1.*fc*(1./(rho_a*Q*1E-3))**tmp1) &
     216    ) * 1E-12
    188217
    189218    elseif (abs(p2+1) > 1E-8) then
    190219
    191 !     // ld has been specified, determine N0
    192       ld = p2
    193 
    194       if (scaled .eqv. .false.) then
     220    !     used default value for lambda
     221        ld = p2
    195222
    196223        fc = (ld*1E6)**(1.+bpm)/(apm*gamma(1+bpm))* &
    197224             exp(-1.*(ld*1E6)*(D*1E-6))*1E-12
    198         scaled = .true.
    199 
    200       endif
    201 
    202       N = fc*rho_a*(Q*1E-3)
    203 
    204     else
    205 
    206 !     // ld will be determined from temperature, then N0 follows
     225     
     226        N = fc*rho_a*(Q*1E-3)
     227
     228    else
     229
     230      !  ld "parameterized" from temperature (carry over from original Quickbeam).
    207231      ld = 1220*10.**(-0.0245*tc)*1E-6
    208232      N0 = ((ld*1E6)**(1+bpm)*Q*1E-3*rho_a)/(apm*gamma(1+bpm))
    … …  
    223247
    224248  case(3)
     249 
     250    if(Re>0) then
     251        print *, 'Variable Re not supported for ', &
     252         'Power-Law distribution'
     253    stop
     254    elseif(Np>0) then
     255        print *, 'Variable Np not supported for ', &
     256         'Power-Law distribution'
     257    stop
     258    endif
    225259
    226260!   :: br parameter
    … …  
    257291
    258292!   :: commented lines are original method with constant density
    259       ! rc = 500.               ! (kg/m^3)
     293      ! rc = 500.       ! (kg/m^3)
    260294      ! tmp1 = 6*rho_a*(bhp+4)
    261295      ! tmp2 = pi*rc*(dmax_mm**(bhp+4))*(1-(dmin_mm/dmax_mm)**(bhp+4))
    … …  
    273307      do k=lidx,uidx
    274308 
    275         N(k) = ( &
     309        N(k) = ( &
    276310        ahp*(D(k)*1E-3)**bhp &
    277         ) * 1E-12   
     311    ) * 1E-12   
    278312
    279313      enddo
    280314
    281         ! print *,'test=',ahp,bhp,ahp/(rho_a*Q),D(100),N(100),bpm,dmin_mm,dmax_mm
     315    ! print *,'test=',ahp,bhp,ahp/(rho_a*Q),D(100),N(100),bpm,dmin_mm,dmax_mm
    282316
    283317! ---------------------------------------------------------!
    … …  
    288322  case(4)
    289323 
    290     if (scaled .eqv. .false.) then
    291    
    292       D0 = p1
     324    if (Re>0) then
     325        D0 = Re
     326    else
     327        D0 = p1
     328    endif
     329   
    293330      rho_e = (6/pi)*apm*(D0*1E-6)**(bpm-3)
    294331      fc(1) = (6./(pi*D0**3*rho_e))*1E12
    295       scaled = .true.
    296      
    297     endif
    298    
    299     N(1) = fc(1)*rho_a*(Q*1E-3)
     332      N(1) = fc(1)*rho_a*(Q*1E-3)
    300333   
    301334! ---------------------------------------------------------!
    … …  
    308341
    309342  case(5)
    310     if (abs(p1+1) < 1E-8) then
     343    if (abs(p1+1) < 1E-8 .or. Re>0 ) then
    311344
    312345!     // rg, log_sigma_g are given
    … …  
    314347      tmp2 = (bpm*log_sigma_g)**2.
    315348      if(Re.le.0) then
    316         rg = p2
    317       else
    318         rg =Re*exp(-2.5*(log_sigma_g**2))
     349        rg = p2
     350      else
     351    rg =Re*exp(-2.5*(log_sigma_g**2))
    319352      endif
    320353 
    321       if (scaled .eqv. .false.) then
    322            
    323         fc = 0.5 * ( &
    324              (1./((2.*rg*1E-6)**(bpm)*apm*(2.*pi)**(0.5) * &
    325              log_sigma_g*D*0.5*1E-6)) * &
    326              exp(-0.5*((log(0.5*D/rg)/log_sigma_g)**2.+tmp2)) &
    327              ) * 1E-12
    328         scaled = .true.
    329              
     354         fc = 0.5 * ( &
     355         (1./((2.*rg*1E-6)**(bpm)*apm*(2.*pi)**(0.5) * &
     356         log_sigma_g*D*0.5*1E-6)) * &
     357         exp(-0.5*((log(0.5*D/rg)/log_sigma_g)**2.+tmp2)) &
     358         ) * 1E-12
     359               
     360      N = fc*rho_a*(Q*1E-3)
     361     
     362    elseif (abs(p2+1) < 1E-8 .or. Np>0) then
     363
     364!     // Np, log_sigma_g are given   
     365      if(Np>0) then
     366        local_Np=Np
     367      else
     368        local_Np = p1
    330369      endif
    331                
    332       N = fc*rho_a*(Q*1E-3)
    333      
    334     elseif (abs(p2+1) < 1E-8) then
    335 
    336 !     // N0, log_sigma_g are given   
    337       Np = p1
     370     
    338371      log_sigma_g = p3
    339       N0 = np*rho_a
     372      N0 = local_np*rho_a
    340373      tmp1 = (rho_a*(Q*1E-3))/(2.**bpm*apm*N0)
    341374      tmp2 = exp(0.5*bpm**2.*(log_sigma_g))**2.     
    … …  
    344377      N = 0.5*( &
    345378        N0 / ((2.*pi)**(0.5)*log_sigma_g*D*0.5*1E-6) * &
    346         exp((-0.5*(log(0.5*D/rg)/log_sigma_g)**2.)) &
    347         ) * 1E-12     
    348      
    349     else
    350 
    351 !     // vu isn't given
     379    exp((-0.5*(log(0.5*D/rg)/log_sigma_g)**2.)) &
     380    ) * 1E-12     
     381     
     382    else
     383
    352384      print *, 'Error: Must specify a value for sigma_g'
    353385      stop
Note: See TracChangeset for help on using the changeset viewer.