80 Agenda& iy_cloudbox_agenda,
84 cloudbox_limits.resize ( 0 );
85 iy_cloudbox_agenda =
Agenda();
86 iy_cloudbox_agenda.
set_name (
"iy_cloudbox_agenda" );
98 const Index& atmosphere_dim,
103 const Tensor4& massdensity_field,
105 const Numeric& cloudbox_margin,
119 Index type_flag=0, i=0, j=0, k=0, l=0;
128 chk_atm_grids ( atmosphere_dim, p_grid, lat_grid, lon_grid );
133 cloudbox_limits.resize ( atmosphere_dim*2 );
136 for ( l=0; l<part_species.
nelem(); l++ )
146 if ( part_type ==
"LWC" ) type_flag = 0;
147 else if ( part_type ==
"IWC" ) type_flag = 1;
148 else if ( part_type ==
"Rain" ) type_flag = 2;
149 else if ( part_type ==
"Snow" ) type_flag = 3;
166 if ( atmosphere_dim == 1 )
173 if ( cloudbox_margin == -1 )
175 cloudbox_limits[0] = 0;
182 for ( i=0; i<hydro_p.
nelem(); i++ )
184 if ( hydro_p[i] != 0.0 )
199 for ( j=hydro_p.
nelem()-1; j>=i; j-- )
201 if ( hydro_p[j] != 0.0 )
289 while ( p_grid[k+1] >= p_margin1 && k+1 < p_grid.
nelem() ) k++;
290 cloudbox_limits[0]= k;
294 p2 =
min(p2+1, massdensity_field.
npages()-1);
297 if ( p2 >= massdensity_field.
npages()-1)
300 out2<<
"The cloud reaches to TOA!\n"
301 <<
"Check massdensity_field data, if realistic!\n";
303 cloudbox_limits[1] = p2;
314 out0<<
"Cloudbox is switched off!\n";
323 assert ( p_grid[p1] > p_grid[p2] );
325 assert ( p_grid[p1] > p_grid[p_grid.
nelem()-1] );
327 assert ( p_grid[p2] < p_grid[0] );
329 if ( atmosphere_dim >= 2 )
332 assert ( lat_grid[lat2] > lat_grid[lat1] );
334 assert ( lat_grid[lat1] >= lat_grid[1] );
336 assert ( lat_grid[lat2] <= lat_grid[lat_grid.
nelem()-2] );
338 if ( atmosphere_dim == 3 )
341 assert ( lon_grid[lon2] > lon_grid[lon1] );
343 assert ( lon_grid[lon1] >= lon_grid[1] );
345 assert ( lon_grid[lon2] <= lon_grid[lon_grid.
nelem()-2] );
354 const Index& atmosphere_dim,
373 throw runtime_error(
"The pressure in *p1* must be bigger than the "
374 "pressure in *p2*." );
375 if( p1 <= p_grid[p_grid.
nelem()-1] )
376 throw runtime_error(
"The pressure in *p1* must be larger than the "
377 "last value in *p_grid*." );
378 if( p2 >= p_grid[0] )
379 throw runtime_error(
"The pressure in *p2* must be smaller than the "
380 "first value in *p_grid*." );
381 if( atmosphere_dim >= 2 )
384 throw runtime_error(
"The latitude in *lat2* must be bigger than the "
385 "latitude in *lat1*.");
386 if( lat1 < lat_grid[1] )
387 throw runtime_error(
"The latitude in *lat1* must be >= the "
388 "second value in *lat_grid*." );
389 if( lat2 > lat_grid[lat_grid.
nelem()-2] )
390 throw runtime_error(
"The latitude in *lat2* must be <= the "
391 "next to last value in *lat_grid*." );
393 if( atmosphere_dim == 3 )
396 throw runtime_error(
"The longitude in *lon2* must be bigger than the "
397 "longitude in *lon1*.");
398 if( lon1 < lon_grid[1] )
399 throw runtime_error(
"The longitude in *lon1* must be >= the "
400 "second value in *lon_grid*." );
401 if( lon2 > lon_grid[lon_grid.
nelem()-2] )
402 throw runtime_error(
"The longitude in *lon2* must be <= the "
403 "next to last value in *lon_grid*." );
410 cloudbox_limits.resize( atmosphere_dim*2 );
415 cloudbox_limits[0] = 0;
419 for( cloudbox_limits[0]=1; p_grid[cloudbox_limits[0]+1]>=p1;
420 cloudbox_limits[0]++ ) {}
422 if( p2 < p_grid[p_grid.
nelem()-2] )
424 cloudbox_limits[1] = p_grid.
nelem() - 1;
428 for( cloudbox_limits[1]=p_grid.
nelem()-2;
429 p_grid[cloudbox_limits[1]-1]<=p2; cloudbox_limits[1]-- ) {}
433 if( atmosphere_dim >= 2 )
435 for( cloudbox_limits[2]=1; lat_grid[cloudbox_limits[2]+1]<=lat1;
436 cloudbox_limits[2]++ ) {}
437 for( cloudbox_limits[3]=lat_grid.
nelem()-2;
438 lat_grid[cloudbox_limits[3]-1]>=lat2; cloudbox_limits[3]-- ) {}
442 if( atmosphere_dim == 3 )
444 for( cloudbox_limits[4]=1; lon_grid[cloudbox_limits[4]+1]<=lon1;
445 cloudbox_limits[4]++ ) {}
446 for( cloudbox_limits[5]=lon_grid.
nelem()-2;
447 lon_grid[cloudbox_limits[5]-1]>=lon2; cloudbox_limits[5]-- ) {}
457 const Index& atmosphere_dim,
475 throw runtime_error(
"The altitude in *z1* must be smaller than the "
476 "altitude in *z2*." );
485 if( atmosphere_dim == 3 )
488 throw runtime_error(
"The latitude in *lat2* must be bigger than the "
489 " latitude in *lat1*.");
490 if( lat1 < lat_grid[1] )
491 throw runtime_error(
"The latitude in *lat1* must be >= the "
492 "second value in *lat_grid*." );
493 if( lat2 > lat_grid[lat_grid.
nelem()-2] )
494 throw runtime_error(
"The latitude in *lat2* must be <= the "
495 "next to last value in *lat_grid*." );
497 throw runtime_error(
"The longitude in *lon2* must be bigger than the "
498 "longitude in *lon1*.");
499 if( lon1 < lon_grid[1] )
500 throw runtime_error(
"The longitude in *lon1* must be >= the "
501 "second value in *lon_grid*." );
502 if( lon2 > lon_grid[lon_grid.
nelem()-2] )
503 throw runtime_error(
"The longitude in *lon2* must be <= the "
504 "next to last value in *lon_grid*." );
511 cloudbox_limits.resize( atmosphere_dim*2 );
514 if( z1 < z_field(1, 0, 0) )
516 cloudbox_limits[0] = 0;
520 for( cloudbox_limits[0]=1; z_field(cloudbox_limits[0]+1, 0, 0) <= z1;
521 cloudbox_limits[0]++ ) {}
523 if( z2 > z_field(z_field.
npages()-2, 0, 0) )
525 cloudbox_limits[1] = z_field.
npages() - 1;
529 for( cloudbox_limits[1]=z_field.
npages()- 2;
530 z_field(cloudbox_limits[1]-1, 0, 0) >= z2; cloudbox_limits[1]-- ) {}
534 if( atmosphere_dim >= 2 )
536 for( cloudbox_limits[2]=1; lat_grid[cloudbox_limits[2]+1]<=lat1;
537 cloudbox_limits[2]++ ) {}
538 for( cloudbox_limits[3]=lat_grid.
nelem()-2;
539 lat_grid[cloudbox_limits[3]-1]>=lat2; cloudbox_limits[3]-- ) {}
543 if( atmosphere_dim == 3 )
545 for( cloudbox_limits[4]=1; lon_grid[cloudbox_limits[4]+1]<=lon1;
546 cloudbox_limits[4]++ ) {}
547 for( cloudbox_limits[5]=lon_grid.
nelem()-2;
548 lon_grid[cloudbox_limits[5]-1]>=lon2; cloudbox_limits[5]-- ) {}
556 const Index& basics_checked,
557 const Index& atmosphere_dim,
564 const Index& cloudbox_on,
571 if( !basics_checked )
572 throw runtime_error(
"The atmosphere and basic control varaibles must be "
573 "flagged to have passed a consistency check (basics_checked=1)." );
582 ow <<
"The scattering methods are not (yet?) handling winds. For this\n"
583 <<
"reason, the WSVs for wind fields must all be empty with an\n."
584 <<
"active cloudbox.";
585 if( wind_w_field.
npages() > 0 )
586 {
throw runtime_error( ow.str() ); }
587 if( wind_v_field.
npages() > 0 )
588 {
throw runtime_error( ow.str() ); }
589 if( atmosphere_dim > 2 && wind_u_field.
npages() > 0 )
590 {
throw runtime_error( ow.str() ); }
594 if( cloudbox_limits.
nelem() != atmosphere_dim*2 )
597 os <<
"The array *cloudbox_limits* has incorrect length.\n"
598 <<
"For atmospheric dim. = " << atmosphere_dim
599 <<
" the length shall be " << atmosphere_dim*2
600 <<
" but it is " << cloudbox_limits.
nelem() <<
".";
601 throw runtime_error( os.str() );
603 if( cloudbox_limits[1]<=cloudbox_limits[0] || cloudbox_limits[0]<0 ||
604 cloudbox_limits[1]>=p_grid.
nelem() )
607 os <<
"Incorrect value(s) for cloud box pressure limit(s) found."
608 <<
"\nValues are either out of range or upper limit is not "
609 <<
"greater than lower limit.\nWith present length of "
610 <<
"*p_grid*, OK values are 0 - " << p_grid.
nelem()-1
611 <<
".\nThe pressure index limits are set to "
612 << cloudbox_limits[0] <<
" - " << cloudbox_limits[1] <<
".";
613 throw runtime_error( os.str() );
615 if( atmosphere_dim >= 2 )
618 if( cloudbox_limits[3]<=cloudbox_limits[2] || cloudbox_limits[2]<1 ||
619 cloudbox_limits[3]>=n-1 )
622 os <<
"Incorrect value(s) for cloud box latitude limit(s) found."
623 <<
"\nValues are either out of range or upper limit is not "
624 <<
"greater than lower limit.\nWith present length of "
625 <<
"*lat_grid*, OK values are 1 - " << n-2
626 <<
".\nThe latitude index limits are set to "
627 << cloudbox_limits[2] <<
" - " << cloudbox_limits[3] <<
".";
628 throw runtime_error( os.str() );
630 if( ( lat_grid[cloudbox_limits[2]] - lat_grid[0] < llmin ) &&
631 ( atmosphere_dim==2 ||
632 ( atmosphere_dim==3 && lat_grid[0]>-90) ) )
635 os <<
"Too small distance between cloudbox and lower end of\n"
636 <<
"latitude grid. This distance must be " << llmin
637 <<
" degrees. Cloudbox ends at " << lat_grid[cloudbox_limits[2]]
638 <<
" and latitude grid starts at " << lat_grid[0] <<
".";
639 throw runtime_error( os.str() );
641 if( ( lat_grid[n-1] - lat_grid[cloudbox_limits[3]] < llmin ) &&
642 ( atmosphere_dim==2 ||
643 (atmosphere_dim==3 && lat_grid[n-1]<90) ) )
646 os <<
"Too small distance between cloudbox and upper end of\n"
647 <<
"latitude grid. This distance must be " << llmin
648 <<
"degrees. Cloudbox ends at " << lat_grid[cloudbox_limits[3]]
649 <<
" and latitude grid ends at " << lat_grid[n-1] <<
".";
650 throw runtime_error( os.str() );
653 if( atmosphere_dim >= 3 )
656 if( cloudbox_limits[5]<=cloudbox_limits[4] || cloudbox_limits[4]<1 ||
657 cloudbox_limits[5]>=n-1 )
660 os <<
"Incorrect value(s) for cloud box longitude limit(s) found"
661 <<
".\nValues are either out of range or upper limit is not "
662 <<
"greater than lower limit.\nWith present length of "
663 <<
"*lon_grid*, OK values are 1 - " << n-2
664 <<
".\nThe longitude limits are set to "
665 << cloudbox_limits[4] <<
" - " << cloudbox_limits[5] <<
".";
666 throw runtime_error( os.str() );
668 if( lon_grid[n-1] - lon_grid[0] < 360 )
670 const Numeric latmax =
max(
abs(lat_grid[cloudbox_limits[2]]),
671 abs(lat_grid[cloudbox_limits[3]]) );
673 if( lon_grid[cloudbox_limits[4]]-lon_grid[0] < llmin/lfac )
676 os <<
"Too small distance between cloudbox and lower end of\n"
677 <<
"longitude grid. This distance must here be "
678 << llmin/lfac <<
" degrees.";
679 throw runtime_error( os.str() );
681 if( lon_grid[n-1] - lon_grid[cloudbox_limits[5]] < llmin/lfac )
684 os <<
"Too small distance between cloudbox and upper end of\n"
685 <<
"longitude grid. This distance must here be "
686 << llmin/lfac <<
" degrees.";
687 throw runtime_error( os.str() );
694 cloudbox_checked = 1;
703 const Numeric& massdensity_threshold,
708 for (
Index i=0; i<massdensity_field.
nbooks(); i++ )
710 for (
Index j=0; j<massdensity_field.
npages(); j++ )
712 for (
Index k=0; k<massdensity_field.
nrows(); k++ )
714 for (
Index l=0; l<massdensity_field.
ncols(); l++ )
716 if ( massdensity_field ( i,j,k,l ) < massdensity_threshold )
718 massdensity_field ( i,j,k,l ) = 0.0;
731 part_species.resize ( 0 );
744 part_species.resize ( names.
nelem() );
746 part_species = names;
749 out3 <<
" Defined particle settings: ";
750 for (
Index i=0; i<part_species.
nelem(); ++i )
752 out3 <<
"\n " << i <<
": "<<part_species[i];
764 scat_data_raw.reserve ( 20 );
765 pnd_field_raw.reserve ( 20 );
774 const Index& atmosphere_dim,
781 const String& filename_scat_data,
782 const String& pnd_field_file,
791 chk_atm_grids ( atmosphere_dim, p_grid, lat_grid, lon_grid );
794 if ( cloudbox_limits.
nelem() != 2*atmosphere_dim )
795 throw runtime_error (
796 "*cloudbox_limits* is a vector which contains"
797 "the upper and lower\n"
798 "limit of the cloud for all atmospheric dimensions.\n"
799 "So its length must be 2 x *atmosphere_dim*" );
801 if ( f_grid.
nelem() == 0 )
802 throw runtime_error (
"The frequency grid is empty." );
809 scat_data_raw.resize ( data_files.
nelem() );
811 for (
Index i = 0; i<data_files.
nelem(); i++ )
814 out2 <<
" Read single scattering data\n";
818 data_files[i], f_grid,
823 out2 <<
" Read particle number density data \n";
827 pnd_field_file, atmosphere_dim, p_grid, lat_grid,
828 lon_grid, cloudbox_limits, verbosity);
837 const String& filename_scat_data,
838 const String& filename_scat_meta_data,
849 scat_data_raw.resize ( data_files.nelem() );
851 for (
Index i = 0; i<data_files.nelem(); i++ )
853 out3 <<
" Read single scattering data\n";
857 data_files[i], f_grid,
864 scat_data_meta_array.resize ( meta_data_files.
nelem() );
866 for (
Index i = 0; i<meta_data_files.
nelem(); i++ )
869 out3 <<
" Read scattering meta data\n";
873 meta_data_files[i], verbosity );
879 scat_data_meta_array, verbosity );
898 Index intarr_total = 0;
905 scat_data_nelem.resize( part_species.
nelem() );
908 selected.resize(scat_data_meta_array_tmp.
nelem());
911 for (
Index k=0; k<part_species.
nelem(); k++ )
921 if ( part_type ==
"IWC" || part_type==
"Snow" ) type =
"Ice";
922 else if ( part_type==
"LWC" || part_type ==
"Rain" ) type =
"Water";
929 for (
Index j=0; j<scat_data_meta_array_tmp.
nelem(); j++ )
932 if ( scat_data_meta_array_tmp[j].type == type )
937 pow ( 3./4. * scat_data_meta_array_tmp[j].V * 1e18 /
PI , 1./3. );
942 if ( r_particle >= sizemin &&
943 ( sizemax >= r_particle || sizemax < 0. ))
946 intarr.push_back ( j );
949 out3 <<
"Selecting particle " << j+1 <<
"/" << scat_data_meta_array_tmp.
nelem()
950 <<
" (" << scat_data_meta_array_tmp[j].type <<
")\n";
955 scat_data_nelem[k] = intarr.
nelem() - intarr_total;
956 intarr_total = intarr.
nelem();
959 if ( !intarr.
nelem() )
962 os <<
"The selection in " << part_species <<
963 " is NOT choosing any of the given Scattering Data.\n"
964 "--> Does the selection in *part_species* fit any of the "
965 "Single Scattering Data input? \n";
966 throw runtime_error ( os.str() );
973 out1 <<
"WARNING! Ignored SMD[" << j <<
"] (" << scat_data_meta_array_tmp[j].type <<
")!\n";
979 scat_data_raw.resize ( intarr.
nelem() );
980 scat_data_meta_array.resize ( intarr.
nelem() );
985 scat_data_meta_array[j] = scat_data_meta_array_tmp[intarr[j]] ;
986 scat_data_raw[j] = scat_data_raw_tmp[intarr[j]] ;
998 const Index& atmosphere_dim,
1005 const String& scat_data_file,
1006 const String& pnd_field_file,
1015 chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1018 if ( cloudbox_limits.
nelem() != 2*atmosphere_dim )
1019 throw runtime_error(
1020 "*cloudbox_limits* is a vector which contains"
1021 "the upper and lower\n"
1022 "limit of the cloud for all atmospheric dimensions.\n"
1023 "So its length must be 2 x *atmosphere_dim*" );
1025 if( f_grid.
nelem() == 0 )
1026 throw runtime_error(
"The frequency grid is empty." );
1034 scat_data_raw.push_back(scat_data);
1037 pnd_field_raw.push_back(pnd_field_data);
1039 out2 <<
" Read single scattering data\n";
1044 scat_data_file, f_grid, verbosity);
1046 out2 <<
" Read particle number density field\n";
1047 if (pnd_field_file.
nelem() < 1)
1050 out1 <<
"Warning: No pnd_field_file specified. Ignored. \n";
1058 pnd_field_file, atmosphere_dim, p_grid, lat_grid,
1059 lon_grid, cloudbox_limits, verbosity);
1072 const Index& atmosphere_dim,
1080 if (pnd_field_raw.
nelem() == 0)
1081 throw runtime_error(
1082 "No particle number density data given. Please\n"
1083 "use WSMs *ParticleTypeInit* and \n"
1084 "*ParticleTypeAdd(All)* for reading cloud particle\n"
1088 chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1089 if ( cloudbox_limits.
nelem()!= 2*atmosphere_dim)
1090 throw runtime_error(
1091 "*cloudbox_limits* is a vector which contains the"
1092 "upper and lower limit of the cloud for all "
1093 "atmospheric dimensions. So its dimension must"
1094 "be 2 x *atmosphere_dim*");
1099 for (
Index d = 0; d < atmosphere_dim; d++)
1101 for (
Index i = 0; i < pnd_field_raw.
nelem(); i++)
1103 if (pnd_field_raw[i].get_grid_size(d) < 2)
1106 os <<
"Error in pnd_field_raw data. ";
1107 os <<
"Dimension " << d <<
" (name: \"";
1108 os << pnd_field_raw[i].get_grid_name(d);
1109 os <<
"\") has only ";
1110 os << pnd_field_raw[i].get_grid_size(d);
1112 os << ((pnd_field_raw[i].get_grid_size(d)==1) ?
"" :
"s");
1113 os <<
". Must be at least 2.";
1114 throw runtime_error(os.str());
1118 const Index Np_cloud = cloudbox_limits[1]-cloudbox_limits[0]+1;
1124 p_grid, lat_grid, lon_grid,
1128 if ( atmosphere_dim == 1)
1131 pnd_field.
resize(pnd_field_raw.
nelem(), Np_cloud, 1, 1 );
1138 for (
Index i = 0; i < pnd_field_raw.
nelem(); ++ i)
1150 pnd_field_raw[i].data(
joker,0,0), gp_p );
1154 else if(atmosphere_dim == 2)
1156 const Index Nlat_cloud = cloudbox_limits[3]-cloudbox_limits[2]+1;
1159 lat_grid[
Range(cloudbox_limits[2],Nlat_cloud)];
1162 pnd_field.
resize( pnd_field_raw.
nelem(), Np_cloud, Nlat_cloud, 1 );
1170 for (
Index i = 0; i < pnd_field_raw.
nelem(); ++ i)
1179 Tensor3 itw( Np_cloud, Nlat_cloud, 4 );
1184 pnd_field_raw[i].data(
joker,
joker,0), gp_p, gp_lat );
1189 const Index Nlat_cloud = cloudbox_limits[3]-cloudbox_limits[2]+1;
1190 const Index Nlon_cloud = cloudbox_limits[5]-cloudbox_limits[4]+1;
1193 lat_grid[
Range(cloudbox_limits[2],Nlat_cloud)];
1195 lon_grid[
Range(cloudbox_limits[4],Nlon_cloud)];
1198 pnd_field.
resize( pnd_field_raw.
nelem(), Np_cloud, Nlat_cloud,
1208 for (
Index i = 0; i < pnd_field_raw.
nelem(); ++ i)
1219 Tensor4 itw( Np_cloud, Nlat_cloud, Nlon_cloud, 8 );
1224 pnd_field_raw[i].data, gp_p, gp_lat, gp_lon );
1232 const Index& atmosphere_dim,
1233 const Index& cloudbox_checked,
1234 const Index& cloudbox_on,
1239 if( !cloudbox_checked )
1240 throw runtime_error(
"The cloudbox must be flagged to have passed a "
1241 "consistency check (cloudbox_checked=1)." );
1243 if( atmosphere_dim == 1 )
1244 {
throw runtime_error(
"No use in calling this method for 1D." ); }
1246 {
throw runtime_error(
1247 "No use in calling this method with cloudbox off." ); }
1250 {
throw runtime_error(
"The argument *nzero must be > 0." ); }
1254 const Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1255 const Index nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1257 if( atmosphere_dim == 3 )
1258 { nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1; }
1260 if( pnd_field.
npages() != np || pnd_field.
nrows() != 1 ||
1261 pnd_field.
ncols() != 1 )
1262 {
throw runtime_error(
"The input *pnd_field* is either not 1D or does not "
1263 "match pressure size of cloudbox." );}
1269 pnd_field.
resize( npart, np, nlat, nlon );
1272 for(
Index ilon=nzero; ilon<nlon-nzero; ilon++ )
1274 for(
Index ilat=nzero; ilat<nlat-nzero; ilat++ )
1276 for(
Index ip=0; ip<np; ip++ )
1278 for(
Index is=0; is<npart; is++ )
1279 { pnd_field(is,ip,ilat,ilon) = pnd_temp(is,ip,0,0); }
1297 if (lat_grid.
nelem()>1)
1312 scat_data_raw.resize(1);
1314 scat_data_raw[0].description =
" ";
1316 nlinspace(scat_data_raw[0].f_grid, 1e9, 3.848043e+13 , 5);
1317 nlinspace(scat_data_raw[0].T_grid, 0, 400, 5);
1318 nlinspace(scat_data_raw[0].za_grid, 0, 180, 5);
1319 nlinspace(scat_data_raw[0].aa_grid, 0, 360, 5);
1321 scat_data_raw[0].pha_mat_data.resize(5,5,5,1,1,1,6);
1322 scat_data_raw[0].pha_mat_data = 0.;
1323 scat_data_raw[0].ext_mat_data.resize(5,5,1,1,1);
1324 scat_data_raw[0].ext_mat_data = 0.;
1325 scat_data_raw[0].abs_vec_data.resize(5,5,1,1,1);
1326 scat_data_raw[0].abs_vec_data = 0.;
1334 const Index& atmosphere_dim,
1335 const Index& cloudbox_on,
1337 const Tensor4& massdensity_field,
1348 pnd_field.
resize(0, 0, 0, 0);
1354 Index p_cbstart = 0;
1356 Index lat_cbstart = 0;
1357 Index lat_cbend = 1;
1358 Index lon_cbstart = 0;
1359 Index lon_cbend = 1;
1361 p_cbstart = cloudbox_limits[0];
1362 p_cbend = cloudbox_limits[1]+1;
1365 if ( atmosphere_dim >= 2 )
1367 lat_cbstart = cloudbox_limits[2];
1368 lat_cbend = cloudbox_limits[3]+1;
1371 if ( atmosphere_dim == 3 )
1373 lon_cbstart = cloudbox_limits[4];
1374 lon_cbend = cloudbox_limits[5]+1;
1380 if ((p_cbend > massdensity_field.
npages()) ||
1381 (p_cbend > t_field.
npages()) ||
1382 (lat_cbend > massdensity_field.
nrows()) ||
1383 (lat_cbend > t_field.
nrows()) ||
1384 (lon_cbend > massdensity_field.
ncols()) ||
1385 (lon_cbend > t_field.
ncols()))
1388 os <<
"Cloudbox out of bounds compared to fields. "
1389 <<
"Upper limits: (p, lat, lon): "
1390 <<
"(" << p_cbend <<
", " << lat_cbend <<
", " << lon_cbend <<
"). "
1391 <<
"*massdensity_field*: "
1392 <<
"(" << massdensity_field.
npages() <<
", "
1393 << massdensity_field.
nrows() <<
", "
1394 << massdensity_field.
ncols() <<
"). "
1396 <<
"(" << t_field.
npages() <<
", "
1397 << t_field.
nrows() <<
", "
1398 << t_field.
ncols() <<
").";
1399 throw runtime_error(os.str());
1403 pnd_field.
resize ( scat_data_meta_array.
nelem(), p_cbend-p_cbstart,
1404 lat_cbend-lat_cbstart, lon_cbend-lon_cbstart );
1405 Index scat_data_start = 0;
1411 for (
Index k=0; k<part_species.
nelem(); k++ )
1420 Vector vol_unsorted ( scat_data_nelem[k], 0.0 );
1421 Vector d_max_unsorted (scat_data_nelem[k], 0.0);
1422 Vector vol ( scat_data_nelem[k], 0.0 );
1423 Vector dm ( scat_data_nelem[k], 0.0 );
1424 Vector r ( scat_data_nelem[k], 0.0 );
1425 Vector rho ( scat_data_nelem[k], 0.0 );
1426 Vector pnd ( scat_data_nelem[k], 0.0 );
1428 Vector dN ( scat_data_nelem[k], 0.0 );
1434 if ( psd_param ==
"MH97" )
1437 for (
Index i=0; i < scat_data_nelem[k]; i++ )
1440 vol_unsorted[i] = ( scat_data_meta_array[i+scat_data_start].V );
1453 for (
Index i=0; i< scat_data_nelem[k]; i++ )
1455 vol[i] = ( scat_data_meta_array[intarr[i]+scat_data_start].V );
1458 ( (6*scat_data_meta_array[intarr[i]+scat_data_start].V) /
PI ),
1461 rho[i] = scat_data_meta_array[intarr[i]+scat_data_start].density * 1000;
1464 if ( scat_data_meta_array[intarr[i]+scat_data_start].type !=
"Ice" )
1466 throw runtime_error (
"The particle phase is unequal 'Ice'.\n"
1467 "MH97 can only be applied to ice particles.\n"
1468 "Check ScatteringMetaData!" );
1474 for (
Index p=p_cbstart; p<p_cbend; p++ )
1476 for (
Index lat=lat_cbstart; lat<lat_cbend; lat++ )
1478 for (
Index lon=lon_cbstart; lon<lon_cbend; lon++ )
1486 t_field ( p, lat, lon ), rho[i] );
1503 chk_pndsum ( pnd, IWC_field ( p,lat,lon ), vol, rho, p, lat, lon, verbosity );
1507 for (
Index i = 0; i< scat_data_nelem[k]; i++ )
1509 pnd_field ( intarr[i]+scat_data_start, p-p_cbstart,
1510 lat-lat_cbstart, lon-lon_cbstart ) = pnd[i];
1520 else if ( psd_param ==
"H11" )
1525 for (
Index i=0; i < scat_data_nelem[k]; i++ )
1528 d_max_unsorted[i] = ( scat_data_meta_array[i+scat_data_start].d_max );
1535 if ( part_type ==
"IWC" )
1543 else if ( part_type ==
"Snow" )
1552 for (
Index i=0; i< scat_data_nelem[k]; i++ )
1554 vol[i]= scat_data_meta_array[intarr[i]+scat_data_start].V;
1557 dm[i] = scat_data_meta_array[intarr[i]+scat_data_start].d_max;
1560 rho[i] = scat_data_meta_array[intarr[i]+scat_data_start].density * 1000;
1563 if ( scat_data_meta_array[intarr[i]+scat_data_start].type !=
"Ice" )
1565 throw runtime_error (
1566 "The particle phase is unequal 'Ice'.\n"
1567 "H11 can only be applied to ice/snow particles.\n"
1568 "Check ScatteringMetaData!" );
1572 for (
Index p=p_cbstart; p<p_cbend; p++ )
1574 for (
Index lat=lat_cbstart; lat<lat_cbend; lat++ )
1576 for (
Index lon=lon_cbstart; lon<lon_cbend; lon++ )
1583 dN[i] =
psd_H11 ( X_field ( p, lat, lon), dm[i], t_field ( p, lat, lon ) );
1598 scale_H11 ( pnd, X_field ( p,lat,lon ), vol, rho );
1601 chk_pndsum ( pnd, X_field ( p,lat,lon ), vol, rho, p, lat, lon, verbosity );
1606 for (
Index i =0; i< scat_data_nelem[k]; i++ )
1608 pnd_field ( intarr[i]+scat_data_start, p-p_cbstart,
1609 lat-lat_cbstart, lon-lon_cbstart ) = pnd[i];
1619 else if ( psd_param ==
"liquid" )
1621 for (
Index i=0; i < scat_data_nelem[k]; i++ )
1624 vol_unsorted[i] = ( scat_data_meta_array[i+scat_data_start].V );
1636 for (
Index i=0; i< scat_data_nelem[k]; i++ )
1638 vol[i]= scat_data_meta_array[intarr[i]+scat_data_start].V;
1641 ( 6*scat_data_meta_array[intarr[i]+scat_data_start].V/
PI ),
1646 rho[i] = scat_data_meta_array[intarr[i]+scat_data_start].density * 1000;
1649 if ( scat_data_meta_array[intarr[i]+scat_data_start].type !=
"Water" )
1651 throw runtime_error (
1652 "The particle phase is unequal 'Water'.\n"
1653 "All particles must be of liquid phase to apply this PSD.\n"
1654 "Check ScatteringMetaData!" );
1659 for (
Index p=p_cbstart; p<p_cbend; p++ )
1661 for (
Index lat=lat_cbstart; lat<lat_cbend; lat++ )
1663 for (
Index lon=lon_cbstart; lon<lon_cbend; lon++ )
1670 dN[i] =
LWCtopnd ( LWC_field ( p,lat,lon ), r[i] );
1690 chk_pndsum ( pnd, LWC_field ( p,lat,lon ), vol, rho, p, lat, lon, verbosity );
1696 for (
Index i =0; i< scat_data_nelem[k]; i++ )
1698 pnd_field ( intarr[i]+scat_data_start, p-p_cbstart,
1699 lat-lat_cbstart, lon-lon_cbstart ) = pnd[i];
1716 scat_data_start = scat_data_start + scat_data_nelem[k];