Changeset 5613 for LMDZ6/branches
- Timestamp:
- Apr 14, 2025, 6:16:16 PM (3 months ago)
- Location:
- LMDZ6/branches/contrails
- Files:
-
- 5 edited
-
DefLists/file_def_histday_lmdz.xml (modified) (1 diff)
-
DefLists/file_def_histhf_lmdz.xml (modified) (1 diff)
-
libf/phylmd/lmdz_aviation.f90 (modified) (1 diff)
-
libf/phylmd/lmdz_lscp_condensation.f90 (modified) (15 diffs)
-
libf/phylmd/lmdz_lscp_precip.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/contrails/DefLists/file_def_histday_lmdz.xml
r5573 r5613 504 504 <!-- VARS 3D --> 505 505 <field_group operation="average" grid_ref="grid_out_presnivs"> 506 <field field_ref="flightdist" level="1" />507 <field field_ref="cfseri" level="1" />508 <field field_ref="qissr" level="1" />509 <field field_ref="qcld" level="1" />510 <field field_ref="subfra" level="1" />511 <field field_ref="issrfra" level="1" />512 <field field_ref="dqvcmix" level="1" />513 <field field_ref="dqimix" level="1" />514 <field field_ref="dcfmix" level="1" />515 506 <field field_ref="CO2" level="5" /> <!-- Added PC --> 516 507 <field field_ref="dCO2_vdf" level="5" /> -
LMDZ6/branches/contrails/DefLists/file_def_histhf_lmdz.xml
r5573 r5613 479 479 <!-- VARS 3D --> 480 480 <field_group operation="average" grid_ref="grid_out_presnivs"> 481 <field field_ref="flightdist" level="1" operation="instant" />482 481 <field field_ref="CO2" level="10" /> <!-- Added PC --> 483 482 <field field_ref="tke" level="10" /> -
LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
r5609 r5613 174 174 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 175 175 pdf_x = qcritcont(i) / qsatl(i) * 100. 176 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 177 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 178 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 179 pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i) 180 pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) & 181 + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100. 176 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 177 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 178 pdf_fra_above_qcritcont = 0. 179 pdf_q_above_qcritcont = 0. 180 ELSEIF ( pdf_y .LT. -10. ) THEN 181 pdf_fra_above_qcritcont = 1. 182 pdf_q_above_qcritcont = qclr(i) / clrfra(i) 183 ELSE 184 pdf_y = EXP( pdf_y ) 185 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 186 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 187 pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i) 188 pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) & 189 + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100. 190 ENDIF 182 191 183 192 pdf_x = qsat(i) / qsatl(i) * 100. 184 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 185 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 186 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 187 pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i) 188 pdf_q_above_qsat = ( pdf_e3 * clrfra(i) & 189 + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100. 193 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 194 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 195 pdf_fra_above_qsat = 0. 196 pdf_q_above_qsat = 0. 197 ELSEIF ( pdf_y .LT. -10. ) THEN 198 pdf_fra_above_qsat = 1. 199 pdf_q_above_qsat = qclr(i) / clrfra(i) 200 ELSE 201 pdf_y = EXP( pdf_y ) 202 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 203 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 204 pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i) 205 pdf_q_above_qsat = ( pdf_e3 * clrfra(i) & 206 + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100. 207 ENDIF 190 208 191 209 potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat) -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
r5609 r5613 239 239 ! 240 240 ! for deposition / sublimation 241 REAL :: pres_sat, kappa_depsub, tau _depsub241 REAL :: pres_sat, kappa_depsub, tauinv_depsub 242 242 REAL :: air_thermal_conduct, water_vapor_diff 243 243 REAL :: N_ice … … 264 264 REAL :: pdf_fra_above_nuc, pdf_q_above_nuc 265 265 REAL :: pdf_fra_above_lim, pdf_q_above_lim 266 REAL :: pdf_fra_below_lim , pdf_q_below_lim266 REAL :: pdf_fra_below_lim 267 267 REAL :: mixed_fraction 268 268 ! … … 545 545 !--If the cloud is initially supersaturated 546 546 !--Exact explicit integration (qvc exact, qice explicit) 547 tau _depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub )547 tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub 548 548 ELSE 549 549 !--If the cloud is initially subsaturated 550 550 !--Exact explicit integration (qvc exact, qice explicit) 551 551 !--Same but depo_coef_cirrus = 1 552 tau _depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub )552 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub 553 553 ENDIF ! qvapincld .GT. qsat 554 qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / tau_depsub )554 qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub ) 555 555 !--If all the ice is sublimated 556 556 IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0. … … 641 641 !--New PDF 642 642 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 643 rhl_clr = MAX(0., MIN(1 000., rhl_clr))643 rhl_clr = MAX(0., MIN(150., rhl_clr)) 644 644 645 645 !--Calculation of the properties of the PDF … … 662 662 pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) 663 663 pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3 664 pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows 664 665 665 666 !IF ( ok_warm_cloud ) THEN … … 671 672 672 673 pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i)) 673 pdf_scale(i) = MAX(eps, pdf_std / SQRT( & 674 !--Barrier to avoid overflows 675 pdf_scale(i) = MAX(0.01, pdf_std / SQRT( & 674 676 GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 )) 675 677 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) … … 680 682 681 683 pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. 682 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 683 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 684 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 685 pdf_fra_above_nuc = EXP( - pdf_y ) 686 pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100. 684 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 685 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 686 pdf_fra_above_nuc = 0. 687 pdf_q_above_nuc = 0. 688 ELSEIF ( pdf_y .LT. -10. ) THEN 689 pdf_fra_above_nuc = 1. 690 pdf_q_above_nuc = qclr(i) / clrfra(i) 691 ELSE 692 pdf_y = EXP( pdf_y ) 693 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 694 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 695 pdf_fra_above_nuc = EXP( - pdf_y ) 696 pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100. 697 ENDIF 687 698 688 699 IF ( pdf_fra_above_nuc .GT. eps ) THEN … … 701 712 qvapincld = gamma_cond(i) * qsat(i) 702 713 qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) 703 tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub ) 704 qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / 2. / tau_depsub ) 714 tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub 715 qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) & 716 * EXP( - dtime / 2. * tauinv_depsub ) 705 717 ELSE 706 718 !--We keep the saturation adjustment hypothesis, and the vapor in the … … 720 732 qclr(i) = qclr(i) - dqt_con 721 733 722 ENDIF ! ok_warm_cloud, cf_cond .GT. eps 723 ENDIF ! ( 1. - cldfra(i) ) .GT. eps 734 ENDIF ! pdf_fra_above_nuc .GT. eps 735 ELSE 736 !--Default value for the clear sky distribution: homogeneous distribution 737 pdf_alpha(i) = 1. 738 pdf_gamma(i) = 1. 739 pdf_scale(i) = 0.01 740 ENDIF ! clrfra(i) .GT. eps 724 741 725 742 … … 826 843 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 827 844 qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix ) 828 tau _depsub = 1. / ( qiceinmix**(1./3.) * kappa_depsub )829 qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime / tau_depsub ) )845 tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub 846 qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) ) 830 847 qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) & 831 848 - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix 832 849 ELSE 833 !--NB. if tau_depsub = 0 , we get the same result as above850 !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above 834 851 qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & 835 852 - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix … … 851 868 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 852 869 pdf_x = qvapinclr_lim / qsatl(i) * 100. 853 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 854 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 855 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 856 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 857 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 858 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 870 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 871 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 872 pdf_fra_above_lim = 0. 873 pdf_q_above_lim = 0. 874 ELSEIF ( pdf_y .LT. -10. ) THEN 875 pdf_fra_above_lim = 1. 876 pdf_q_above_lim = qclr(i) / clrfra(i) 877 ELSE 878 pdf_y = EXP( pdf_y ) 879 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 880 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 881 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 882 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 883 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 884 ENDIF 859 885 860 886 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 861 pdf_q_below_lim = qclr(i) - pdf_q_above_lim862 887 863 888 !--sigma_mix is the ratio of the surroundings of the clouds where mixing … … 1000 1025 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1001 1026 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1002 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 1003 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1004 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1005 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1006 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1007 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1027 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 1028 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1029 pdf_fra_above_lim = 0. 1030 pdf_q_above_lim = 0. 1031 ELSEIF ( pdf_y .LT. -10. ) THEN 1032 pdf_fra_above_lim = 1. 1033 pdf_q_above_lim = qclr(i) / clrfra(i) 1034 ELSE 1035 pdf_y = EXP( pdf_y ) 1036 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1037 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1038 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1039 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1040 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1041 ENDIF 1008 1042 1009 1043 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim … … 1086 1120 qiceincld = qice_sedim(i) / cldfra_above(i) 1087 1121 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 1088 tau _depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub )1089 qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime / tau_depsub ) )1122 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub 1123 qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) 1090 1124 ELSE 1091 1125 qvapinclr_lim = qsat(i) - qiceincld … … 1095 1129 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1096 1130 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1097 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 1098 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1099 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1100 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1101 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1102 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1131 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 1132 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1133 pdf_fra_above_lim = 0. 1134 pdf_q_above_lim = 0. 1135 ELSEIF ( pdf_y .LT. -10. ) THEN 1136 pdf_fra_above_lim = 1. 1137 pdf_q_above_lim = qclr(i) / clrfra(i) 1138 ELSE 1139 pdf_y = EXP( pdf_y ) 1140 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1141 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1142 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1143 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1144 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1145 ENDIF 1103 1146 1104 1147 IF ( pdf_fra_above_lim .GT. eps ) THEN … … 1140 1183 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1141 1184 pdf_x = qsat(i) / qsatl(i) * 100. 1142 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 1143 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1144 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1145 issrfra(i) = EXP( - pdf_y ) * clrfra(i) 1146 qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. 1185 pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) 1186 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1187 issrfra(i) = 0. 1188 qissr(i) = 0. 1189 ELSEIF ( pdf_y .LT. -10. ) THEN 1190 issrfra(i) = 1. 1191 qissr(i) = qclr(i) / clrfra(i) 1192 ELSE 1193 pdf_y = EXP( pdf_y ) 1194 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1195 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1196 issrfra(i) = EXP( - pdf_y ) * clrfra(i) 1197 qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. 1198 ENDIF 1147 1199 ELSE 1148 1200 issrfra(i) = 0. -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90
r5609 r5613 697 697 698 698 iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3 699 tempc = zt(i) - RTT ! celcius temp 700 ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 701 IF ( tempc .GE. -60. ) THEN 702 icevelo = EXP( 1.9714 + 0.00216078 * tempc & 703 - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) 704 ELSE 705 icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 706 ENDIF 707 icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s 699 !tempc = zt(i) - RTT ! celcius temp 700 !! so-called 'empirical parameterization' in Stubenrauch et al. 2019 701 !IF ( tempc .GE. -60. ) THEN 702 ! icevelo = EXP( 1.9714 + 0.00216078 * tempc & 703 ! - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) 704 !ELSE 705 ! icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 706 !ENDIF 707 !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s 708 709 icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16 708 710 709 711 qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i)) 710 712 !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two 711 713 !--times in the same timestep) 712 !qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)714 qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime) 713 715 !--Add tendencies 714 716 zoliqi(i) = zoliqi(i) - qice_sedim … … 1964 1966 rho = pplay(i) / temp(i) / RD 1965 1967 iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3 1966 tempc = temp(i) - RTT ! celcius temp 1967 ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 1968 IF ( tempc .GE. -60. ) THEN 1969 icevelo = EXP( 1.9714 + 0.00216078 * tempc & 1970 - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) 1971 ELSE 1972 icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 1973 ENDIF 1974 icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s 1968 !tempc = temp(i) - RTT ! celcius temp 1969 !! so-called 'empirical parameterization' in Stubenrauch et al. 2019 1970 !IF ( tempc .GE. -60. ) THEN 1971 ! icevelo = EXP( 1.9714 + 0.00216078 * tempc & 1972 ! - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) 1973 !ELSE 1974 ! icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 1975 !ENDIF 1976 !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s 1977 1978 icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16 1975 1979 1976 1980 dz = ( paprsdn(i) - paprsup(i) ) / RG / rho 1977 1981 qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz) 1982 !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two 1983 !--times in the same timestep) 1984 qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime) 1978 1985 !--Add tendencies 1979 1986 qice(i) = qice(i) - qice_sedim
Note: See TracChangeset
for help on using the changeset viewer.
