run/create_p3_lookupTable_2.f90-v5.3

   1 PROGRAM create_p3_lookuptable_2
   2
   3 !______________________________________________________________________________________
   4 !
   5 ! This program creates the lookup table 'p3_lookupTable_2.dat' for inter-category ice-ice
   6 ! interactions for the multi-ice-category configuration of the P3 microphysics scheme.
   7 !
   8 !--------------------------------------------------------------------------------------
   9 ! Version:       5.3
  10 ! Last modified: 2022 June
  11 !______________________________________________________________________________________
  12
  13 !______________________________________________________________________________________
  14 !
  15 ! To generate 'p3_lookupTable_2.dat' using this code, do the following steps :
  16 !
  17 ! 1. Break up this code into two parts (-top.f90 and -bottom.90).  Search for the string
  18 !    'RUNNING IN PARALLEL MODE' and follow instructions.  (In the future, a script
  19 !    will be written to automate this.)
  20 !
  21 ! 2. Copy the 3 pieces of text below to create indivudual executable shell scripts.
  22
  23 ! 3. Run script 1 (./go_1-compile.ksh).  This script will recreate several
  24 !    versions of the full code, concatenating the -top.f90 and the -bottom.90 with
  25 !    the following code (e.g.) in between (looping through values of i1:
  26 !
  27 !    i1  = 1
  28 !
  29 !    Each version of full_code.f90 is then compiled, with a unique executable name.
  30 !    Note, this is done is place the outer i1 loop in order to parallelized .
  31 !
  32 ! 4. Run script 2 (./go_2-submit.csh)  This create temporary work directories,
  33 !    moves each executable to each work directory, and runs the executables
  34 !    simultaneously.  (Note, it is assumed that the machine on which this is done
  35 !    has multiple processors, though it is not necessary.)
  36 !
  37 ! 5. Run script 3 (./go_3-concatenate.ksh).  This concatenates the individual
  38 !    partial tables (the output in each working directory) into a single, final
  39 !    file 'p3_lookupTable_2.dat'  Once it is confirmed that the table is correct,
  40 !    the work directories and their contents can be removed.
  41 !
  42 !  Note:  For testing or to run in serial, compile with double-precision
  43 !         e.g. ifort -r8 create_p3_lookupTable_2.f90
  44 !              gfortran -fdefault-real-8 create_p3_lookupTable_2.f90
  45 !______________________________________________________________________________________
  46
  47 !--------------------------------------------------------------------------------------------
  48 !# Parallel script 1 (of 3):  [copy text below (uncommented) to file 'go_1-compile.ksh']
  49 !#  - creates individual parallel codes, compiles each
  50
  51 !#!/bin/ksh
  52 !
  53 ! for i1 in 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
  54 ! do
  55 !
  56 !  rm cfg_input full_code.f90
  57 !
  58 !  cat > cfg_input << EOF
  59 !   i1  = ${i1}
  60 !EOF
  61 !
  62 !  cat create_p3_lookupTable_2-top.f90 cfg_input create_p3_lookupTable_2-bottom.f90 > full_code.f90
  63 !
  64 !  echo 'Compiling 'exec_${i1}
  65 !  #pgf90 -r8 full_code.f90
  66 !  #gfortran -fdefault-real-8 full_code.f90
  67 !  ifort -r8 full_code.f90
  68 !  mv a.out exec_${i1}
  69 !
  70 ! done
  71 !
  72 ! rm cfg_input full_code.f90
  73
  74 !--------------------------------------------------------------------------------------------
  75 !# Parallel script 2 (of 3):   [copy text below (uncommented) to file 'go_2-submit.ksh']
  76 !#  - creates individual work directories, launches each executable
  77
  78 !#!/bin/ksh
  79
  80 !for exec in `ls exec_*`
  81 !do
  82 !   echo Submitting: ${exec}
  83 !   mkdir ${exec}-workdir
  84 !   mv ${exec} ${exec}-workdir
  85 !   cd ${exec}-workdir
  86 !   ./${exec} > log &
  87 !   cd ..
  88 !done
  89
  90 !--------------------------------------------------------------------------------------------
  91 !# Parallel script 3 (of 3):   [copy text below (uncommented) to file 'go_3-concatenate.ksh]
  92 !#  - concatenates the output of each parallel job into a single output file.
  93
  94 !#!/bin/ksh
  95 !
  96 ! rm lt_total
  97 !
  98 ! for i in `ls exec*/*dat`
  99 ! do
 100 !    echo $i
 101 !    cat lt_total $i > lt_total_tmp
 102 !    mv lt_total_tmp lt_total
 103 ! done
 104 !
 105 ! mv lt_total p3_lookupTable_2.dat
 106 !
 107 ! echo 'Done.  Work directories and contents can now be removed.'
 108 ! echo 'Be sure to re-name the file with the appropriate extension, with the version number
 109 ! echo 'corresponding to that in the header.  (e.g. 'p3_lookupTable_2.dat-v5.0')'
 110
 111 !--------------------------------------------------------------------------------------------
 112
 113  implicit none
 114
 115  character(len=16), parameter :: version = '5.3'
 116
 117  real    :: pi,g,p,t,rho,mu,pgam,ds,cs,bas,aas,dcrit,eii
 118  integer :: k,ii,jj,kk,dumii,j2
 119
 120  integer, parameter :: n_rhor  =  5
 121  integer, parameter :: n_Fr    =  4
 122  integer, parameter :: n_Qnorm = 25
 123
 124  integer, parameter :: num_bins1 =  1000   ! number of bins for numerical integration of fall speeds and total N
 125 !integer, parameter :: num_bins2 =  9000   ! number of bins for solving PSD parameters  [based on Dm_max = 2000.]
 126  integer, parameter :: num_bins2 = 11000   ! number of bins for solving PSD parameters  [based on Dm_max = 400000.]
 127  real,    parameter :: dd        = 20.e-6  ! width of size bin (m)
 128
 129  integer            :: i_rhor,i_rhor1,i_rhor2  ! indices for rho_rime                 [1 .. n_rhor]
 130  integer            :: i_Fr,i_Fr1,i_Fr2        ! indices for rime-mass-fraction loop  [1 .. n_Fr]
 131  integer            :: i,i1,i2                 ! indices for normalized (by N) Q loop [1 .. n_Qnorm]  (i is i_Qnorm in LT1)
 132
 133  real :: N,q,qdum,dum1,dum2,cs1,ds1,lam,n0,lamf,qerror,del0,c0,c1,c2,sum1,sum2,          &
 134          sum3,sum4,xx,a0,b0,a1,b1,dum,bas1,aas1,aas2,bas2,gammq,d1,d2,delu,lamold,       &
 135          cap,lamr,dia,amg,dv,n0dum,sum5,sum6,sum7,sum8,dg,cg,bag,aag,dcritg,dcrits,      &
 136          dcritr,csr,dsr,duml,dum3,rhodep,cgpold,m1,m2,m3,dt,mur,initlamr,lamv,Fr,        &
 137          rdumii,lammin,lammax,cs2,ds2,intgrR1,intgrR2,intgrR3,intgrR4,q_agg,n_agg
 138
 139  real :: diagnostic_mui ! function to return diagnostic value of shape paramter, mu_i
 140  real :: Q_normalized   ! function
 141  real :: lambdai        ! function
 142
 143  real, parameter :: Dm_max = 40000.e-6   ! max. mean ice [m] size for lambda limiter
 144 !real, parameter :: Dm_max =  2000.e-6   ! max. mean ice [m] size for lambda limiter
 145  real, parameter :: Dm_min =     2.e-6   ! min. mean ice [m] size for lambda limiter
 146
 147  real, parameter :: thrd = 1./3.
 148  real, parameter :: sxth = 1./6.
 149
 150 !real, dimension(n_Qnorm,n_Fr,n_rhor)    :: qsave,nsave,qon1
 151  real, dimension(n_Qnorm,n_Fr,n_rhor)    :: qsave
 152  real, dimension(n_rhor,n_Fr)            :: dcrits1,dcritr1,csr1,dsr1,dcrits2,dcritr2,csr2,dsr2
 153  real, dimension(n_rhor,n_Fr,n_Qnorm)    :: n01,mu_i1,lam1,n02,mu_i2,lam2
 154  real, dimension(n_rhor)                 :: cgp,crp
 155  real, dimension(n_rhor,n_Fr)            :: cgp1,cgp2
 156  real, dimension(n_Fr)                   :: Fr_arr
 157  real, dimension(n_rhor,n_Fr,num_bins1)  :: fall1,fall2
 158  real, dimension(num_bins1)              :: num1,num2
 159
 160  logical, dimension(n_rhor,n_Fr,n_Qnorm) :: log_lamIsMax
 161
 162  character(len=2014) :: filename
 163
 164 !-------------------------------------------------------------------------------
 165
 166 ! set constants and parameters
 167
 168 ! assume 600 hPa, 253 K for p and T for fallspeed calcs (for reference air density)
 169  pi  = acos(-1.)
 170  g   = 9.861                           ! gravity
 171  p   = 60000.                          ! air pressure (pa)
 172  t   = 253.15                          ! temp (K)
 173  rho = p/(287.15*t)                    ! air density (kg m-3)
 174  mu  = 1.496E-6*t**1.5/(t+120.)/rho    ! viscosity of air
 175  dv  = 8.794E-5*t**1.81/p              ! diffusivity of water vapor in air
 176  dt  = 10.
 177
 178 ! parameters for surface roughness of ice particle
 179 ! see mitchell and heymsfield 2005
 180  del0 = 5.83
 181  c0   = 0.6
 182  c1   = 4./(del0**2*c0**0.5)
 183  c2   = del0**2/4.
 184
 185 !--- specified mass-dimension relationship (cgs units) for unrimed crystals:
 186 ! ms = cs*D^ds
 187 !
 188 ! for graupel:
 189 ! mg = cg*D^dg     no longer used, since bulk volume is predicted
 190
 191 ! Heymsfield et al. 2006
 192 !      ds=1.75
 193 !      cs=0.0040157+6.06e-5*(-20.)
 194 ! sector-like branches (P1b)
 195 !      ds=2.02
 196 !      cs=0.00142
 197 ! bullet-rosette
 198 !     ds=2.45
 199 !      cs=0.00739
 200 ! side planes
 201 !      ds=2.3
 202 !      cs=0.00419
 203 ! radiating assemblages of plates (mitchell et al. 1990)
 204 !      ds=2.1
 205 !      cs=0.00239
 206 ! aggreagtes of side planes, bullets, etc. (Mitchell 1996)
 207 !      ds=2.1
 208 !      cs=0.0028
 209 ! Brown and Francis (1995)
 210  ds = 1.9
 211 ! cs = 0.01855 ! original, based on assumption of Dmax
 212  cs = 0.0121 ! scaled value based on assumtion of Dmean from Hogan et al. 2012, JAMC
 213
 214 ! note: if using brown and francis, already in mks units!!!!!
 215 ! uncomment line below if using other snow m-D relationships
 216 !      cs=cs*100.**ds/1000.  ! convert from cgs units to mks
 217 !===
 218
 219 ! applicable for prognostic graupel density
 220 !  note:  cg is not constant, due to variable density
 221  dg = 3.
 222
 223
 224 !--- projected area-diam relationship (mks units) for unrimed crystals:
 225 !       note: projected area = aas*D^bas
 226 ! sector-like branches (P1b)
 227 !      bas = 1.97
 228 !      aas = 0.55*100.**bas/(100.**2)
 229 ! bullet-rosettes
 230 !      bas = 1.57
 231 !      aas = 0.0869*100.**bas/(100.**2)
 232 ! graupel (values for hail)
 233 !      bag=2.0
 234 !      aag=0.625*100.**bag/(100.**2)
 235 ! aggreagtes of side planes, bullets, etc.
 236  bas = 1.88
 237  aas = 0.2285*100.**bas/(100.**2)
 238 !===
 239
 240 !--- projected area-diam relationship (mks units) for graupel:
 241 !      (assumed spheres)
 242 !       note: projected area = aag*D^bag
 243  aag = pi*0.25
 244  bag = 2.
 245 !===
 246
 247 ! calculate critical diameter separating small spherical ice from crystalline ice
 248 ! "Dth" in Morrison and Grabowski 2008
 249
 250 ! !  !open file to write to look-up table (which gets used by P3 scheme)
 251 ! !  open(unit=1,file='./p3_lookup_table_2.dat-v4',status='unknown')
 252
 253 !.........................................................
 254
 255 !dcrit = (pi/(6.*cs)*0.9)**(1./(ds-3.))
 256  dcrit = (pi/(6.*cs)*900.)**(1./(ds-3.))
 257 !dcrit=dcrit/100.  ! convert from cm to m
 258
 259 !.........................................................
 260 ! main loop over graupel density
 261
 262 ! 1D array for RIME density (not ice/graupel density)
 263  crp(1) =  50.*pi*sxth
 264  crp(2) = 250.*pi*sxth
 265  crp(3) = 450.*pi*sxth
 266  crp(4) = 650.*pi*sxth
 267  crp(5) = 900.*pi*sxth
 268
 269 ! array for rime fraction, Fr
 270  Fr_arr(1) = 0.
 271  Fr_arr(2) = 0.333
 272  Fr_arr(3) = 0.667
 273  Fr_arr(4) = 1.
 274
 275 !...........................................................................................
 276
 277 ! parameters for category 1
 278
 279 ! main loop over graupel density
 280
 281  i_rhor_loop_1: do i_rhor = 1,n_rhor
 282
 283 !------------------------------------------------------------------------
 284 ! main loops around N, q, Fr for lookup tables
 285
 286 ! find threshold with rimed mass added
 287
 288 ! loop over rimed mass fraction (4 points)
 289
 290     i_Fr_loop_1: do i_Fr = 1,n_Fr   ! loop for rime mass fraction, Fr
 291
 292       ! Rime mass fraction for the lookup table (specific values in model are interpolated between points)
 293        Fr = Fr_arr(i_Fr)
 294
 295 ! ! ! calculate critical dimension separate graupel and nonspherical ice
 296 ! ! ! "Dgr" in morrison and grabowski (2008)
 297 ! !
 298 ! ! !   dcrits = (cs/crp(i_rhor))**(1./(dg-ds)) ! calculated below for variable graupel density
 299 ! ! !   dcrits = dcrits/100.  ! convert from cm to m
 300 ! !
 301 ! ! ! check to make sure projected area at dcrit not greater than than of solid sphere
 302 ! ! ! stop and give warning message if that is the case
 303 ! !
 304 ! ! !    if (pi/4.*dcrit**2.lt.aas*dcrit**bas) then
 305 ! ! !       print*,'STOP, area > area of solid ice sphere, unrimed'
 306 ! ! !       stop
 307 ! ! !    endif
 308 ! ! !    if (pi/4.*dcrits1(i_rhor,i_Fr)**2.lt.aag*dcrits1(i_rhor,i_Fr)**bag) then
 309 ! ! !       print*,'STOP, area > area of solid ice sphere, graupel'
 310 ! ! !       stop
 311 ! ! !    endif
 312 ! !
 313 ! ! !      cg=cg*100.**dg/1000.  ! convert from cgs units to mks
 314 ! !
 315 ! ! !      print*,cg,dg
 316 ! ! !      stop
 317 ! ! !      do jj=1,100
 318 ! ! !         dd=real(jj)*30.e-6
 319 ! ! !         write(6,'5e15.5')dd,aas*dd**bas,pi/4.*dd**2,
 320 ! ! !     1      cs*dd**ds,pi*sxth*917.*dd**3
 321 ! ! !      end do
 322
 323 ! calculate mass-dimension relationship for partially-rimed crystals
 324 ! msr = csr*D^dsr
 325 ! formula from morrison grabowski 2008
 326
 327 ! dcritr is critical size separating graupel from partially-rime crystal
 328 ! same as "Dcr" in morrison and grabowski 2008
 329
 330 ! first guess, set cgp=crp
 331        cgp1(i_rhor,i_Fr) = crp(i_rhor)
 332
 333 ! case of no riming (Fr = 0%), then we need to set dcrits and dcritr to arbitrary large values
 334
 335        if (i_Fr.eq.1) then
 336
 337           dcrits1(i_rhor,i_Fr) = 1.e+6
 338           dcritr1(i_rhor,i_Fr) = dcrits1(i_rhor,i_Fr)
 339           csr1(i_rhor,i_Fr)    = cs
 340           dsr1(i_rhor,i_Fr)    = ds
 341 ! case of partial riming (Fr between 0 and 100%)
 342
 343        elseif (i_Fr.eq.2.or.i_Fr.eq.3) then
 344
 345           do
 346              dcrits1(i_rhor,i_Fr) = (cs/cgp1(i_rhor,i_Fr))**(1./(dg-ds))
 347              dcritr1(i_rhor,i_Fr) = ((1.+Fr/(1.-Fr))*cs/cgp1(i_rhor,i_Fr))**(1./(dg-ds))
 348              csr1(i_rhor,i_Fr)    = cs*(1.+Fr/(1.-Fr))
 349              dsr1(i_rhor,i_Fr)    = ds
 350 ! get mean density of vapor deposition/aggregation grown ice
 351              rhodep = 1./(dcritr1(i_rhor,i_Fr)-dcrits1(i_rhor,i_Fr))*6.*cs/(pi*(ds-2.))* &
 352                       (dcritr1(i_rhor,i_Fr)**(ds-2.)-dcrits1(i_rhor,i_Fr)**(ds-2.))
 353 ! get graupel density as rime mass fraction weighted rime density plus
 354 ! density of vapor deposition/aggregation grown ice
 355              cgpold = cgp1(i_rhor,i_Fr)
 356              cgp1(i_rhor,i_Fr) = crp(i_rhor)*Fr+rhodep*(1.-Fr)*pi*sxth
 357              if (abs((cgp1(i_rhor,i_Fr)-cgpold)/cgp1(i_rhor,i_Fr)).lt.0.01) goto 115
 358           enddo
 359
 360  115  continue
 361
 362 ! case of complete riming (Fr=100%)
 363        else
 364
 365 ! set threshold size for pure graupel arbitrary large
 366           dcrits1(i_rhor,i_Fr) = (cs/cgp1(i_rhor,i_Fr))**(1./(dg-ds))
 367           dcritr1(i_rhor,i_Fr) = 1.e6
 368           csr1(i_rhor,i_Fr)    = cgp1(i_rhor,i_Fr)
 369           dsr1(i_rhor,i_Fr)    = dg
 370
 371        endif  !if i_Fr.eq.1
 372
 373 !---------------------------------------------------------------------------------------
 374 ! set up particle fallspeed arrays
 375 ! fallspeed is a function of mass dimension and projected area dimension relationships
 376 ! following mitchell and heymsfield (2005), jas
 377
 378 ! set up array of particle fallspeed to make computationally efficient
 379 !.........................................................
 380 ! ****
 381 !  note: this part could be incorporated into the longer (every 2 micron) loop
 382 ! ****
 383        jj_loop_1: do jj = 1,num_bins1
 384
 385           d1 = real(jj)*dd - 0.5*dd   !particle size [m]
 386
 387           if (d1.le.dcrit) then
 388              cs1  = pi*sxth*900.
 389              ds1  = 3.
 390              bas1 = 2.
 391              aas1 = pi/4.
 392           else if (d1.gt.dcrit.and.d1.le.dcrits1(i_rhor,i_Fr)) then
 393              cs1  = cs
 394              ds1  = ds
 395              bas1 = bas
 396              aas1 = aas
 397           else if (d1.gt.dcrits1(i_rhor,i_Fr).and.d1.le.dcritr1(i_rhor,i_Fr)) then
 398              cs1  = cgp1(i_rhor,i_Fr)
 399              ds1  = dg
 400              bas1 = bag
 401              aas1 = aag
 402           else if (d1.gt.dcritr1(i_rhor,i_Fr)) then
 403              cs1  = csr1(i_rhor,i_Fr)
 404              ds1  = dsr1(i_rhor,i_Fr)
 405              if (i_Fr.eq.1) then
 406                 aas1 = aas
 407                 bas1 = bas
 408              else
 409 ! for area, keep bas1 constant, but modify aas1 according to rimed fraction
 410                 bas1 = bas
 411                 dum1 = aas*d1**bas
 412                 dum2 = aag*d1**bag
 413 !               dum3 = (1.-Fr)*dum1+Fr*dum2
 414                 m1   = cs1*d1**ds1
 415                 m2   = cs*d1**ds
 416                 m3   = cgp1(i_rhor,i_Fr)*d1**dg
 417                 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)  !linearly interpolate based on particle mass
 418                 aas1 = dum3/(d1**bas)
 419              endif
 420           endif
 421
 422 ! correction for turbulence
 423 !            if (d1.lt.500.e-6) then
 424           a0 = 0.
 425           b0 = 0.
 426 !            else
 427 !               a0=1.7e-3
 428 !               b0=0.8
 429 !            end if
 430
 431 ! fall speed for ice
 432 ! best number
 433           xx = 2.*cs1*g*rho*d1**(ds1+2.-bas1)/(aas1*(mu*rho)**2)
 434 ! drag terms
 435           b1 = c1*xx**0.5/(2.*((1.+c1*xx**0.5)**0.5-1.)*(1.+c1*xx**0.5)**0.5)-a0*b0*xx** &
 436                b0/(c2*((1.+c1*xx**0.5)**0.5-1.)**2)
 437           a1 = (c2*((1.+c1*xx**0.5)**0.5-1.)**2-a0*xx**b0)/xx**b1
 438 ! velocity in terms of drag terms
 439           fall1(i_rhor,i_Fr,jj) = a1*mu**(1.-2.*b1)*(2.*cs1*g/(rho*aas1))**b1*d1**(b1*(ds1-bas1+2.)-1.)
 440
 441        enddo jj_loop_1
 442
 443 !---------------------------------------------------------------------------------
 444 ! main loops around normalized q for lookup table
 445 !
 446 ! q = normalized ice mass mixing ratio = q/N, units are kg^-1
 447
 448        i_Qnorm_loop_1: do i = 1,n_Qnorm
 449
 450           q = Q_normalized(i)
 451
 452           print*,'i,Fr,i_rhor ',i,i_Fr,i_rhor
 453           print*,'q* ',q
 454
 455           qerror = 1.e+20   !initialize qerror to arbitrarily large value
 456
 457 !.....................................................................................
 458 ! find parameters for gamma distribution
 459
 460 ! size distribution for ice is assumed to be
 461 ! N(D) = n0 * D^pgam * exp(-lam*D)
 462
 463 ! for the given q and N, we need to find n0, pgam, and lam
 464
 465 ! approach for finding lambda:
 466 ! cycle through a range of lambda, find closest lambda to produce correct q
 467
 468 ! start with lam, range of lam from 100 to 5 x 10^6 is large enough to
 469 ! cover full range over mean size from 2 to 5000 micron
 470
 471           ii_loop: do ii = 1,num_bins2
 472
 473              lam1(i_rhor,i_Fr,i) = lambdai(ii)
 474              mu_i1(i_rhor,i_Fr,i) = diagnostic_mui(lam1(i_rhor,i_Fr,i),q,cgp1(i_rhor,i_Fr),Fr,pi)
 475
 476 ! set min,max lam corresponding to Dm_max,Dm_min:
 477              lam1(i_rhor,i_Fr,i) = max(lam1(i_rhor,i_Fr,i),(mu_i1(i_rhor,i_Fr,i)+1.)/Dm_max)
 478              lam1(i_rhor,i_Fr,i) = min(lam1(i_rhor,i_Fr,i),(mu_i1(i_rhor,i_Fr,i)+1.)/Dm_min)
 479 ! get normalized n0 = n0/N
 480              n01(i_rhor,i_Fr,i) = lam1(i_rhor,i_Fr,i)**(mu_i1(i_rhor,i_Fr,i)+1.)/(gamma(mu_i1(i_rhor,i_Fr,i)+1.))
 481
 482 ! calculate integral for each of the 4 parts of the size distribution
 483 ! check difference with respect to q
 484
 485 ! set up m-D relationship for solid ice with D < Dcrit
 486              cs1  = pi*sxth*900.
 487              ds1  = 3.
 488
 489              call intgrl_section(lam1(i_rhor,i_Fr,i),mu_i1(i_rhor,i_Fr,i), ds1,ds,dg,       &
 490                                  dsr1(i_rhor,i_Fr),dcrit,dcrits1(i_rhor,i_Fr),              &
 491                                  dcritr1(i_rhor,i_Fr),intgrR1,intgrR2,intgrR3,intgrR4)
 492
 493 ! intgrR1 is integral from 0 to dcrit (solid ice)
 494 ! intgrR2 is integral from dcrit to dcrits (snow)
 495 ! intgrR3 is integral from dcrits to dcritr (graupel)
 496 ! intgrR4 is integral from dcritr to inf (rimed snow)
 497
 498 ! sum of the integrals from the 4 regions of the size distribution
 499              qdum = n01(i_rhor,i_Fr,i)*(cs1*intgrR1 + cs*intgrR2 + cgp1(i_rhor,i_Fr)*intgrR3 +   &
 500                      csr1(i_rhor,i_Fr)*intgrR4)
 501
 502              if (ii.eq.1) then
 503                 qerror = abs(q-qdum)
 504                 lamf   = lam1(i_rhor,i_Fr,i)
 505              endif
 506
 507 ! find lam with smallest difference between q and estimate of q, assign to lamf
 508              if (abs(q-qdum).lt.qerror) then
 509                 lamf   = lam1(i_rhor,i_Fr,i)
 510                 qerror = abs(q-qdum)
 511              endif
 512
 513           enddo ii_loop
 514
 515 ! check and print relative error in q to make sure it is not too large
 516 ! note: large error is possible if size bounds are exceeded!!!!!!!!!!
 517
 518           print*,'qerror (%)',qerror/q*100.
 519
 520 ! find n0 based on final lam value
 521 ! set final lamf to 'lam' variable
 522 ! this is the value of lam with the smallest qerror
 523           lam1(i_rhor,i_Fr,i) = lamf
 524 ! recalculate mu_i based on final lam
 525           mu_i1(i_rhor,i_Fr,i) = diagnostic_mui(lam1(i_rhor,i_Fr,i),q,cgp1(i_rhor,i_Fr),Fr,pi)
 526
 527 !            n0 = N*lam**(pgam+1.)/(gamma(pgam+1.))
 528
 529 ! find n0 from lam and q
 530 ! this is done instead of finding n0 from lam and N, since N
 531 ! may need to be adjusted to constrain mean size within reasonable bounds
 532
 533           call intgrl_section(lam1(i_rhor,i_Fr,i),mu_i1(i_rhor,i_Fr,i), ds1,ds,dg,        &
 534                                  dsr1(i_rhor,i_Fr),dcrit,dcrits1(i_rhor,i_Fr),            &
 535                                  dcritr1(i_rhor,i_Fr),intgrR1,intgrR2,intgrR3,intgrR4)
 536
 537 ! normalized n0
 538           n01(i_rhor,i_Fr,i) = q/(cs1*intgrR1 + cs*intgrR2 + cgp1(i_rhor,i_Fr)*intgrR3 +       &
 539                                csr1(i_rhor,i_Fr)*intgrR4)
 540           print*,'lam,N0:',lam1(i_rhor,i_Fr,i),n01(i_rhor,i_Fr,i)
 541           print*,'mu_i:',mu_i1(i_rhor,i_Fr,i)
 542           print*,'mean size:',(mu_i1(i_rhor,i_Fr,i)+1.)/lam1(i_rhor,i_Fr,i)
 543
 544 ! At this point, we have solve for all of the size distribution parameters
 545
 546 ! NOTE: In the code it is assumed that mean size and number have already been
 547 ! adjusted, so that mean size will fall within allowed bounds. Thus, we do
 548 ! not apply a lambda limiter here.
 549
 550        enddo i_Qnorm_loop_1
 551
 552     enddo i_Fr_loop_1
 553  enddo i_rhor_loop_1
 554
 555 !--------------------------------------------------------------------
 556 !.....................................................................................
 557 ! now calculate parameters for category 2
 558
 559  i_rhor_loop_2: do i_rhor = 1,n_rhor
 560
 561     i_Fr_loop_2: do i_Fr = 1,n_Fr
 562
 563        Fr = Fr_arr(i_Fr)
 564
 565 ! calculate mass-dimension relationship for partially-rimed crystals
 566 ! msr = csr*D^dsr
 567 ! formula from morrison grabowski 2008
 568
 569 ! dcritr is critical size separating graupel from partially-rime crystal
 570 ! same as "Dcr" in morrison and grabowski 2008
 571
 572 ! first guess, set cgp=crp
 573        cgp2(i_rhor,i_Fr) = crp(i_rhor)
 574
 575 ! case of no riming (Fr = 0%), then we need to set dcrits and dcritr to arbitrary large values
 576
 577        if (i_Fr.eq.1) then
 578
 579           dcrits2(i_rhor,i_Fr) = 1.e+6
 580           dcritr2(i_rhor,i_Fr) = dcrits2(i_rhor,i_Fr)
 581           csr2(i_rhor,i_Fr)    = cs
 582           dsr2(i_rhor,i_Fr)    = ds
 583
 584 ! case of partial riming (Fr between 0 and 100%)
 585        elseif (i_Fr.eq.2.or.i_Fr.eq.3) then
 586
 587           do
 588              dcrits2(i_rhor,i_Fr) = (cs/cgp2(i_rhor,i_Fr))**(1./(dg-ds))
 589              dcritr2(i_rhor,i_Fr) = ((1.+Fr/(1.-Fr))*cs/cgp2(i_rhor,i_Fr))**(1./(dg-ds))
 590              csr2(i_rhor,i_Fr)    = cs*(1.+Fr/(1.-Fr))
 591              dsr2(i_rhor,i_Fr)    = ds
 592 ! get mean density of vapor deposition/aggregation grown ice
 593              rhodep = 1./(dcritr2(i_rhor,i_Fr)-dcrits2(i_rhor,i_Fr))*6.*cs/(pi*(ds-2.))* &
 594                       (dcritr2(i_rhor,i_Fr)**(ds-2.)-dcrits2(i_rhor,i_Fr)**(ds-2.))
 595 ! get graupel density as rime mass fraction weighted rime density plus
 596 ! density of vapor deposition/aggregation grown ice
 597              cgpold = cgp2(i_rhor,i_Fr)
 598              cgp2(i_rhor,i_Fr) = crp(i_rhor)*Fr+rhodep*(1.-Fr)*pi*sxth
 599              if (abs((cgp2(i_rhor,i_Fr)-cgpold)/cgp2(i_rhor,i_Fr)).lt.0.01) goto 116
 600           enddo
 601
 602  116  continue
 603
 604 ! case of complete riming (Fr=100%)
 605        else
 606
 607 ! set threshold size for pure graupel arbitrary large
 608           dcrits2(i_rhor,i_Fr) = (cs/cgp2(i_rhor,i_Fr))**(1./(dg-ds))
 609           dcritr2(i_rhor,i_Fr) = 1.e+6
 610           csr2(i_rhor,i_Fr)    = cgp2(i_rhor,i_Fr)
 611           dsr2(i_rhor,i_Fr)    = dg
 612
 613        endif
 614
 615 !---------------------------------------------------------------------------------------
 616 ! set up particle fallspeed arrays
 617 ! fallspeed is a function of mass dimension and projected area dimension relationships
 618 ! following mitchell and heymsfield (2005), jas
 619
 620 ! set up array of particle fallspeed to make computationally efficient
 621
 622        jj_loop_2: do jj = 1,num_bins1
 623
 624           d1 = real(jj)*dd - 0.5*dd   !particle size [m]
 625
 626           if (d1.le.dcrit) then
 627              cs1  = pi*sxth*900.
 628              ds1  = 3.
 629              bas1 = 2.
 630              aas1 = pi/4.
 631           else if (d1.gt.dcrit.and.d1.le.dcrits2(i_rhor,i_Fr)) then
 632              cs1  = cs
 633              ds1  = ds
 634              bas1 = bas
 635              aas1 = aas
 636           else if (d1.gt.dcrits2(i_rhor,i_Fr).and.d1.le.dcritr2(i_rhor,i_Fr)) then
 637              cs1  = cgp2(i_rhor,i_Fr)
 638              ds1  = dg
 639              bas1 = bag
 640              aas1 = aag
 641           else if (d1.gt.dcritr2(i_rhor,i_Fr)) then
 642              cs1  = csr2(i_rhor,i_Fr)
 643              ds1  = dsr2(i_rhor,i_Fr)
 644              if (i_Fr.eq.1) then
 645                 aas1 = aas
 646                 bas1 = bas
 647              else
 648 ! for area, keep bas1 constant, but modify aas1 according to rimed fraction
 649                 bas1 = bas
 650                 dum1 = aas*d1**bas
 651                 dum2 = aag*d1**bag
 652 !               dum3 = (1.-Fr)*dum1+Fr*dum2
 653                 m1   = cs1*d1**ds1
 654                 m2   = cs*d1**ds
 655                 m3   = cgp2(i_rhor,i_Fr)*d1**dg
 656                 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)  !linearly interpolate based on particle mass
 657                 aas1 = dum3/(d1**bas)
 658              endif
 659           endif
 660
 661 ! correction for turbulence
 662 !            if (d1.lt.500.e-6) then
 663           a0 = 0.
 664           b0 = 0.
 665 !            else
 666 !               a0=1.7e-3
 667 !               b0=0.8
 668 !            end if
 669
 670 ! fall speed for ice
 671 ! best number
 672           xx = 2.*cs1*g*rho*d1**(ds1+2.-bas1)/(aas1*(mu*rho)**2)
 673
 674 ! drag terms
 675           b1 = c1*xx**0.5/(2.*((1.+c1*xx**0.5)**0.5-1.)*(1.+c1*xx**0.5)**0.5)-a0*b0*xx** &
 676                b0/(c2*((1.+c1*xx**0.5)**0.5-1.)**2)
 677
 678           a1 = (c2*((1.+c1*xx**0.5)**0.5-1.)**2-a0*xx**b0)/xx**b1
 679
 680 ! velocity in terms of drag terms
 681           fall2(i_rhor,i_Fr,jj) = a1*mu**(1.-2.*b1)*(2.*cs1*g/(rho*aas1))**b1*d1**(b1*(ds1-bas1+2.)-1.)
 682
 683 !---------------------------------------------------------------
 684        enddo jj_loop_2
 685
 686
 687 !---------------------------------------------------------------------------------
 688
 689 ! q = total ice mixing ratio (vapor dep. plus rime mixing ratios), normalized by N
 690
 691        i_Qnorm_loop_2: do i = 1,n_Qnorm
 692
 693           lamold = 0.
 694
 695           q = Q_normalized(i)
 696
 697           print*,'&&&&&&&&&&&i_rhor',i_rhor
 698           print*,'***************',i,k
 699           print*,'Fr',Fr
 700           print*,'q,N',q,N
 701
 702           qerror = 1.e+20   ! initialize qerror to arbitrarily large value
 703
 704 !.....................................................................................
 705 ! find parameters for gamma distribution
 706
 707 ! size distribution for ice is assumed to be
 708 ! N(D) = n0 * D^pgam * exp(-lam*D)
 709
 710 ! for the given q and N, we need to find n0, pgam, and lam
 711
 712 ! approach for finding lambda:
 713 ! cycle through a range of lambda, find closest lambda to produce correct q
 714
 715 ! start with lam, range of lam from 100 to 5 x 10^6 is large enough to
 716 ! cover full range over mean size from 2 to 5000 micron
 717
 718           ii_loop_2: do ii = 1,num_bins2
 719
 720              lam2(i_rhor,i_Fr,i) = lambdai(ii)
 721              mu_i2(i_rhor,i_Fr,i) = diagnostic_mui(lam2(i_rhor,i_Fr,i),q,cgp2(i_rhor,i_Fr),Fr,pi)
 722
 723             !set min,max lam corresponding to Dm_max, Dm_min
 724              lam2(i_rhor,i_Fr,i) = max(lam2(i_rhor,i_Fr,i),(mu_i2(i_rhor,i_Fr,i)+1.)/Dm_max)
 725              lam2(i_rhor,i_Fr,i) = min(lam2(i_rhor,i_Fr,i),(mu_i2(i_rhor,i_Fr,i)+1.)/Dm_min)
 726
 727             !get n0, note this is normalized
 728              n02(i_rhor,i_Fr,i) = lam2(i_rhor,i_Fr,i)**(mu_i2(i_rhor,i_Fr,i)+1.)/ &
 729                    (gamma(mu_i2(i_rhor,i_Fr,i)+1.))
 730
 731 ! calculate integral for each of the 4 parts of the size distribution
 732 ! check difference with respect to q
 733
 734 ! set up m-D relationship for solid ice with D < Dcrit
 735              cs1  = pi*sxth*900.
 736              ds1  = 3.
 737
 738              call intgrl_section(lam2(i_rhor,i_Fr,i),mu_i2(i_rhor,i_Fr,i), ds1,ds,dg,      &
 739                                  dsr2(i_rhor,i_Fr),dcrit,dcrits2(i_rhor,i_Fr),             &
 740                                  dcritr2(i_rhor,i_Fr),intgrR1,intgrR2,intgrR3,intgrR4)
 741
 742 ! sum of the integrals from the 4 regions of the size distribution
 743              qdum = n02(i_rhor,i_Fr,i)*(cs1*intgrR1 + cs*intgrR2 + cgp2(i_rhor,i_Fr)*intgrR3 +  &
 744                     csr2(i_rhor,i_Fr)*intgrR4)
 745
 746              if (ii.eq.1) then
 747                 qerror = abs(q-qdum)
 748                 lamf   = lam2(i_rhor,i_Fr,i)
 749              endif
 750
 751 ! find lam with smallest difference between q and estimate of q, assign to lamf
 752              if (abs(q-qdum).lt.qerror) then
 753                 lamf   = lam2(i_rhor,i_Fr,i)
 754                 qerror = abs(q-qdum)
 755              endif
 756
 757           enddo ii_loop_2
 758
 759 ! check and print relative error in q to make sure it is not too large
 760 ! note: large error is possible if size bounds are exceeded!!!!!!!!!!
 761
 762           print*,'qerror (%)',qerror/q*100.
 763
 764 ! find n0 based on final lam value
 765 ! set final lamf to 'lam' variable
 766 ! this is the value of lam with the smallest qerror
 767           lam2(i_rhor,i_Fr,i) = lamf
 768
 769 ! recalculate mu based on final lam
 770 !         mu_i2(i_rhor,i_Fr,i) = diagnostic_mui(log_diagmu_orig,lam2(i_rhor,i_Fr,i),q,cgp2(i_rhor,i_Fr),Fr,pi)
 771           mu_i2(i_rhor,i_Fr,i) = diagnostic_mui(lam2(i_rhor,i_Fr,i),q,cgp2(i_rhor,i_Fr),Fr,pi)
 772
 773 !            n0 = N*lam**(pgam+1.)/(gamma(pgam+1.))
 774
 775 ! find n0 from lam and q
 776 ! this is done instead of finding n0 from lam and N, since N
 777 ! may need to be adjusted to constrain mean size within reasonable bounds
 778
 779           call intgrl_section(lam2(i_rhor,i_Fr,i),mu_i2(i_rhor,i_Fr,i), ds1,ds,dg,         &
 780                                  dsr2(i_rhor,i_Fr),dcrit,dcrits2(i_rhor,i_Fr),             &
 781                                  dcritr2(i_rhor,i_Fr),intgrR1,intgrR2,intgrR3,intgrR4)
 782
 783           n02(i_rhor,i_Fr,i) = q/(cs1*intgrR1 + cs*intgrR2 + cgp2(i_rhor,i_Fr)*intgrR3 +        &  ! n0 is normalized
 784                                csr2(i_rhor,i_Fr)*intgrR4)
 785           print*,'lam,N0:',lam2(i_rhor,i_Fr,i),n02(i_rhor,i_Fr,i)
 786           print*,'pgam:',mu_i2(i_rhor,i_Fr,i)
 787
 788           qsave(i,i_Fr,i_rhor) = q      ! q is normalized (Q/N)
 789
 790           log_lamIsMax(i_rhor,i_Fr,i) = abs(lam2(i_rhor,i_Fr,i)-lamold) .lt. 1.e-8
 791           lamold = lam2(i_rhor,i_Fr,i)
 792
 793        enddo i_Qnorm_loop_2
 794
 795     enddo i_Fr_loop_2
 796  enddo i_rhor_loop_2
 797
 798
 799 ! At this point, we have solve for all of the size distribution parameters
 800
 801 ! NOTE: In the code it is assumed that mean size and number have already been
 802 ! adjusted, so that mean size will fall within allowed bounds. Thus, we do
 803 ! not apply a lambda limiter here.
 804
 805 !--------------------------------------------------------------------
 806
 807 !                            RUNNING IN PARALLEL MODE:
 808 !
 809 !------------------------------------------------------------------------------------
 810 ! CODE ABOVE HERE IS FOR THE "TOP" OF THE BROKEN UP CODE (for running in parallel)
 811 !
 812 !   Before running ./go_1-compile.ksh, delete all lines below this point and
 813 !   and save as 'create_p3_lookupTable_2-top.f90'
 814 !------------------------------------------------------------------------------------
 815
 816 ! For testing single values, uncomment the following:
 817 ! i1 = 1
 818
 819 !------------------------------------------------------------------------------------
 820 ! CODE BELOW HERE IS FOR THE "BOTTOM" OF THE BROKEN UP CODE (for running in parallel)
 821 !
 822 !   Before running ./go_1-compile.ksh, delete all lines below this point and
 823 !   and save as 'create_p3_lookupTable_2-bottom.f90'
 824 !------------------------------------------------------------------------------------
 825
 826 !.....................................................................................
 827 ! begin category process interaction calculations for the lookup table
 828
 829 !.....................................................................................
 830 ! collection of category 1 by category 2
 831 !.....................................................................................
 832
 833  write (filename, "(A12,I0.2,A4)") "lookupTable_2-",i1,".dat"
 834  filename = trim(filename)
 835  open(unit=1, file=filename, status='unknown')
 836
 837  !header:
 838  if (i1==1) then
 839     write(1,*) 'LOOKUP_TABLE_2-version:  ',trim(version)
 840     write(1,*)
 841  endif
 842
 843
 844  ! Note: i1 loop (do/enddo statements) is commented out for parallelization; i1 gets initizatized there
 845  ! - to run in serial, uncomment the 'Qnorm_loop_3: do i1' statement and the corresponding 'enddo'
 846
 847 !Qnorm_loop_3: do i1 = 1,n_Qnorm    ! COMMENTED OUT FOR PARALLELIZATION
 848    do i_Fr1 = 1,n_Fr
 849      do i_rhor1 = 1,n_rhor
 850        do i2 = 1,n_Qnorm
 851          do i_Fr2 = 1,n_Fr
 852            do i_rhor2 = 1,n_rhor
 853
 854                 lamIsMax: if (.not. log_lamIsMax(i_rhor2,i_Fr2,i2)) then
 855
 856                 sum1 = 0.
 857                 sum2 = 0.
 858
 859               !set up binned distribution of ice from categories 1 and 2 (distributions normalized by N)
 860                 do jj = 1,num_bins1
 861                    d1 = real(jj)*dd - 0.5*dd
 862                    num1(jj) = n01(i_rhor1,i_Fr1,i1)*d1**mu_i1(i_rhor1,i_Fr1,i1)*         &
 863                               exp(-lam1(i_rhor1,i_Fr1,i1)*d1)*dd
 864                    num2(jj) = n02(i_rhor2,i_Fr2,i2)*d1**mu_i2(i_rhor2,i_Fr2,i2)*         &
 865                               exp(-lam2(i_rhor2,i_Fr2,i2)*d1)*dd
 866                 enddo
 867
 868 ! loop over size distribution
 869 ! note: collection of ice within the same bin is neglected
 870
 871 ! loop over particle 1
 872                 jj_loop_3: do jj = num_bins1,1,-1
 873
 874                    d1 = real(jj)*dd - 0.5*dd   !particle size [m]
 875
 876                    if (d1.le.dcrit) then
 877                       cs1  = pi*sxth*900.
 878                       ds1  = 3.
 879                       bas1 = 2.
 880                       aas1 = pi*0.25
 881                    elseif (d1.gt.dcrit.and.d1.le.dcrits1(i_rhor1,i_Fr1)) then
 882                       cs1  = cs
 883                       ds1  = ds
 884                       bas1 = bas
 885                       aas1 = aas
 886                    else if (d1.gt.dcrits1(i_rhor1,i_Fr1).and.d1.le.dcritr1(i_rhor1,i_Fr1)) then
 887                       cs1  = cgp1(i_rhor1,i_Fr1)
 888                       ds1  = dg
 889                       bas1 = bag
 890                       aas1 = aag
 891                    else if (d1.gt.dcritr1(i_rhor1,i_Fr1)) then
 892                       cs1 = csr1(i_rhor1,i_Fr1)
 893                       ds1 = dsr1(i_rhor1,i_Fr1)
 894                       if (i_Fr1.eq.1) then
 895                          aas1 = aas
 896                          bas1 = bas
 897                       else
 898 ! for area, keep bas1 constant, but modify aas1 according to rimed fraction
 899                          bas1 = bas
 900                          dum1 = aas*d1**bas
 901                          dum2 = aag*d1**bag
 902                          m1   = cs1*d1**ds1
 903                          m2   = cs*d1**ds
 904                          m3   = cgp1(i_rhor1,i_Fr1)*d1**dg
 905 ! linearly interpolate based on particle mass
 906                          dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
 907                          aas1 = dum3/(d1**bas)
 908                       endif
 909                    endif
 910
 911 ! loop over particle 2
 912                    kk_loop: do kk = num_bins1,1,-1
 913
 914                       d2 = real(kk)*dd - 0.5*dd   !particle size [m]
 915
 916 ! parameters for particle 2
 917                       if (d2.le.dcrit) then
 918 !                        cs2  = pi*sxth*900.
 919 !                        ds2  = 3.
 920                          bas2 = 2.
 921                          aas2 = pi*0.25
 922                       elseif (d2.gt.dcrit.and.d2.le.dcrits2(i_rhor2,i_Fr2)) then
 923 !                        cs2  = cs
 924 !                        ds2  = ds
 925                          bas2 = bas
 926                          aas2 = aas
 927                       else if (d2.gt.dcrits2(i_rhor2,i_Fr2).and.d2.le.dcritr2(i_rhor2,i_Fr2)) then
 928 !                        cs2  = cgp2(i_rhor2,i_Fr2)
 929 !                        ds2  = dg
 930                          bas2 = bag
 931                          aas2 = aag
 932                       else if (d2.gt.dcritr2(i_rhor2,i_Fr2)) then
 933                          cs2 = csr2(i_rhor2,i_Fr2)
 934                          ds2 = dsr2(i_rhor2,i_Fr2)
 935                          if (i_Fr2.eq.1) then
 936                             aas2 = aas
 937                             bas2 = bas
 938                          else
 939 ! for area, keep bas1 constant, but modify aas1 according to rimed fraction
 940                             bas2 = bas
 941                             dum1 = aas*d2**bas
 942                             dum2 = aag*d2**bag
 943                             m1   = cs2*d2**ds2
 944                             m2   = cs*d2**ds
 945                             m3   = cgp2(i_rhor2,i_Fr2)*d2**dg
 946 ! linearly interpolate based on particle mass
 947                             dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
 948                             aas2 = dum3/(d2**bas)
 949                          endif
 950                       endif
 951
 952 ! absolute value, differential fallspeed
 953 !                   delu = abs(fall2(i_rhor2,i_Fr2,kk)-fall1(i_rhor1,i_Fr1,jj))
 954
 955 ! calculate collection of category 1 by category 2, which occurs
 956 ! in fallspeed of particle in category 2 is greater than category 1
 957
 958                       if (fall2(i_rhor2,i_Fr2,kk).gt.fall1(i_rhor1,i_Fr1,jj)) then
 959
 960                          delu = fall2(i_rhor2,i_Fr2,kk)-fall1(i_rhor1,i_Fr1,jj)
 961
 962 ! note: in micro code we have to multiply by air density
 963 ! correction factor for fallspeed, and collection efficiency
 964
 965 ! sum for integral
 966
 967 ! sum1 = # of collision pairs
 968 ! the assumption is that each collision pair reduces crystal
 969 ! number mixing ratio by 1 kg^-1 s^-1 per kg/m^3 of air (this is
 970 ! why we need to multiply by air density, to get units of
 971 ! 1/kg^-1 s^-1)
 972 ! NOTE: For consideration of particle depletion, air density is assumed to be 1 kg m-3
 973 ! This problem could be avoided by using number concentration instead of number mixing ratio
 974 ! for the lookup table calculations, and then not multipling process rate by air density
 975 ! in the P3 code... TO BE FIXED IN THE FUTURE
 976
 977 !                   sum1 = sum1+min((aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk),   &
 978 !                          num(kk)/dt)
 979
 980 ! set collection efficiency
 981 !                  eii = 0.1
 982
 983 ! accretion of number
 984                          sum1 = sum1 + (aas1*d1**bas1 + aas2*d2**bas2)*delu*num1(jj)*num2(kk)
 985 ! accretion of mass
 986                          sum2 = sum2 + cs1*d1**ds1*(aas1*d1**bas1 + aas2*d2**bas2)*delu*num1(jj)*num2(kk)
 987
 988 ! remove collected particles from distribution over time period dt, update num1
 989 !  note -- dt is time scale for removal, not necessarily the model time step
 990
 991                       endif ! fall2(i_rhor2,i_Fr2,kk) > fall1(i_rhor1,i_Fr1,jj)
 992
 993                    enddo kk_loop
 994                 enddo jj_loop_3
 995
 996 ! save for output
 997                 n_agg = sum1
 998                 q_agg = sum2
 999
