forked from ESCOMP/CARMA_base
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcarma_mod.F90
1514 lines (1309 loc) · 59.5 KB
/
carma_mod.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!! The CARMA module contains an interface to the Community Aerosol and Radiation
!! Model for Atmospheres (CARMA) bin microphysical model [Turco et al. 1979;
!! Toon et al. 1988]. This implementation has been customized to work within
!! other model frameworks, so although it can be provided with an array of
!! columns, it does not do horizontal transport and just does independent 1-D
!! calculations upon each column.
!!
!! The typical usage for the CARMA and CARMASTATE objects within a model would be:
!!>
!! ! This first section of code is done during the parent model's initialzation,
!! ! and there should be a unique CARMA object created for each thread of
!! ! execution.
!!
!! ! Create the CARMA object.
!! call CARMA_Create(carma, ...)
!!
!! ! Define the microphysical components.
!! call CARMAGROUP_Create(carma, ...) ! One or more calls
!!
!! call CARMAELEMENT_Create(carma, ...) ! One or more calls
!!
!! call CARMASOLUTE_Create(carma, ...) ! Zero or more calls
!!
!! call CARMAGAS_Create(carma, ...) ! Zero or more calls
!!
!! ! Define the relationships for the microphysical processes.
!! call CARMA_AddCoagulation(carma, ...) ! Zero or more calls
!! call CARMA_AddGrowth(carma, ...) ! Zero or more calls
!! call CARMA_AddNucleation(carma, ...) ! Zero or more calls
!!
!! ! Initialize things that are state and timestep independent.
!! call CARMA_Initialize(carma, ...)
!!
!! ...
!!
!! ! This section of code is within the parent model's timing loop.
!! !
!! ! NOTE: If using OPEN/MP, then each thread will execute one of
!! ! of these loops per column of data. To avoid having to destroy
!! ! the CARMASTATE object, a pool of CARMASTATE objects could be
!! ! created so that there is one per thread and then the
!! ! CARMA_Destroy() could be called after all columns have been
!! ! processed.
!!
!! ! Initialize CARMA for this model state and timestep.
!! call CARMASTATE_Create(cstate, carma, ...)
!!
!! ! Set the model state for each bin and gas.
!! call CARMASTATE_SetBin(cstate, ...) ! One call for each bin
!! call CARMASTATE_SetGas(cstate, ...) ! One call for each gas
!!
!! ! Calculate the new state
!! call CARMASTATE_Step(cstate, ...)
!!
!! ! Get the results to return back to the parent model.
!! call CARMASTATE_GetBin(cstate, ...) ! One call for each Bin
!! call CARMASTATE_GetGas(cstate, ...) ! One call for each gas
!! call CARMASTATE_GetState(cstate, ...) ! Zero or one calls
!!
!! ! (optional) Deallocate arrays that are not needed beyond this timestep.
!! call CARMASTATE_Destroy(cstate)
!!
!! ...
!!
!! ! This section of code is done during the parent model's cleanup.
!!
!! ! Deallocate all arrays.
!! call CARMA_Destroy(carma)
!!<
!!
!! @version Feb-2009
!! @author Chuck Bardeen, Pete Colarco, Jamie Smith
!
! NOTE: Documentation for this code can be generated automatically using f90doc,
! which is freely available from:
! http://erikdemaine.org/software/f90doc/
! Comment lines with double comment characters are processed by f90doc, and there are
! some special characters added to the comments to control the documentation process.
! In addition to the special characters mentioned in the f990doc documentation, html
! formatting tags (e.g. <i></i>, <sup></sup>, ...) can also be added to the f90doc
! comments.
module carma_mod
! This module maps the parents models constants into the constants need by CARMA. NOTE: CARMA
! constants are in CGS units, while the parent models are typically in MKS units.
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
! CARMA explicitly declares all variables.
implicit none
! All CARMA variables and procedures are private except those explicitly declared to be public.
private
! Declare the public methods.
public CARMA_AddCoagulation
public CARMA_AddGrowth
public CARMA_AddNucleation
public CARMA_Create
public CARMA_Destroy
public CARMA_Get
public CARMA_Initialize
contains
! These are the methods that provide the interface between the parent model and the CARMA
! microphysical model. There are many other methods that are not in this file that are
! used to implement the microphysical calculations needed by the CARMA model. These other
! methods are in effect private methods of the CARMA module, but are in individual files
! since that is the way that CARMA has traditionally been structured and where users may
! want to extend or replace code to affect the microphysics.
!! Creates the CARMA object and allocates arrays to store configuration information
!! that will follow from the CARMA_AddXXX() methods. When the CARMA object is no longer
!! needed, the CARMA_Destroy() method should be used to clean up any allocations
!! that have happened. If LUNOPRT is specified, then the logical unit should be open and
!! ready for output. The caller is responsible for closing the LUNOPRT logical unit
!! after the CARMA object has been destroyed.
!!
!! @version Feb-2009
!! @author Chuck Bardeen
subroutine CARMA_Create(carma, NBIN, NELEM, NGROUP, NSOLUTE, NGAS, NWAVE, rc, &
LUNOPRT, wave, dwave, do_wave_emit, NREFIDX)
type(carma_type), intent(out) :: carma !! the carma object
integer, intent(in) :: NBIN !! number of radius bins per group
integer, intent(in) :: NELEM !! total number of elements
integer, intent(in) :: NGROUP !! total number of groups
integer, intent(in) :: NSOLUTE !! total number of solutes
integer, intent(in) :: NGAS !! total number of gases
integer, intent(in) :: NWAVE !! number of wavelengths
integer, intent(out) :: rc !! return code, negative indicates failure
integer, intent(in), optional :: LUNOPRT !! logical unit number for output
real(kind=f), intent(in), optional :: wave(NWAVE) !! wavelength centers (cm)
real(kind=f), intent(in), optional :: dwave(NWAVE) !! wavelength width (cm)
logical, intent(in), optional :: do_wave_emit(NWAVE) !! do emission in band?
integer, intent(in), optional :: NREFIDX !! number of refractive indices per wavelength
! Local Varaibles
integer :: ier
! Assume success.
rc = RC_OK
! Save off the logic unit used for output if one was provided. If one was provided,
! then assume that CARMA can print output.
if (present(LUNOPRT)) then
carma%f_LUNOPRT = LUNOPRT
carma%f_do_print = .TRUE.
end if
! Save the defintion of the number of comonents involved in the microphysics.
carma%f_NGROUP = NGROUP
carma%f_NELEM = NELEM
carma%f_NBIN = NBIN
carma%f_NGAS = NGAS
carma%f_NSOLUTE = NSOLUTE
carma%f_NWAVE = NWAVE
if (present(NREFIDX)) then
carma%f_NREFIDX = NREFIDX
else
carma%f_NREFIDX = 0
end if
! Allocate tables for the groups.
allocate( &
carma%f_group(NGROUP), &
carma%f_icoag(NGROUP, NGROUP), &
carma%f_inucgas(NGROUP), &
stat=ier)
if(ier /= 0) then
if (carma%f_do_print) write(carma%f_LUNOPRT, *) "CARMA_Create: ERROR allocating groups, NGROUP=", &
carma%f_NGROUP, ", status=", ier
rc = RC_ERROR
return
endif
! Initialize
carma%f_icoag(:, :) = 0
carma%f_inucgas(:) = 0
! Allocate tables for the elements.
allocate( &
carma%f_element(NELEM), &
carma%f_igrowgas(NELEM), &
carma%f_inuc2elem(NELEM, NELEM), &
carma%f_inucproc(NELEM, NELEM), &
carma%f_ievp2elem(NELEM), &
carma%f_nnuc2elem(NELEM), &
carma%f_nnucelem(NELEM), &
carma%f_inucelem(NELEM,NELEM*NGROUP), &
carma%f_if_nuc(NELEM,NELEM), &
carma%f_rlh_nuc(NELEM, NELEM), &
carma%f_icoagelem(NELEM, NGROUP), &
carma%f_icoagelem_cm(NELEM, NGROUP), &
stat=ier)
if(ier /= 0) then
if (carma%f_do_print) write(carma%f_LUNOPRT, *) "CARMA_Create: ERROR allocating elements, NELEM=", &
carma%f_NELEM, ", status=", ier
rc = RC_ERROR
return
endif
! Initialize
carma%f_igrowgas(:) = 0
carma%f_inuc2elem(:,:) = 0
carma%f_inucproc(:,:) = 0
carma%f_ievp2elem(:) = 0
carma%f_nnuc2elem(:) = 0
carma%f_nnucelem(:) = 0
carma%f_inucelem(:,:) = 0
carma%f_if_nuc(:,:) = .FALSE.
carma%f_rlh_nuc(:,:) = 0._f
carma%f_icoagelem(:,:) = 0
carma%f_icoagelem_cm(:,:) = 0
! Allocate tables for the bins.
allocate( &
carma%f_inuc2bin(NBIN,NGROUP,NGROUP), &
carma%f_ievp2bin(NBIN,NGROUP,NGROUP), &
carma%f_nnucbin(NGROUP,NBIN,NGROUP), &
carma%f_inucbin(NBIN*NGROUP,NGROUP,NBIN,NGROUP), &
carma%f_diffmass(NBIN, NGROUP, NBIN, NGROUP), &
carma%f_volx(NGROUP,NGROUP,NGROUP,NBIN,NBIN), &
carma%f_ilow(NGROUP,NBIN,NBIN*NBIN), &
carma%f_jlow(NGROUP,NBIN,NBIN*NBIN), &
carma%f_iup(NGROUP,NBIN,NBIN*NBIN), &
carma%f_jup(NGROUP,NBIN,NBIN*NBIN), &
carma%f_npairl(NGROUP,NBIN), &
carma%f_npairu(NGROUP,NBIN), &
carma%f_iglow(NGROUP,NBIN,NBIN*NBIN), &
carma%f_jglow(NGROUP,NBIN,NBIN*NBIN), &
carma%f_igup(NGROUP,NBIN,NBIN*NBIN), &
carma%f_jgup(NGROUP,NBIN,NBIN*NBIN), &
carma%f_kbin(NGROUP,NGROUP,NGROUP,NBIN,NBIN), &
carma%f_pkernel(NBIN,NBIN,NGROUP,NGROUP,NGROUP,6), &
carma%f_pratt(3,NBIN,NGROUP), &
carma%f_prat(4,NBIN,NGROUP), &
carma%f_pden1(NBIN,NGROUP), &
carma%f_palr(4,NGROUP), &
stat=ier)
if(ier /= 0) then
if (carma%f_do_print) write(carma%f_LUNOPRT, *) "CARMA_Create: ERROR allocating bins, NBIN=", &
carma%f_NBIN, ", status=", ier
rc = RC_ERROR
return
endif
! Initialize
carma%f_inuc2bin(:,:,:) = 0
carma%f_ievp2bin(:,:,:) = 0
carma%f_nnucbin(:,:,:) = 0
carma%f_inucbin(:,:,:,:) = 0
carma%f_diffmass(:, :, :, :) = 0._f
carma%f_volx(:,:,:,:,:) = 0._f
carma%f_ilow(:,:,:) = 0
carma%f_jlow(:,:,:) = 0
carma%f_iup(:,:,:) = 0
carma%f_jup(:,:,:) = 0
carma%f_npairl(:,:) = 0
carma%f_npairu(:,:) = 0
carma%f_iglow(:,:,:) = 0
carma%f_jglow(:,:,:) = 0
carma%f_igup(:,:,:) = 0
carma%f_jgup(:,:,:) = 0
carma%f_kbin(:,:,:,:,:) = 0._f
carma%f_pkernel(:,:,:,:,:,:) = 0._f
carma%f_pratt(:,:,:) = 0._f
carma%f_prat(:,:,:) = 0._f
carma%f_pden1(:,:) = 0._f
carma%f_palr(:,:) = 0._f
! Allocate tables for solutes, if any are needed.
if (NSOLUTE > 0) then
allocate( &
carma%f_solute(NSOLUTE), &
stat=ier)
if(ier /= 0) then
if (carma%f_do_print) write(carma%f_LUNOPRT, *) "CARMA_Create: ERROR allocating solutes, NSOLUTE=", &
carma%f_NSOLUTE, ", status=", ier
rc = RC_ERROR
return
endif
end if
! Allocate tables for gases, if any are needed.
if (NGAS > 0) then
allocate( &
carma%f_gas(NGAS), &
stat=ier)
if(ier /= 0) then
if (carma%f_do_print) write(carma%f_LUNOPRT, *) "CARMA_Create: ERROR allocating gases, NGAS=", &
carma%f_NGAS, ", status=", ier
rc = RC_ERROR
return
endif
end if
! Allocate tables for optical properties, if any are needed.
if (NWAVE > 0) then
allocate( &
carma%f_wave(NWAVE), &
carma%f_dwave(NWAVE), &
carma%f_do_wave_emit(NWAVE), &
stat=ier)
if(ier /= 0) then
if (carma%f_do_print) write(carma%f_LUNOPRT, *) "CARMA_Create: ERROR allocating wavelengths, NWAVE=", &
carma%f_NWAVE, ", status=", ier
rc = RC_ERROR
return
endif
! Initialize
carma%f_do_wave_emit(:) = .TRUE.
if (present(wave)) carma%f_wave(:) = wave(:)
if (present(dwave)) carma%f_dwave(:) = dwave(:)
if (present(do_wave_emit)) carma%f_do_wave_emit(:) = do_wave_emit(:)
end if
return
end subroutine CARMA_Create
!! Called after the CARMA object has been created and the microphysics description has been
!! configured. The optional flags control which microphysical processes are enabled and all of
!! them default to FALSE. For a microphysical process to be active it must have been both
!! configured (using a CARMA_AddXXX() method) and enabled here.
!!
!! NOTE: After initialization, the structure of the particle size bins is determined, and
!! the resulting r, dr, rmass and dm can be retrieved with the CARMA_GetGroup() method.
!!
!! @version Feb-2009
!! @author Chuck Bardeen
subroutine CARMA_Initialize(carma, rc, do_cnst_rlh, do_coag, do_detrain, do_fixedinit, &
do_grow, do_incloud, do_explised, do_print_init, do_substep, do_thermo, do_vdiff, &
do_vtran, do_drydep, vf_const, minsubsteps, maxsubsteps, maxretries, conmax, &
do_pheat, do_pheatatm, dt_threshold, cstick, gsticki, gstickl, tstick, do_clearsky, &
do_partialinit, do_coremasscheck, sulfnucl_method )
type(carma_type), intent(inout) :: carma !! the carma object
integer, intent(out) :: rc !! return code, negative indicates failure
logical, intent(in), optional :: do_cnst_rlh !! use constant values for latent heats
!! (instead of varying with temperature)?
logical, intent(in), optional :: do_coag !! do coagulation?
logical, intent(in), optional :: do_detrain !! do detrainement?
logical, intent(in), optional :: do_fixedinit !! do initialization from reference atm?
logical, intent(in), optional :: do_grow !! do nucleation, growth and evaporation?
logical, intent(in), optional :: do_incloud !! do incloud growth and coagulation?
logical, intent(in), optional :: do_explised !! do sedimentation with substepping
logical, intent(in), optional :: do_substep !! do substepping
logical, intent(in), optional :: do_print_init !! do prinit initializtion information
logical, intent(in), optional :: do_thermo !! do thermodynamics
logical, intent(in), optional :: do_vdiff !! do Brownian diffusion
logical, intent(in), optional :: do_vtran !! do sedimentation
logical, intent(in), optional :: do_drydep !! do dry deposition
real(kind=f), intent(in), optional :: vf_const !! if specified and non-zero,
!! constant fall velocity for all particles [cm/s]
integer, intent(in), optional :: minsubsteps !! minimum number of substeps, default = 1
integer, intent(in), optional :: maxsubsteps !! maximum number of substeps, default = 1
integer, intent(in), optional :: maxretries !! maximum number of substep retries, default = 5
real(kind=f), intent(in), optional :: conmax !! minimum relative concentration to consider, default = 1e-1
logical, intent(in), optional :: do_pheat !! do particle heating
logical, intent(in), optional :: do_pheatatm !! do particle heating of atmosphere
real(kind=f), intent(in), optional :: dt_threshold !! convergence criteria for temperature [fraction]
real(kind=f), intent(in), optional :: cstick !! accommodation coefficient - coagulation, default = 1.0
real(kind=f), intent(in), optional :: gsticki !! accommodation coefficient - growth (ice), default = 0.93
real(kind=f), intent(in), optional :: gstickl !! accommodation coefficient - growth (liquid), default = 1.0
real(kind=f), intent(in), optional :: tstick !! accommodation coefficient - temperature, default = 1.0
logical, intent(in), optional :: do_clearsky !! do clear sky growth and coagulation?
logical, intent(in), optional :: do_partialinit !! do initialization of coagulation from reference atm (requires do_fixedinit)?
logical, intent(in), optional :: do_coremasscheck
character(len=*), intent(in), optional :: sulfnucl_method
! Assume success.
rc = RC_OK
! Set default values for control flags.
carma%f_do_cnst_rlh = .FALSE.
carma%f_do_coag = .FALSE.
carma%f_do_detrain = .FALSE.
carma%f_do_fixedinit = .FALSE.
carma%f_do_grow = .FALSE.
carma%f_do_incloud = .FALSE.
carma%f_do_explised = .FALSE.
carma%f_do_pheat = .FALSE.
carma%f_do_pheatatm = .FALSE.
carma%f_do_print_init = .FALSE.
carma%f_do_substep = .FALSE.
carma%f_do_thermo = .FALSE.
carma%f_do_vdiff = .FALSE.
carma%f_do_vtran = .FALSE.
carma%f_do_drydep = .FALSE.
carma%f_dt_threshold = 0._f
carma%f_cstick = 1._f
carma%f_gsticki = 0.93_f
carma%f_gstickl = 1._f
carma%f_tstick = 1._f
carma%f_do_clearsky = .FALSE.
carma%f_do_partialinit = .FALSE.
carma%f_do_coremasscheck = .FALSE.
! Store off any control flag values that have been supplied.
if (present(do_cnst_rlh)) carma%f_do_cnst_rlh = do_cnst_rlh
if (present(do_coag)) carma%f_do_coag = do_coag
if (present(do_detrain)) carma%f_do_detrain = do_detrain
if (present(do_fixedinit)) carma%f_do_fixedinit = do_fixedinit
if (present(do_grow)) carma%f_do_grow = do_grow
if (present(do_incloud)) carma%f_do_incloud = do_incloud
if (present(do_explised)) carma%f_do_explised = do_explised
if (present(do_pheat)) carma%f_do_pheat = do_pheat
if (present(do_pheatatm)) carma%f_do_pheatatm = do_pheatatm
if (present(do_print_init)) carma%f_do_print_init = (do_print_init .and. carma%f_do_print)
if (present(do_substep)) carma%f_do_substep = do_substep
if (present(do_thermo)) carma%f_do_thermo = do_thermo
if (present(do_vdiff)) carma%f_do_vdiff = do_vdiff
if (present(do_vtran)) carma%f_do_vtran = do_vtran
if (present(do_drydep)) carma%f_do_drydep = do_drydep
if (present(dt_threshold)) carma%f_dt_threshold = dt_threshold
if (present(cstick)) carma%f_cstick = cstick
if (present(gsticki)) carma%f_gsticki = gsticki
if (present(gstickl)) carma%f_gstickl = gstickl
if (present(tstick)) carma%f_tstick = tstick
if (present(do_clearsky)) carma%f_do_clearsky = do_clearsky
if (present(do_partialinit)) carma%f_do_partialinit = do_partialinit
if (present(do_coremasscheck)) carma%f_do_coremasscheck = do_coremasscheck
if (present(sulfnucl_method)) carma%sulfnucl_method = trim(sulfnucl_method)
! Setup the bin structure.
call setupbins(carma, rc)
if (rc < 0) return
! Substepping
carma%f_minsubsteps = 1 ! minimum number of substeps
carma%f_maxsubsteps = 1 ! maximum number of substeps
carma%f_maxretries = 1 ! maximum number of retries
carma%f_conmax = 1.e-1_f
if (present(minsubsteps)) carma%f_minsubsteps = minsubsteps
if (present(maxsubsteps)) carma%f_maxsubsteps = maxsubsteps
if (present(maxretries)) carma%f_maxretries = maxretries
if (present(conmax)) carma%f_conmax = conmax
carma%f_do_step = .TRUE.
! Calculate the Optical Properties
!
! NOTE: This is only needed by CARMA if particle heating is being used. For
! fractal particle the optics can be very slow, so only do it if necessary,
if (carma%f_do_pheat) then
call CARMA_InitializeOptics(carma, rc)
if (rc < 0) return
end if
! If any of the processes have initialization that can be done without the state
! information, then perform that now. This will mostly be checking the configuration
! and setting up any tables based upon the configuration.
if (carma%f_do_vtran .or. carma%f_do_coag) then
call CARMA_InitializeVertical(carma, rc, vf_const)
if (rc < 0) return
end if
if (carma%f_do_coag) then
call setupcoag(carma, rc)
if (rc < 0) return
end if
if (carma%f_do_grow) then
call CARMA_InitializeGrowth(carma, rc)
if (rc < 0) return
end if
if (carma%f_do_thermo) then
call CARMA_InitializeThermo(carma, rc)
if (rc < 0) return
end if
return
end subroutine CARMA_Initialize
subroutine CARMA_InitializeGrowth(carma, rc)
type(carma_type), intent(inout) :: carma
integer, intent(out) :: rc
! Local Variables
integer :: i
logical :: bad_grid
integer :: igroup ! group index
integer :: igas ! gas index
integer :: isol ! solute index
integer :: ielem ! element index
integer :: ibin ! bin index
integer :: igfrom
integer :: igto
integer :: ibto
integer :: ieto
integer :: ifrom
integer :: iefrom
integer :: jefrom
integer :: ip
integer :: jcore
integer :: iecore
integer :: im
integer :: jnucelem
integer :: inuc2
integer :: neto
integer :: jfrom
integer :: j
integer :: nnucb
! Define formats
1 format(a,': ',12i6)
2 format(/,a,': ',i6)
3 format(a,a)
4 format(a,': ',1pe12.3)
5 format(/,'Particle nucleation mapping arrays (setupnuc):')
7 format(/,'Warning: nucleation cannot occur from group',i3, &
' bin',i3,' into group',i3,' (<inuc2bin> is zero)')
! Assume success.
rc = RC_OK
! Compute radius-dependent terms used in PPM advection scheme
do igroup = 1, carma%f_NGROUP
do i = 2,carma%f_NBIN-1
carma%f_pratt(1,i,igroup) = carma%f_group(igroup)%f_dm(i) / &
( carma%f_group(igroup)%f_dm(i-1) + carma%f_group(igroup)%f_dm(i) + carma%f_group(igroup)%f_dm(i+1) )
carma%f_pratt(2,i,igroup) = ( 2._f*carma%f_group(igroup)%f_dm(i-1) + carma%f_group(igroup)%f_dm(i) ) / &
( carma%f_group(igroup)%f_dm(i+1) + carma%f_group(igroup)%f_dm(i) )
carma%f_pratt(3,i,igroup) = ( 2._f*carma%f_group(igroup)%f_dm(i+1) + carma%f_group(igroup)%f_dm(i) ) / &
( carma%f_group(igroup)%f_dm(i-1) + carma%f_group(igroup)%f_dm(i) )
enddo
do i = 2,carma%f_NBIN-2
carma%f_prat(1,i,igroup) = carma%f_group(igroup)%f_dm(i) / &
( carma%f_group(igroup)%f_dm(i) + carma%f_group(igroup)%f_dm(i+1) )
carma%f_prat(2,i,igroup) = 2._f * carma%f_group(igroup)%f_dm(i+1) * carma%f_group(igroup)%f_dm(i) / &
( carma%f_group(igroup)%f_dm(i) + carma%f_group(igroup)%f_dm(i+1) )
carma%f_prat(3,i,igroup) = ( carma%f_group(igroup)%f_dm(i-1) + carma%f_group(igroup)%f_dm(i) ) / &
( 2._f*carma%f_group(igroup)%f_dm(i) + carma%f_group(igroup)%f_dm(i+1) )
carma%f_prat(4,i,igroup) = ( carma%f_group(igroup)%f_dm(i+2) + carma%f_group(igroup)%f_dm(i+1) ) / &
( 2._f*carma%f_group(igroup)%f_dm(i+1) + carma%f_group(igroup)%f_dm(i) )
carma%f_pden1(i,igroup) = carma%f_group(igroup)%f_dm(i-1) + carma%f_group(igroup)%f_dm(i) + &
carma%f_group(igroup)%f_dm(i+1) + carma%f_group(igroup)%f_dm(i+2)
enddo
if( carma%f_NBIN .gt. 1 )then
carma%f_palr(1,igroup) = &
(carma%f_group(igroup)%f_rmassup(1)-carma%f_group(igroup)%f_rmass(1)) / &
(carma%f_group(igroup)%f_rmass(2)-carma%f_group(igroup)%f_rmass(1))
carma%f_palr(2,igroup) = &
(carma%f_group(igroup)%f_rmassup(1)/carma%f_group(igroup)%f_rmrat-carma%f_group(igroup)%f_rmass(1)) / &
(carma%f_group(igroup)%f_rmass(2)-carma%f_group(igroup)%f_rmass(1))
carma%f_palr(3,igroup) = &
(carma%f_group(igroup)%f_rmassup(carma%f_NBIN-1)-carma%f_group(igroup)%f_rmass(carma%f_NBIN-1)) &
/ (carma%f_group(igroup)%f_rmass(carma%f_NBIN)-carma%f_group(igroup)%f_rmass(carma%f_NBIN-1))
carma%f_palr(4,igroup) = &
(carma%f_group(igroup)%f_rmassup(carma%f_NBIN)-carma%f_group(igroup)%f_rmass(carma%f_NBIN-1)) &
/ (carma%f_group(igroup)%f_rmass(carma%f_NBIN)-carma%f_group(igroup)%f_rmass(carma%f_NBIN-1))
endif
end do
! Check the nucleation mapping.
!
! NOTE: This code was moved from setupnuc, because it is not dependent on the model's
! state. A small part of setupnuc which deals with scrit is state specific, and that was
! left in setupnuc.
! Bin mapping for nucleation : nucleation would transfer mass from particles
! in <ifrom,igfrom> into target bin <inuc2bin(ifrom,igfrom,igto)> in group
! <igto>. The target bin is the smallest bin in the target size grid with
! mass exceeding that of nucleated particle.
do igfrom = 1,carma%f_NGROUP ! nucleation source group
do igto = 1,carma%f_NGROUP ! nucleation target group
do ifrom = 1,carma%f_NBIN ! nucleation source bin
carma%f_inuc2bin(ifrom,igfrom,igto) = 0
do ibto = carma%f_NBIN,1,-1 ! nucleation target bin
if( carma%f_group(igto)%f_rmass(ibto) .ge. carma%f_group(igfrom)%f_rmass(ifrom) )then
carma%f_inuc2bin(ifrom,igfrom,igto) = ibto
endif
enddo
enddo
enddo
enddo
! Mappings for nucleation sources:
!
! <nnucelem(ielem)> is the number of particle elements that nucleate to
! particle element <ielem>.
!
! <inuc2elem(jefrom,ielem)> are the particle elements that
! nucleate to particle element <ielem>, where
! jefrom = 1,nnucelem(ielem).
!
! <if_nuc(iefrom,ieto)> is true if nucleation transfers mass from element
! <iefrom> to element <ieto>.
!
! <nnucbin(igfrom,ibin,igroup)> is the number of particle bins that nucleate
! to particles in bin <ibin,igroup> from group <igfrom>.
!
! <inucbin(jfrom,igfrom,ibin,igto)> are the particle bins
! that nucleate to particles in bin <ibin,igto>, where
! jfrom = 1,nnucbin(igfrom,ibin,igto).
!
!
! First, calculate <nnucelem(ielem)> and <if_nuc(iefrom,ieto)>
! based on <inucelem(jefrom,ielem)>
do iefrom = 1,carma%f_NELEM
do ieto = 1,carma%f_NELEM
carma%f_if_nuc(iefrom,ieto) = .false.
enddo
enddo
do ielem = 1,carma%f_NELEM
carma%f_nnuc2elem(ielem) = 0
do jefrom = 1,carma%f_NGROUP
if( carma%f_inuc2elem(jefrom,ielem) .ne. 0 ) then
carma%f_nnuc2elem(ielem) = carma%f_nnuc2elem(ielem) + 1
carma%f_if_nuc(ielem,carma%f_inuc2elem(jefrom,ielem)) = .true.
! Also check for cases where neither the source or destinaton don't have cores (e.g.
! melting ice to water drops).
if ((carma%f_group(carma%f_element(ielem)%f_igroup)%f_ncore .eq. 0) .and. &
(carma%f_group(carma%f_element(carma%f_inuc2elem(jefrom,ielem))%f_igroup)%f_ncore .eq. 0)) then
! For particle concentration target elements, only count source elements
! that are also particle concentrations.
carma%f_nnucelem(carma%f_inuc2elem(jefrom,ielem)) = carma%f_nnucelem(carma%f_inuc2elem(jefrom,ielem)) + 1
carma%f_inucelem(carma%f_nnucelem(carma%f_inuc2elem(jefrom,ielem)),carma%f_inuc2elem(jefrom,ielem)) = ielem
end if
endif
enddo
enddo
! Next, enumerate and count elements that nucleate to cores.
do igroup = 1,carma%f_NGROUP
ip = carma%f_group(igroup)%f_ienconc ! target particle number concentration element
do jcore = 1,carma%f_group(igroup)%f_ncore
iecore = carma%f_group(igroup)%f_icorelem(jcore) ! target core element
! carma%f_nnucelem(iecore) = 0
do iefrom = 1,carma%f_NELEM
if( carma%f_if_nuc(iefrom,iecore) ) then
carma%f_nnucelem(iecore) = carma%f_nnucelem(iecore) + 1
carma%f_inucelem(carma%f_nnucelem(iecore),iecore) = iefrom
endif
enddo ! iefrom=1,NELEM
enddo ! jcore=1,ncore
enddo ! igroup=1,NGROUP
! Now enumerate and count elements nucleating to particle concentration
! (itype=I_INVOLATILE and itype=I_VOLATILE) and core second moment
! (itype=I_COREMASS). Elements with itype = I_VOLATILE are special because all
! nucleation sources for core elements in same group are also sources
! for the itype = I_VOLATILE element.
do igroup = 1,carma%f_NGROUP
ip = carma%f_group(igroup)%f_ienconc ! target particle number concentration element
im = carma%f_group(igroup)%f_imomelem ! target core second moment element
! carma%f_nnucelem(ip) = 0
! if( im .ne. 0 )then
! carma%f_nnucelem(im) = 0
! endif
do jcore = 1,carma%f_group(igroup)%f_ncore
iecore = carma%f_group(igroup)%f_icorelem(jcore) ! target core mass element
do jnucelem = 1,carma%f_nnucelem(iecore) ! elements nucleating to cores
iefrom = carma%f_inucelem(jnucelem,iecore) ! source
! For particle concentration target elements, only count source elements
! that are also particle concentrations.
carma%f_nnucelem(ip) = carma%f_nnucelem(ip) + 1
carma%f_inucelem(carma%f_nnucelem(ip),ip) = carma%f_group(carma%f_element(iefrom)%f_igroup)%f_ienconc
if( im .ne. 0 )then
carma%f_nnucelem(im) = carma%f_nnucelem(im) + 1
carma%f_inucelem(carma%f_nnucelem(im),im) = iefrom
endif
enddo
enddo ! jcore=1,ncore
enddo ! igroup=1,NGROUP
! Now enumerate and count nucleating bins.
do igroup = 1,carma%f_NGROUP ! target group
do ibin = 1,carma%f_NBIN ! target bin
do igfrom = 1,carma%f_NGROUP ! source group
carma%f_nnucbin(igfrom,ibin,igroup) = 0
do ifrom = 1,carma%f_NBIN ! source bin
if( carma%f_inuc2bin(ifrom,igfrom,igroup) .eq. ibin ) then
carma%f_nnucbin(igfrom,ibin,igroup) = carma%f_nnucbin(igfrom,ibin,igroup) + 1
carma%f_inucbin(carma%f_nnucbin(igfrom,ibin,igroup),igfrom,ibin,igroup) = ifrom
endif
enddo
enddo ! igfrom=1,NGROUP
enddo ! ibin=1,NBIN=1,NGROUP
enddo ! igroup=1,NGROUP
if (carma%f_do_print_init) then
! Report nucleation mapping arrays (should be 'write' stmts, of course)
write(carma%f_LUNOPRT,*) ' '
write(carma%f_LUNOPRT,*) 'Nucleation mapping arrays (setupnuc):'
write(carma%f_LUNOPRT,*) ' '
write(carma%f_LUNOPRT,*) 'Elements mapping:'
do ielem = 1,carma%f_NELEM
write(carma%f_LUNOPRT,*) 'ielem,nnucelem=',ielem,carma%f_nnucelem(ielem)
if(carma%f_nnucelem(ielem) .gt. 0) then
do jfrom = 1,carma%f_nnucelem(ielem)
write(carma%f_LUNOPRT,*) 'jfrom,inucelem= ',jfrom,carma%f_inucelem(jfrom,ielem)
enddo
endif
enddo
write(carma%f_LUNOPRT,*) ' '
write(carma%f_LUNOPRT,*) 'Bin mapping:'
do igfrom = 1,carma%f_NGROUP
do igroup = 1,carma%f_NGROUP
write(carma%f_LUNOPRT,*) ' '
write(carma%f_LUNOPRT,*) 'Groups (from, to) = ', igfrom, igroup
do ibin = 1,carma%f_NBIN
nnucb = carma%f_nnucbin(igfrom,ibin,igroup)
if(nnucb .eq. 0) write(carma%f_LUNOPRT,*) ' None for bin ',ibin
if(nnucb .gt. 0) then
write(carma%f_LUNOPRT,*) ' ibin,nnucbin=',ibin,nnucb
write(carma%f_LUNOPRT,*) ' inucbin=',(carma%f_inucbin(j,igfrom,ibin,igroup),j=1,nnucb)
endif
enddo
enddo
enddo
endif
! Check that values are valid.
do ielem = 1, carma%f_NELEM
if( carma%f_element(ielem)%f_isolute .gt. carma%f_NSOLUTE )then
if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::ERROR - component of isolute > NSOLUTE'
rc = RC_ERROR
return
endif
if( carma%f_ievp2elem(ielem) .gt. carma%f_NELEM )then
if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::ERROR - component of ievp2elem > NELEM'
rc = RC_ERROR
return
endif
! Check that <isolute> is consistent with <ievp2elem>.
if( carma%f_ievp2elem(ielem) .ne. 0 .and. carma%f_element(ielem)%f_itype .eq. I_COREMASS )then
if( carma%f_element(ielem)%f_isolute .ne. carma%f_element(carma%f_ievp2elem(ielem))%f_isolute)then
if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::ERROR - isolute and ievp2elem are inconsistent'
rc = RC_ERROR
return
endif
endif
! Check that <isolute> is consistent with <inucgas>.
! igas = carma%f_inucgas( carma%f_element(ielem)%f_igroup )
! if( igas .ne. 0 )then
! if( carma%f_element(ielem)%f_itype .eq. I_COREMASS .and. carma%f_element(ielem)%f_isolute .eq. 0 )then
! if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::ERROR - inucgas ne 0 but isolute eq 0'
! rc = RC_ERROR
! return
! endif
! endif
enddo
do ielem = 1, carma%f_NELEM
if( carma%f_nnuc2elem(ielem) .gt. 0 ) then
do inuc2 = 1, carma%f_nnuc2elem(ielem)
if( carma%f_inuc2elem(inuc2,ielem) .gt. carma%f_NELEM )then
if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::ERROR - component of inuc2elem > NELEM'
rc = RC_ERROR
return
endif
enddo
endif
enddo
! Particle grids are incompatible if there is no target bin with enough
! mass to accomodate nucleated particle.
bad_grid = .false.
do iefrom = 1,carma%f_NELEM ! source element
igfrom = carma%f_element(iefrom)%f_igroup
neto = carma%f_nnuc2elem(iefrom)
if( neto .gt. 0 )then
do inuc2 = 1,neto
ieto = carma%f_inuc2elem(inuc2,iefrom)
igto = carma%f_element(ieto)%f_igroup
do ifrom = 1,carma%f_NBIN ! source bin
if( carma%f_inuc2bin(ifrom,igfrom,igto) .eq. 0 )then
if ((carma%f_do_print) .and. (carma%f_do_print_init)) write(carma%f_LUNOPRT,7) igfrom,ifrom,igto
bad_grid = .true.
endif
enddo
enddo
endif
enddo
if (carma%f_do_print_init) then
if( bad_grid )then
if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::Warning - incompatible grids for nucleation'
endif
! Report some initialization values!
write(carma%f_LUNOPRT,5)
write(carma%f_LUNOPRT,1) 'inucgas ',(carma%f_inucgas(i),i=1,carma%f_NGROUP)
write(carma%f_LUNOPRT,1) 'inuc2elem',(carma%f_inuc2elem(1,i),i=1,carma%f_NELEM)
write(carma%f_LUNOPRT,1) 'ievp2elem',(carma%f_ievp2elem(i),i=1,carma%f_NELEM)
write(carma%f_LUNOPRT,1) 'isolute ',(carma%f_element(i)%f_isolute,i=1,carma%f_NELEM)
do isol = 1,carma%f_NSOLUTE
write(carma%f_LUNOPRT,2) 'solute number ',isol
write(carma%f_LUNOPRT,3) 'solute name: ',carma%f_solute(isol)%f_name
write(carma%f_LUNOPRT,4) 'molecular weight',carma%f_solute(isol)%f_wtmol
write(carma%f_LUNOPRT,4) 'mass density ',carma%f_solute(isol)%f_rho
enddo
endif
! Initialize indexes for the gases and check to make sure if H2SO4 is used
! that it occurs after H2O. This is necessary for supersaturation calculations.
carma%f_igash2o = 0
carma%f_igash2so4 = 0
carma%f_igasso2 = 0
do igas = 1, carma%f_NGAS
if (carma%f_gas(igas)%f_icomposition == I_GCOMP_H2O) then
carma%f_igash2o = igas
else if (carma%f_gas(igas)%f_icomposition == I_GCOMP_H2SO4) then
carma%f_igash2so4 = igas
else if (carma%f_gas(igas)%f_icomposition == I_GCOMP_SO2) then
carma%f_igasso2 = igas
end if
end do
if ((carma%f_igash2so4 /= 0) .and. (carma%f_igash2o > carma%f_igash2so4)) then
if (carma%f_do_print) write(carma%f_LUNOPRT,*) 'CARMA_InitializeGrowth::ERROR - H2O gas must come before H2SO4.'
rc = RC_ERROR
return
end if
return
end subroutine CARMA_InitializeGrowth
!! Calculate the optical properties for each particle bin at each of
!! the specified wavelengths. The optical properties include the
!! extinction efficiency, the single scattering albedo and the
!! asymmetry factor.
!!
!! NOTE: For these calculations, the particles are assumed to be spheres and
!! Mie code is used to calculate the optical properties.
!!
!! @author Chuck Bardeen
!! @version May-2009
subroutine CARMA_InitializeOptics(carma, rc)
type(carma_type), intent(inout) :: carma
integer, intent(out) :: rc
integer :: igroup ! group index
integer :: iwave ! wavelength index
integer :: ibin ! bin index
real(kind=f) :: Qext
real(kind=f) :: Qsca
real(kind=f) :: asym
! Assume success.
rc = RC_OK
! Were any wavelengths specified?
do iwave = 1, carma%f_NWAVE
do igroup = 1, carma%f_NGROUP
! Should we calculate mie properties for this group?
if (carma%f_group(igroup)%f_do_mie) then
do ibin = 1, carma%f_NBIN
! Assume the particle is homogeneous (no core).
!
! Refractive index comes from the concentration element.
!
! NOTE: The miess does not converge over as broad a
! range of input parameters as bhmie, but it can handle
! coated spheres.
call mie(carma, &
carma%f_group(igroup)%f_imiertn, &
carma%f_group(igroup)%f_r(ibin), &
carma%f_wave(iwave), &
carma%f_group(igroup)%f_nmon(ibin), &
carma%f_group(igroup)%f_df(ibin), &
carma%f_group(igroup)%f_rmon, &
carma%f_group(igroup)%f_falpha, &
carma%f_element(carma%f_group(igroup)%f_ienconc)%f_refidx(iwave, 1), &
0.0_f, &
0.0_f, &
Qext, &
Qsca, &
asym, &
rc)
if (rc < RC_OK) then
if (carma%f_do_print) then
write(carma%f_LUNOPRT, *) "CARMA_InitializeOptics::&
&Mie failed for (band, wavelength, group, bin)", &
iwave, carma%f_wave(iwave), igroup, ibin
end if
return
end if
carma%f_group(igroup)%f_qext(iwave, ibin) = Qext
carma%f_group(igroup)%f_ssa(iwave, ibin) = Qsca / Qext
carma%f_group(igroup)%f_asym(iwave, ibin) = asym
end do
end if
end do
end do
return
end subroutine CARMA_InitializeOptics
!! Perform initialization of variables related to thermodynamical calculations that
!! are not dependent on the model state.
!!
!! @author Chuck Bardeen
!! @version May-2009
subroutine CARMA_InitializeThermo(carma, rc)
type(carma_type), intent(inout) :: carma
integer, intent(out) :: rc
! Assume success.
rc = RC_OK
return
end subroutine CARMA_InitializeThermo
!! Perform initialization of variables related to vertical transport that are not dependent
!! on the model state.
!!
!! @author Chuck Bardeen
!! @version May-2009
subroutine CARMA_InitializeVertical(carma, rc, vf_const)
type(carma_type), intent(inout) :: carma
integer, intent(out) :: rc
real(kind=f), intent(in), optional :: vf_const
! Assume success.
rc = RC_OK
! Was a constant vertical velocity specified?
carma%f_ifall = 1
carma%f_vf_const = 0._f
if (present(vf_const)) then
if (vf_const /= 0._f) then
carma%f_ifall = 0
carma%f_vf_const = vf_const
end if
end if