| 1 | % This script plots horizontal migration distances predicted by |
|---|
| 2 | % SPARTACUS's "explicit entrapment" option, but it requires (1) the |
|---|
| 3 | % whole code to have been compiled with "make |
|---|
| 4 | % PRINT_ENTRAPMENT_DATA=1", and (2) "make i3rc_print_entrapment" to |
|---|
| 5 | % have been run in this directory. This produces the fort.10[12] files |
|---|
| 6 | % containing the migration distances. |
|---|
| 7 | |
|---|
| 8 | xdata = load('fort.101'); |
|---|
| 9 | fdata = load('fort.102'); |
|---|
| 10 | in = loadnc('i3rc_mls_cumulus_sza.nc'); |
|---|
| 11 | out= loadnc('i3rc_mls_cumulus_mono_entr_out.nc'); |
|---|
| 12 | |
|---|
| 13 | z=in.height_hl(1:end-1,1)./1000; |
|---|
| 14 | zh=in.height_hl(:,1)./1000; |
|---|
| 15 | [nlev,ncol] = size(in.cloud_fraction); |
|---|
| 16 | |
|---|
| 17 | x_region_direct = flip(reshape(xdata(:,3:5),nlev,ncol,3),1); |
|---|
| 18 | x_region_diffuse = flip(reshape(xdata(:,6:8),nlev,ncol,3),1); |
|---|
| 19 | n_scat_region = flip(reshape(xdata(:,9:11),nlev,ncol,3),1); |
|---|
| 20 | |
|---|
| 21 | x_diffuse_scaling = sqrt(2.0 * (n_scat_region+exp(-n_scat_region)-1.0) + 1.0e-16) ./ (n_scat_region+1.0e-8); |
|---|
| 22 | |
|---|
| 23 | tot_region_diffuse = sqrt(2.0).*x_region_diffuse.*x_diffuse_scaling; |
|---|
| 24 | tot_region_direct = sqrt(x_region_direct.^2 + (x_region_diffuse.*x_diffuse_scaling).^2); |
|---|
| 25 | |
|---|
| 26 | flux_region_direct = reshape(fdata(:,3:5),nlev,ncol,3); |
|---|
| 27 | flux_region_dn_diffuse = reshape(fdata(:,6:8),nlev,ncol,3); |
|---|
| 28 | |
|---|
| 29 | flux_direct = squeeze(sum(flux_region_direct,3)); |
|---|
| 30 | flux_dn_diffuse = squeeze(sum(flux_region_dn_diffuse,3)); |
|---|
| 31 | |
|---|
| 32 | ind = find(flux_region_dn_diffuse(:,1,1) < 1.0e-8); |
|---|
| 33 | flux_region_dn_diffuse(ind,:,1) = NaN; |
|---|
| 34 | |
|---|
| 35 | xdirect = squeeze(sum(flux_region_direct.*tot_region_direct,3) ... |
|---|
| 36 | ./sum(flux_region_direct,3)); |
|---|
| 37 | xdiffuse = squeeze(sum(flux_region_dn_diffuse.*tot_region_diffuse,3) ... |
|---|
| 38 | ./sum(flux_region_dn_diffuse,3)); |
|---|
| 39 | |
|---|
| 40 | sza = acosd(in.cos_solar_zenith_angle); |
|---|
| 41 | |
|---|
| 42 | isza = [1 16 31 41]; |
|---|
| 43 | isza = [1 23 36]; |
|---|
| 44 | |
|---|
| 45 | xmax = 12; |
|---|
| 46 | zmax = 4; |
|---|
| 47 | |
|---|
| 48 | figure(2) |
|---|
| 49 | clf |
|---|
| 50 | for ii = 1:length(isza) |
|---|
| 51 | subplot(2,2,ii) |
|---|
| 52 | %plot(squeeze(x_region_direct(:,isza(ii),:)),z); |
|---|
| 53 | %plot(squeeze(x_region_diffuse(:,isza(ii),:)),z,'--'); |
|---|
| 54 | |
|---|
| 55 | plot(xdirect(:,isza(ii))./1000,z,'r') |
|---|
| 56 | hold on |
|---|
| 57 | plot(xdiffuse(:,isza(ii))./1000,z,'b'); |
|---|
| 58 | ylim([0 zmax]); |
|---|
| 59 | xlim([0 xmax]) |
|---|
| 60 | xlabel('Horizontal distance travelled (km)'); |
|---|
| 61 | ylabel('Height (km)'); |
|---|
| 62 | title(['(' 'a'-1+ii ') SZA=' num2str(sza(isza(ii))) '\circ']) |
|---|
| 63 | if ii == 1 |
|---|
| 64 | legend('Direct','Diffuse','location','southeast') |
|---|
| 65 | end |
|---|
| 66 | grid on |
|---|
| 67 | end |
|---|
| 68 | |
|---|
| 69 | subplot(2,2,4) |
|---|
| 70 | zz = 0.5.*(zh(1:end-1)+zh(2:end)); |
|---|
| 71 | plot(in.cloud_fraction, zz,'k'); |
|---|
| 72 | ylim([0 zmax]) |
|---|
| 73 | xlabel('Cloud fraction'); |
|---|
| 74 | ylabel('Height (km)'); |
|---|
| 75 | grid on |
|---|