1000              else
1001
1002                 print*,'&&&&&&, skip'
1003                 n_agg = -999.
1004                 q_agg = -999.
1005
1006              endif lamIsMax
1007
1008              n_agg = dim(n_agg, 1.e-50)
1009              q_agg = dim(q_agg, 1.e-50)
1010
1011             !write(6,'(a5,6i5,2e15.5)') 'index:',i1,i_Fr1,i_rhor1,i2,i_Fr2,i_rhor2,n_agg,q_agg
1012              write(1,'(6i5,2e15.6)') i1,i_Fr1,i_rhor1,i2,i_Fr2,i_rhor2,n_agg,q_agg
1013
1014
1015            enddo   ! i_rhor2 loop
1016          enddo   ! i_Fr2 loop
1017        enddo   ! i2 loop
1018      enddo   ! i_rhor1 loop
1019    enddo   ! i_Fr1
1020 !enddo Qnorm_loop_3  ! i1 loop  (Qnorm)     ! COMMENTED OUT FOR PARALLELIZATION
1021
1022  close(1)
1023
1024 END PROGRAM create_p3_lookuptable_2
1025 !______________________________________________________________________________________
1026
1027 ! Incomplete gamma function
1028 ! from Numerical Recipes in Fortran 77: The Art of Scientific Computing
1029
1030       function gammq(a,x)
1031
1032       real a,gammq,x
1033
1034 ! USES gcf,gser
1035 ! Returns the incomplete gamma function Q(a,x) = 1-P(a,x)
1036
1037       real gammcf,gammser,gln
1038       if (x.lt.0..or.a.le.0) print*, 'bad argument in gammq'
1039       if (x.lt.a+1.) then
1040          call gser(gamser,a,x,gln)
1041          gammq=1.-gamser
1042       else
1043          call gcf(gammcf,a,x,gln)
1044          gammq=gammcf
1045       end if
1046       return
1047       end
1048
1049 !-------------------------------------
1050
1051       subroutine gser(gamser,a,x,gln)
1052       integer itmax
1053       real a,gamser,gln,x,eps
1054       parameter(itmax=100,eps=3.e-7)
1055       integer n
1056       real ap,del,sum,gamma
1057       gln = log(gamma(a))
1058       if (x.le.0.) then
1059          if (x.lt.0.) print*, 'x < 0 in gser'
1060          gamser = 0.
1061          return
1062       end if
1063       ap=a
1064       sum=1./a
1065       del=sum
1066       do n=1,itmax
1067          ap=ap+1.
1068          del=del*x/ap
1069          sum=sum+del
1070          if (abs(del).lt.abs(sum)*eps) goto 1
1071       end do
1072       print*, 'a too large, itmax too small in gser'
1073  1    gamser=sum*exp(-x+a*log(x)-gln)
1074       return
1075       end
1076
1077 !-------------------------------------
1078
1079       subroutine gcf(gammcf,a,x,gln)
1080       integer itmax
1081       real a,gammcf,gln,x,eps,fpmin
1082       parameter(itmax=100,eps=3.e-7,fpmin=1.e-30)
1083       integer i
1084       real an,b,c,d,del,h,gamma
1085       gln=log(gamma(a))
1086       b=x+1.-a
1087       c=1./fpmin
1088       d=1./b
1089       h=d
1090       do i=1,itmax
1091          an=-i*(i-a)
1092          b=b+2.
1093          d=an*d+b
1094          if(abs(d).lt.fpmin) d=fpmin
1095          c=b+an/c
1096          if(abs(c).lt.fpmin) c=fpmin
1097          d=1./d
1098          del=d*c
1099          h = h*del
1100          if(abs(del-1.).lt.eps)goto 1
1101       end do
1102 !     pause 'a too large, itmax too small in gcf'
1103       print*, 'a too large, itmax too small in gcf'
1104  1    gammcf=exp(-x+a*log(x)-gln)*h
1105       return
1106       end
1107
1108 !______________________________________________________________________________________
1109
1110  real function diagnostic_mui(lam,q,cgp,Fr,pi)
1111
1112 !----------------------------------------------------------!
1113 ! Compute mu_i diagnostically.
1114 !----------------------------------------------------------!
1115
1116  implicit none
1117
1118 !Arguments:
1119  real :: lam,q,cgp,Fr,pi
1120
1121 ! Local variables:
1122  real            :: mu_i,dum1,dum2,dum3
1123  real, parameter :: mu_i_min =  0.
1124  real, parameter :: mu_i_max = 20.  !note: for orig 2-mom, mu_i_max=6.
1125  real, parameter :: Di_thres = 0.2  !diameter threshold [mm]
1126
1127 !-- original formulation: (from Heymsfield, 2003)
1128 !  mu_i = 0.076*(lam/100.)**0.8-2.   ! /100 is to convert m-1 to cm-1
1129 !  mu_i = max(mu_i,0.)  ! make sure mu_i >= 0, otherwise size dist is infinity at D = 0
1130 !  mu_i = min(mu_i,6.)
1131
1132
1133 !-- formulation based on 3-moment results (see 2021 JAS article)
1134  dum1 = (q/cgp)**(1./3)*1000.              ! estimated Dmvd [mm], assuming spherical
1135  if (dum1<=Di_thres) then
1136     !diagnostic mu_i, original formulation: (from Heymsfield, 2003)
1137     mu_i = 0.076*(lam*0.01)**0.8-2.        ! /100 is to convert m-1 to cm-1
1138     mu_i = min(mu_i,6.)
1139  else
1140     dum2 = (6./pi)*cgp                     ! mean density (total)
1141     dum3 = max(1., 1.+0.00842*(dum2-400.)) ! adjustment factor for density
1142     mu_i = 0.25*(dum1-Di_thres)*dum3*Fr
1143  endif
1144  mu_i = max(mu_i, mu_i_min)
1145  mu_i = min(mu_i, mu_i_max)
1146
1147  diagnostic_mui = mu_i
1148
1149  end function diagnostic_mui
1150
1151 !______________________________________________________________________________________
1152
1153  subroutine intgrl_section(lam,mu, d1,d2,d3,d4, Dcrit1,Dcrit2,Dcrit3,    &
1154                            intsec_1,intsec_2,intsec_3,intsec_4)
1155  !-----------------
1156  ! Computes and returns partial integrals (partial moments) of ice PSD.
1157  !-----------------
1158
1159  implicit none
1160
1161 !Arguments:
1162  real, intent(in)  :: lam,mu, d1,d2,d3,d4, Dcrit1,Dcrit2,Dcrit3
1163  real, intent(out) :: intsec_1,intsec_2,intsec_3,intsec_4
1164
1165 !Local:
1166  real :: dum,gammq
1167 !-----------------
1168
1169  !Region I -- integral from 0 to Dcrit1  (small spherical ice)
1170  intsec_1 = lam**(-d1-mu-1.)*gamma(mu+d1+1.)*(1.-gammq(mu+d1+1.,Dcrit1*lam))
1171
1172  !Region II -- integral from Dcrit1 to Dcrit2  (non-spherical unrimed ice)
1173  intsec_2 = lam**(-d2-mu-1.)*gamma(mu+d2+1.)*(gammq(mu+d2+1.,Dcrit1*lam))
1174  dum      = lam**(-d2-mu-1.)*gamma(mu+d2+1.)*(gammq(mu+d2+1.,Dcrit2*lam))
1175  intsec_2 = intsec_2-dum
1176
1177  !Region III -- integral from Dcrit2 to Dcrit3  (fully rimed spherical ice)
1178  intsec_3 = lam**(-d3-mu-1.)*gamma(mu+d3+1.)*(gammq(mu+d3+1.,Dcrit2*lam))
1179  dum      = lam**(-d3-mu-1.)*gamma(mu+d3+1.)*(gammq(mu+d3+1.,Dcrit3*lam))
1180  intsec_3 = intsec_3-dum
1181
1182  !Region IV -- integral from Dcrit3 to infinity  (partially rimed ice)
1183  intsec_4 = lam**(-d4-mu-1.)*gamma(mu+d4+1.)*(gammq(mu+d4+1.,Dcrit3*lam))
1184
1185  return
1186
1187  end subroutine intgrl_section
1188
1189
1190 !______________________________________________________________________________________
1191
1192  real function Q_normalized(i_qnorm)
1193
1194  !-----------------
1195  ! Computes the normalized Q (qitot/nitot) based on a given index (used in the Qnorm loops)
1196  !-----------------
1197
1198  implicit none
1199
1200 !arguments:
1201  integer :: i_qnorm
1202
1203 !Q_normalized = 261.7**((i_qnorm+5)*0.2)*1.e-18     ! v4, v5.0 (range of mean mass diameter from ~ 1 micron to 1 cm)   (for n_Qnorm = 25)
1204 !Q_normalized = 800.**((i_qnorm+10)*0.1)*1.e-18     ! based on LT1-v5.2 (range for Dm_max = 400000.)                   (for n_Qnorm = 50)
1205  Q_normalized = 800.**(0.2*(i_qnorm+5))*1.e-18      ! (range for Dm_max = 400000.)                                     (for n_Qnorm = 25)
1206
1207  !--- from LT1:  [TO BE DELTED]
1208  !       !q = 261.7**((i_Qnorm+10)*0.1)*1.e-18     ! old (strict) lambda limiter
1209  !         q = 800.**((i_Qnorm+10)*0.1)*1.e-18     ! new lambda limiter               (for n_Qnorm = 50)
1210  !===
1211
1212  return
1213
1214  end function Q_normalized
1215
1216
1217 !______________________________________________________________________________________
1218
1219  real function lambdai(i)
1220
1221  !-----------------
1222  ! Computes the lambda_i for a given index (used in loops to solve for PSD parameters)
1223  !-----------------
1224
1225  implicit none
1226
1227 !arguments:
1228  integer, intent(in) :: i
1229
1230 !lambdai = 1.0013**i*100.   !used with Dm_max =   2000.
1231  lambdai = 1.0013**i*10.    !used with Dm_max = 400000.
1232
1233  return
1234
1235  end function lambdai
1236