Changeset 5413 for LMDZ6


Ignore:
Timestamp:
Dec 17, 2024, 9:48:07 AM (2 months ago)
Author:
evignon
Message:

fin de l'externalisation des precips dans lscp.

  1. Borella, E Vignon
Location:
LMDZ6/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90

    r5412 r5413  
    105105USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
    106106USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
    107 USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb
    108 USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
     107USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
    109108USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
    110 USE lmdz_lscp_precip, ONLY : histprecip_precld, poprecip_precld, poprecip_postcld
     109USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
     110USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
    111111
    112112! Use du module lmdz_lscp_ini contenant les constantes
    113113USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
    114 USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, DDT0, rain_int_min
    115 USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
    116 USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
    117 USE lmdz_lscp_ini, ONLY : cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
    118 USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
    119 USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
    120 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
     114USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
     115USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
     116USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
     117USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
     118USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
    121119USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
    122120USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
    … …  
    154152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
    155153  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
    156   REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
    157   REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
     154  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
     155  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
    158156
    159157  ! INPUT/OUTPUT variables
    … …  
    268266  !----------------
    269267  REAL,DIMENSION(klon) :: qice_ini
    270   REAL :: zct, zcl,zexpo
    271268  REAL, DIMENSION(klon,klev) :: ctot
    272269  REAL, DIMENSION(klon,klev) :: ctot_vol
    273270  REAL, DIMENSION(klon) :: zqs, zdqs
    274   REAL :: zdelta, zcor, zcvm5
     271  REAL :: zdelta
    275272  REAL, DIMENSION(klon) :: zdqsdT_raw
    276273  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
    … …  
    278275  REAL :: num,denom
    279276  REAL :: cste
    280   REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
    281   REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
    282   REAL :: erf
    283277  REAL, DIMENSION(klon) :: zfice_th
    284278  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
    285   REAL, DIMENSION(klon) :: zrfl
    286   REAL, DIMENSION(klon) :: zifl, ziflprev
     279  REAL, DIMENSION(klon) :: zrfl, zifl
    287280  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
    288281  REAL, DIMENSION(klon) :: zoliql, zoliqi
    289282  REAL, DIMENSION(klon) :: zt
    290   REAL, DIMENSION(klon,klev) :: zrho
    291   REAL, DIMENSION(klon) :: zdz,iwc
    292   REAL :: zchau,zfroi
    293283  REAL, DIMENSION(klon) :: zfice,zneb
    294   REAL :: zrain,zsnow,zprecip
    295284  REAL, DIMENSION(klon) :: dzfice
    296285  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
    297   REAL :: zsolid
    298286  REAL, DIMENSION(klon) :: qtot, qzero
    299   REAL :: smallestreal
    300   !  Variables for Bergeron process
    301   REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
    302   REAL, DIMENSION(klon) :: zqpreci, zqprecl
    303287  ! Variables precipitation energy conservation
    304288  REAL, DIMENSION(klon) :: zmqc
    … …  
    309293  REAL, DIMENSION(klon) :: ztupnew
    310294  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
    311   REAL :: zm_solid         ! for liquid -> solid conversion
    312295  REAL, DIMENSION(klon) :: zrflclr, zrflcld
    313   REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
    314   REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
    315296  REAL, DIMENSION(klon) :: ziflclr, ziflcld
    316297  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
    317   REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
    318   REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
    319   REAL, DIMENSION(klon,klev) :: velo
    320   REAL :: vr, ffallv
     298  REAL, DIMENSION(klon) :: tot_zneb
    321299  REAL :: qlmpc, qimpc, rnebmpc
    322   REAL, DIMENSION(klon,klev) :: radocondi, radocondl
    323   REAL :: effective_zneb
    324300  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
    325301  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
     302
     303  ! for quantity of condensates seen by radiation
     304  REAL, DIMENSION(klon) :: zradocond, zradoice
     305  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
    326306 
    327307  ! for condensation and ice supersaturation
    … …  
    336316  REAL :: min_qParent, min_ratio
    337317
    338   INTEGER i, k, n, kk, iter
     318  INTEGER i, k, kk, iter
    339319  INTEGER, DIMENSION(klon) :: n_i
    340320  INTEGER ncoreczq
    … …  
    365345ctot_vol(1:klon,1:klev)=0.0
    366346rneblsvol(1:klon,1:klev)=0.0
    367 smallestreal=1.e-9
    368347znebprecip(:)=0.0
    369348znebprecipclr(:)=0.0
    … …  
    398377rain(:) = 0.0
    399378snow(:) = 0.0
    400 zoliq(:)=0.0
    401379zfice(:)=0.0
    402380dzfice(:)=0.0
    403381zfice_turb(:)=0.0
    404382dzfice_turb(:)=0.0
    405 zqprecl(:)=0.0
    406 zqpreci(:)=0.0
    407383zrfl(:) = 0.0
    408384zifl(:) = 0.0
    409 ziflprev(:)=0.0
    410385zneb(:) = seuil_neb
    411386zrflclr(:) = 0.0
    … …  
    414389ziflcld(:) = 0.0
    415390tot_zneb(:) = 0.0
    416 tot_znebn(:) = 0.0
    417 d_tot_zneb(:) = 0.0
    418391qzero(:) = 0.0
    419392zdistcltop(:)=0.0
    … …  
    913886    ! ----------------------------------------------------------------
    914887
     888    ! Initiate the quantity of liquid and solid condensates
     889    ! Note that in the following, zcond is the total amount of condensates
     890    ! including newly formed precipitations (i.e., condensates formed by the
     891    ! condensation process), while zoliq is the total amount of condensates in
     892    ! the cloud (i.e., on which precipitation processes have applied)
     893    DO i = 1, klon
     894      zoliq(i)  = zcond(i)
     895      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
     896      zoliqi(i) = zcond(i) * zfice(i)
     897      ! c_iso : initialisation of zoliq* also for isotopes
     898    ENDDO
     899
    915900    !================================================================
    916901    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
    917902    IF (ok_poprecip) THEN
    918 
    919       DO i = 1, klon
    920         zoliql(i) = zcond(i) * ( 1. - zfice(i) )
    921         zoliqi(i) = zcond(i) * zfice(i)
    922       ENDDO
    923903
    924904      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    … …  
    934914                            )
    935915
    936       DO i = 1, klon
    937         zoliq(i) = zoliql(i) + zoliqi(i)
    938         IF ( zoliq(i) .GT. 0. ) THEN
    939                 zfice(i) = zoliqi(i)/zoliq(i)
    940         ELSE 
    941                 zfice(i) = 0.0
    942         ENDIF
    943 
    944         ! calculation of specific content of condensates seen by radiative scheme
    945         IF (ok_radocond_snow) THEN
    946            radocond(i,k) = zoliq(i)
    947            radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
    948            radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k)
    949         ELSE
    950            radocond(i,k) = zoliq(i)
    951            radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
    952            radocondi(i,k)= radocond(i,k)*zfice(i)
    953         ENDIF
    954       ENDDO
    955 
    956916    !================================================================
    957917    ELSE
    958918
    959     ! LTP:
    960     IF (iflag_evap_prec .GE. 4) THEN
    961 
    962         !Partitionning between precipitation coming from clouds and that coming from CS
    963 
    964         !0) Calculate tot_zneb, total cloud fraction above the cloud
    965         !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
    966        
    967         DO i=1, klon
    968                 tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
    969                         /(1.-min(zneb(i),1.-smallestreal))
    970                 d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
    971                 tot_zneb(i) = tot_znebn(i)
    972 
    973 
    974                 !1) Cloudy to clear air
    975                 d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
    976                 IF (znebprecipcld(i) .GT. 0.) THEN
    977                         d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
    978                         d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
    979                 ELSE
    980                         d_zrfl_cld_clr(i) = 0.
    981                         d_zifl_cld_clr(i) = 0.
    982                 ENDIF
    983 
    984                 !2) Clear to cloudy air
    985                 d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
    986                 IF (znebprecipclr(i) .GT. 0) THEN
    987                         d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
    988                         d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
    989                 ELSE
    990                         d_zrfl_clr_cld(i) = 0.
    991                         d_zifl_clr_cld(i) = 0.
    992                 ENDIF
    993 
    994                 !Update variables
    995                 znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
    996                 znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
    997                 zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
    998                 ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
    999                 zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
    1000                 ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
    1001 
    1002                 ! c_iso  do the same thing for isotopes precip
    1003         ENDDO
    1004     ENDIF
    1005 
    1006 
    1007     ! Autoconversion
    1008     ! -------------------------------------------------------------------------------
    1009 
    1010 
    1011     ! Initialisation of zoliq and radocond variables
    1012 
    1013     DO i = 1, klon
    1014             zoliq(i) = zcond(i)
    1015             zoliqi(i)= zoliq(i)*zfice(i)
    1016             zoliql(i)= zoliq(i)*(1.-zfice(i))
    1017             ! c_iso : initialisation of zoliq* also for isotopes
    1018             zrho(i,k)  = pplay(i,k) / zt(i) / RD
    1019             zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
    1020             iwc(i)   = 0.
    1021             zneb(i)  = MAX(rneb(i,k),seuil_neb)
    1022             radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
    1023             radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
    1024             radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
    1025     ENDDO
    1026 
    1027        
    1028     DO n = 1, niter_lscp
    1029 
    1030         DO i=1,klon
    1031             IF (rneb(i,k).GT.0.0) THEN
    1032                 iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
    1033             ENDIF
    1034         ENDDO
    1035 
    1036         CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
    1037 
    1038         DO i = 1, klon
    1039 
    1040             IF (rneb(i,k).GT.0.0) THEN
    1041 
    1042                 ! Initialization of zrain, zsnow and zprecip:
    1043                 zrain=0.
    1044                 zsnow=0.
    1045                 zprecip=0.
    1046                 ! c_iso same init for isotopes. Externalisation?
    1047 
    1048                 IF (zneb(i).EQ.seuil_neb) THEN
    1049                     zprecip = 0.0
    1050                     zsnow = 0.0
    1051                     zrain= 0.0
    1052                 ELSE
    1053 
    1054                     IF (ptconv(i,k)) THEN ! if convective point
    1055                         zcl=cld_lc_con
    1056                         zct=1./cld_tau_con
    1057                         zexpo=cld_expo_con
    1058                         ffallv=ffallv_con
    1059                     ELSE
    1060                         zcl=cld_lc_lsc
    1061                         zct=1./cld_tau_lsc
    1062                         zexpo=cld_expo_lsc
    1063                         ffallv=ffallv_lsc
    1064                     ENDIF
    1065 
    1066 
    1067                     ! if vertical heterogeneity is taken into account, we use
    1068                     ! the "true" volume fraction instead of a modified
    1069                     ! surface fraction (which is larger and artificially
    1070                     ! reduces the in-cloud water).
    1071 
    1072                     ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
    1073                     ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
    1074                     !.........................................................
    1075                     IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
    1076 
    1077                     ! if vertical heterogeneity is taken into account, we use
    1078                     ! the "true" volume fraction instead of a modified
    1079                     ! surface fraction (which is larger and artificially
    1080                     ! reduces the in-cloud water).
    1081                         effective_zneb=ctot_vol(i,k)
    1082                     ELSE
    1083                         effective_zneb=zneb(i)
    1084                     ENDIF
    1085 
    1086 
    1087                     IF (iflag_autoconversion .EQ. 2) THEN
    1088                     ! two-steps resolution with niter_lscp=1 sufficient
    1089                     ! we first treat the second term (with exponential) in an explicit way
    1090                     ! and then treat the first term (-q/tau) in an exact way
    1091                         zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
    1092                         *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
    1093                     ELSE
    1094                     ! old explicit resolution with subtimesteps
    1095                         zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
    1096                         *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
    1097                     ENDIF
    1098 
    1099 
    1100                     ! Ice water quantity to remove (Zender & Kiehl, 1997)
    1101                     ! dqice/dt=1/rho*d(rho*wice*qice)/dz
    1102                     !....................................
    1103                     IF (iflag_autoconversion .EQ. 2) THEN
    1104                     ! exact resolution, niter_lscp=1 is sufficient but works only
    1105                     ! with iflag_vice=0
    1106                        IF (zoliqi(i) .GT. 0.) THEN
    1107                           zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
    1108                           +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
    1109                        ELSE
    1110                           zfroi=0.
    1111                        ENDIF
    1112                     ELSE
    1113                     ! old explicit resolution with subtimesteps
    1114                        zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k)
    1115                     ENDIF
    1116 
    1117                     zrain   = MIN(MAX(zchau,0.0),zoliql(i))
    1118                     zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
    1119                     zprecip = MAX(zrain + zsnow,0.0)
    1120 
    1121                 ENDIF
    1122 
    1123 
    1124                 IF (iflag_autoconversion .GE. 1) THEN
    1125                    ! debugged version with phase conservation through the autoconversion process
    1126                    zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
    1127                    zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
    1128                    zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
    1129                 ELSE
    1130                    ! bugged version with phase resetting
    1131                    zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
    1132                    zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
    1133                    zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
    1134                 ENDIF
    1135 
    1136                 !--Diagnostics
    1137                 dqrauto(i,k) = dqrauto(i,k) + zrain / dtime
    1138                 dqsauto(i,k) = dqsauto(i,k) + zsnow / dtime
    1139 
    1140 
    1141                 ! c_iso: call isotope_conversion (for readibility) or duplicate
    1142 
    1143                 radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
    1144                 radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
    1145                 radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
    1146 
    1147             ENDIF ! rneb >0
    1148 
    1149         ENDDO  ! i = 1,klon
    1150 
    1151     ENDDO      ! n = 1,niter
    1152 
    1153     ! Precipitation flux calculation
    1154 
    1155     DO i = 1, klon
    1156            
    1157             IF (iflag_evap_prec.GE.4) THEN
    1158                 ziflprev(i)=ziflcld(i)
    1159             ELSE
    1160                 ziflprev(i)=zifl(i)*zneb(i)
    1161             ENDIF
    1162 
    1163             IF (rneb(i,k) .GT. 0.0) THEN
    1164 
    1165                ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
    1166                ! If T<0C, liquid precip are converted into ice, which leads to an increase in
    1167                ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
    1168                ! taken into account through a linearization of the equations and by approximating
    1169                ! the liquid precipitation process with a "threshold" process. We assume that
    1170                ! condensates are not modified during this operation. Liquid precipitation is
    1171                ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
    1172                ! Water vapor increases as well
    1173                ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
    1174 
    1175                IF (ok_bug_phase_lscp) THEN
    1176                     zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
    1177                     zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
    1178                ELSE
    1179                     zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i)
    1180                     zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i)
    1181                ENDIF
    1182                     zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
    1183                     coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
    1184                     ! Computation of DT if all the liquid precip freezes
    1185                     DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
    1186                     ! T should not exceed the freezing point
    1187                     ! that is Delta > RTT-zt(i)
    1188                     DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
    1189                     zt(i)  = zt(i) + DeltaT
    1190                     ! water vaporization due to temp. increase
    1191                     Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
    1192                     ! we add this vaporized water to the total vapor and we remove it from the precip
    1193                     zq(i) = zq(i) +  Deltaq
    1194                     ! The three "max" lines herebelow protect from rounding errors
    1195                     zcond(i) = max( zcond(i)- Deltaq, 0. )
    1196                     ! liquid precipitation converted to ice precip
    1197                     Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
    1198                     zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
    1199                     ! iced water budget
    1200                     zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
    1201 
    1202                     ! c_iso : mv here condensation of isotopes + redispatchage en precip
    1203 
    1204                 IF (iflag_evap_prec.GE.4) THEN
    1205                     zrflcld(i) = zrflcld(i)+zqprecl(i) &
    1206                     *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    1207                     ziflcld(i) = ziflcld(i)+ zqpreci(i) &
    1208                         *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    1209                     znebprecipcld(i) = rneb(i,k)
    1210                     zrfl(i) = zrflcld(i) + zrflclr(i)
    1211                     zifl(i) = ziflcld(i) + ziflclr(i)
    1212                 ELSE
    1213                     zrfl(i) = zrfl(i)+ zqprecl(i)        &
    1214                     *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    1215                     zifl(i) = zifl(i)+ zqpreci(i)        &
    1216                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    1217                 ENDIF
    1218                 ! c_iso : same for isotopes
    1219  
    1220            ENDIF ! rneb>0
    1221 
    1222    ENDDO
    1223 
    1224     ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
    1225     ! if iflag_evap_prec>=4
    1226     IF (iflag_evap_prec.GE.4) THEN
    1227 
    1228         DO i=1,klon
    1229 
    1230             IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
    1231                 znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
    1232                 (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
    1233             ELSE
    1234                 znebprecipclr(i)=0.0
    1235             ENDIF
    1236 
    1237             IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
    1238                 znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
    1239                 (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
    1240             ELSE
    1241                 znebprecipcld(i)=0.0
    1242             ENDIF
    1243         ENDDO
    1244 
    1245     ENDIF
    1246 
     919      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
     920                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
     921                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
     922                            rneb(:,k), znebprecipclr, znebprecipcld, &
     923                            zneb, tot_zneb, zrho_up, zvelo_up, &
     924                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
     925                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
     926                            )
    1247927
    1248928    ENDIF ! ok_poprecip
    … …  
    1250930    ! End of precipitation processes after cloud formation
    1251931    ! ----------------------------------------------------
    1252 
    1253932
    1254933    !----------------------------------------------------------------------
    … …  
    1256935    !----------------------------------------------------------------------
    1257936
    1258     ! Calculation of the concentration of condensates seen by the radiation scheme
    1259     ! for liquid phase, we take radocondl
    1260     ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
    1261     ! we recalculate radocondi to account for contributions from the precipitation flux
    1262     ! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
    1263 
    1264     IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN
    1265         ! for the solid phase (crystals + snowflakes)
    1266         ! we recalculate radocondi by summing
    1267         ! the ice content calculated in the mesh
    1268         ! + the contribution of the non-evaporated snowfall
    1269         ! from the overlying layer
    1270         DO i=1,klon
    1271            IF (ziflprev(i) .NE. 0.0) THEN
    1272               radocondi(i,k)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
    1273            ELSE
    1274               radocondi(i,k)=zoliqi(i)+zqpreci(i)
    1275            ENDIF
    1276            radocond(i,k)=radocondl(i,k)+radocondi(i,k)
    1277         ENDDO
    1278     ENDIF
    1279 
    1280     ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
    1281937    DO i=1,klon
    1282         IF (radocond(i,k) .GT. 0.) THEN
    1283             radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
    1284         ENDIF
     938
     939      IF (ok_poprecip) THEN
     940        IF (ok_radocond_snow) THEN
     941          radocond(i,k) = zoliq(i)
     942          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
     943        ELSE
     944          radocond(i,k) = zoliq(i)
     945          zradoice(i) = zoliqi(i)
     946        ENDIF
     947      ELSE
     948        radocond(i,k) = zradocond(i)
     949      ENDIF
     950
     951      ! calculate the percentage of ice in "radocond" so cloud+precip seen
     952      ! by the radiation scheme
     953      IF (radocond(i,k) .GT. 0.) THEN
     954        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
     955      ENDIF
    1285956    ENDDO
    1286957
    … …  
    14081079
    14091080
    1410 ENDDO
     1081ENDDO ! loop on k from top to bottom
    14111082
    14121083

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90

    r5411 r5413  
    295295    ELSE
    296296      ! if no precip, we reinitialize the cloud fraction used for the precip to 0
    297       ! AB note that this assignment is useless, as znebprecip is not re-used
    298297      znebprecip(i)=0.
    299298
    … …  
    305304
    306305END SUBROUTINE histprecip_precld
     306
     307
     308SUBROUTINE histprecip_postcld( &
     309        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
     310        zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
     311        rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
     312        zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
     313        zradocond, zradoice, dqrauto, dqsauto &
     314        )
     315
     316USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
     317USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, &
     318                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, &
     319                          seuil_neb, rain_int_min, cice_velo, dice_velo
     320USE lmdz_lscp_ini, ONLY : iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, &
     321                          iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
     322                          niter_lscp
     323USE lmdz_lscp_tools, ONLY : fallice_velocity
     324
     325IMPLICIT NONE
     326
     327
     328INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
     329REAL,    INTENT(IN)                     :: dtime          !--time step [s]
     330LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
     331
     332REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
     333REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
     334REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
     335REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
     336LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
     337REAL,    INTENT(IN),    DIMENSION(klon) :: zdqsdT_raw     !--derivative of qsat w.r.t. temperature times L/Cp [SI]
     338
     339REAL,    INTENT(INOUT), DIMENSION(klon) :: zt             !--current temperature [K]
     340REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity [kg/kg]
     341REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliq          !--current liquid+ice water specific humidity [kg/kg]
     342REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliql         !--current liquid water specific humidity [kg/kg]
     343REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliqi         !--current ice water specific humidity [kg/kg]
     344REAL,    INTENT(INOUT), DIMENSION(klon) :: zcond          !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
     345REAL,    INTENT(IN),    DIMENSION(klon) :: zfice          !--ice fraction AFTER CONDENSATION [-]
     346REAL,    INTENT(IN),    DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
     347
     348REAL,    INTENT(IN),    DIMENSION(klon) :: rneb           !--cloud fraction [-]
     349REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
     350REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipcld  !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
     351REAL,    INTENT(INOUT), DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
     352REAL,    INTENT(INOUT), DIMENSION(klon) :: tot_zneb       !--total cloud cover above the current mesh [-]
     353REAL,    INTENT(INOUT), DIMENSION(klon) :: zrho_up        !--air density IN THE LAYER ABOVE [kg/m3]
     354REAL,    INTENT(INOUT), DIMENSION(klon) :: zvelo_up       !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]
     355
     356REAL,    INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean [kg/s/m2]
     357REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky [kg/s/m2]
     358REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
     359REAL,    INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean [kg/s/m2]
     360REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
     361REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
     362
     363REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
     364REAL,    INTENT(OUT),   DIMENSION(klon) :: zradoice       !--condensed ice water used in the radiation scheme [kg/kg]
     365REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
     366REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
     367
     368
     369! Local variables for precip fraction update
     370REAL :: smallestreal
     371REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
     372REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
     373REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
     374REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
     375
     376! Local variables for autoconversion
     377REAL :: zct, zcl, zexpo, ffallv
     378REAL :: zchau, zfroi
     379REAL :: zrain, zsnow, zprecip
     380REAL :: effective_zneb
     381REAL, DIMENSION(klon) :: zrho, zvelo
     382REAL, DIMENSION(klon) :: zdz, iwc
     383
     384! Local variables for Bergeron process
     385REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
     386REAL, DIMENSION(klon) :: zqpreci, zqprecl
     387
     388! Local variables for calculation of quantity of condensates seen by the radiative scheme
     389REAL, DIMENSION(klon) :: zradoliq
     390REAL, DIMENSION(klon) :: ziflprev
     391
     392! Misc
     393INTEGER :: i, n
     394
     395! Initialisation
     396smallestreal=1.e-9
     397zqprecl(:) = 0.
     398zqpreci(:) = 0.
     399ziflprev(:) = 0.
     400
     401
     402IF (iflag_evap_prec .GE. 4) THEN
     403
     404  !Partitionning between precipitation coming from clouds and that coming from CS
     405
     406  !0) Calculate tot_zneb, total cloud fraction above the cloud
     407  !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
     408 
     409  DO i=1, klon
     410    tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i),zneb(i))) &
     411            /(1.-min(zneb(i),1.-smallestreal))
     412    d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
     413    tot_zneb(i) = tot_znebn(i)
     414
     415
     416    !1) Cloudy to clear air
     417    d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i),znebprecipcld(i))
     418    IF (znebprecipcld(i) .GT. 0.) THEN
     419      d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
     420      d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
     421    ELSE
     422      d_zrfl_cld_clr(i) = 0.
     423      d_zifl_cld_clr(i) = 0.
     424    ENDIF
     425
     426    !2) Clear to cloudy air
     427    d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i)))
     428    IF (znebprecipclr(i) .GT. 0) THEN
     429      d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
     430      d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
     431    ELSE
     432      d_zrfl_clr_cld(i) = 0.
     433      d_zifl_clr_cld(i) = 0.
     434    ENDIF
     435
     436    !Update variables
     437    znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
     438    znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
     439    zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
     440    ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
     441    zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
     442    ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
     443
     444    ! c_iso  do the same thing for isotopes precip
     445  ENDDO
     446ENDIF
     447
     448
     449! Autoconversion
     450! -------------------------------------------------------------------------------
     451
     452! Initialisation
     453DO i = 1, klon
     454  zrho(i)  = pplay(i) / zt(i) / RD
     455  zdz(i)   = (paprsdn(i)-paprsup(i)) / (zrho(i)*RG)
     456  iwc(i)   = 0.
     457  zneb(i)  = MAX(rneb(i),seuil_neb)
     458
     459  ! variables for calculation of quantity of condensates seen by the radiative scheme
     460  ! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow
     461  ! is TRUE
     462  zradocond(i) = zoliq(i)/REAL(niter_lscp+1)
     463  zradoliq(i)  = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
     464  zradoice(i)  = zoliq(i)*zfice(i)/REAL(niter_lscp+1)
     465ENDDO
     466
     467   
     468DO n = 1, niter_lscp
     469
     470  DO i=1,klon
     471    IF (rneb(i).GT.0.0) THEN
     472      iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
     473    ENDIF
     474  ENDDO
     475
     476  CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)
     477
     478  DO i = 1, klon
     479
     480    IF (rneb(i).GT.0.0) THEN
     481
     482      ! Initialization of zrain, zsnow and zprecip:
     483      zrain=0.
     484      zsnow=0.
     485      zprecip=0.
     486      ! c_iso same init for isotopes. Externalisation?
     487
     488      IF (zneb(i).EQ.seuil_neb) THEN
     489        zprecip = 0.0
     490        zsnow = 0.0
     491        zrain= 0.0
     492      ELSE
     493
     494        IF (ptconv(i)) THEN ! if convective point
     495          zcl=cld_lc_con
     496          zct=1./cld_tau_con
     497          zexpo=cld_expo_con
     498          ffallv=ffallv_con
     499        ELSE
     500          zcl=cld_lc_lsc
     501          zct=1./cld_tau_lsc
     502          zexpo=cld_expo_lsc
     503          ffallv=ffallv_lsc
     504        ENDIF
     505
     506
     507        ! if vertical heterogeneity is taken into account, we use
     508        ! the "true" volume fraction instead of a modified
     509        ! surface fraction (which is larger and artificially
     510        ! reduces the in-cloud water).
     511        IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
     512          effective_zneb=ctot_vol(i)
     513        ELSE
     514          effective_zneb=zneb(i)
     515        ENDIF
     516
     517
     518        ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
     519        ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
     520        !.........................................................
     521        IF (iflag_autoconversion .EQ. 2) THEN
     522        ! two-steps resolution with niter_lscp=1 sufficient
     523        ! we first treat the second term (with exponential) in an explicit way
     524        ! and then treat the first term (-q/tau) in an exact way
     525          zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
     526            *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
     527        ELSE
     528        ! old explicit resolution with subtimesteps
     529          zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
     530            *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
     531        ENDIF
     532
     533
     534        ! Ice water quantity to remove (Zender & Kiehl, 1997)
     535        ! dqice/dt=1/rho*d(rho*wice*qice)/dz
     536        !....................................
     537        IF (iflag_autoconversion .EQ. 2) THEN
     538        ! exact resolution, niter_lscp=1 is sufficient but works only
     539        ! with iflag_vice=0
     540          IF (zoliqi(i) .GT. 0.) THEN
     541            zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
     542              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
     543          ELSE
     544            zfroi=0.
     545          ENDIF
     546        ELSE
     547        ! old explicit resolution with subtimesteps
     548          zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*zvelo(i)
     549        ENDIF
     550
     551        zrain   = MIN(MAX(zchau,0.0),zoliql(i))
     552        zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
     553        zprecip = MAX(zrain + zsnow,0.0)
     554
     555      ENDIF
     556
     557
     558      IF (iflag_autoconversion .GE. 1) THEN
     559        ! debugged version with phase conservation through the autoconversion process
     560        zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
     561        zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
     562        zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
     563      ELSE
     564        ! bugged version with phase resetting
     565        zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
     566        zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
     567        zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
     568      ENDIF
     569
     570      ! c_iso: call isotope_conversion (for readibility) or duplicate
     571
     572      ! variables for calculation of quantity of condensates seen by the radiative scheme
     573      zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp+1)
     574      zradoliq(i)  = zradoliq(i) + zoliql(i)/REAL(niter_lscp+1)
     575      zradoice(i)  = zradoice(i) + zoliqi(i)/REAL(niter_lscp+1)
     576
     577      !--Diagnostics
     578      dqrauto(i) = dqrauto(i) + zrain / dtime
     579      dqsauto(i) = dqsauto(i) + zsnow / dtime
     580
     581    ENDIF ! rneb >0
     582
     583  ENDDO  ! i = 1,klon
     584
     585ENDDO      ! n = 1,niter
     586
     587! Precipitation flux calculation
     588
     589DO i = 1, klon
     590       
     591  IF (iflag_evap_prec.GE.4) THEN
     592    ziflprev(i)=ziflcld(i)
     593  ELSE
     594    ziflprev(i)=zifl(i)*zneb(i)
     595  ENDIF
     596
     597  IF (rneb(i) .GT. 0.0) THEN
     598
     599    ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
     600    ! If T<0C, liquid precip are converted into ice, which leads to an increase in
     601    ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
     602    ! taken into account through a linearization of the equations and by approximating
     603    ! the liquid precipitation process with a "threshold" process. We assume that
     604    ! condensates are not modified during this operation. Liquid precipitation is
     605    ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
     606    ! Water vapor increases as well
     607    ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
     608
     609    IF (ok_bug_phase_lscp) THEN
     610      zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
     611      zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
     612    ELSE
     613      zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i)
     614      zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i)
     615    ENDIF
     616    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     617    coef1 = rneb(i)*RLSTT/zcp*zdqsdT_raw(i)
     618    ! Computation of DT if all the liquid precip freezes
     619    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
     620    ! T should not exceed the freezing point
     621    ! that is Delta > RTT-zt(i)
     622    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
     623    zt(i)  = zt(i) + DeltaT
     624    ! water vaporization due to temp. increase
     625    Deltaq = rneb(i)*zdqsdT_raw(i)*DeltaT
     626    ! we add this vaporized water to the total vapor and we remove it from the precip
     627    zq(i) = zq(i) +  Deltaq
     628    ! The three "max" lines herebelow protect from rounding errors
     629    zcond(i) = max( zcond(i)- Deltaq, 0. )
     630    ! liquid precipitation converted to ice precip
     631    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
     632    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
     633    ! iced water budget
     634    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
     635
     636    ! c_iso : mv here condensation of isotopes + redispatchage en precip
     637
     638    IF (iflag_evap_prec.GE.4) THEN
     639      zrflcld(i) = zrflcld(i)+zqprecl(i) &
     640        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     641      ziflcld(i) = ziflcld(i)+ zqpreci(i) &
     642        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     643      znebprecipcld(i) = rneb(i)
     644      zrfl(i) = zrflcld(i) + zrflclr(i)
     645      zifl(i) = ziflcld(i) + ziflclr(i)
     646    ELSE
     647      zrfl(i) = zrfl(i)+ zqprecl(i)        &
     648        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     649      zifl(i) = zifl(i)+ zqpreci(i)        &
     650        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     651    ENDIF
     652    ! c_iso : same for isotopes
     653
     654  ENDIF ! rneb>0
     655
     656ENDDO
     657
     658! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
     659! if iflag_evap_prec>=4
     660IF (iflag_evap_prec.GE.4) THEN
     661
     662  DO i=1,klon
     663
     664    IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
     665      znebprecipclr(i) = min(znebprecipclr(i),max( &
     666        zrflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), &
     667        ziflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
     668    ELSE
     669      znebprecipclr(i)=0.0
     670    ENDIF
     671
     672    IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
     673      znebprecipcld(i) = min(znebprecipcld(i), max( &
     674        zrflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), &
     675        ziflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
     676    ELSE
     677      znebprecipcld(i)=0.0
     678    ENDIF
     679  ENDDO
     680
     681ENDIF
     682
     683! we recalculate zradoice to account for contributions from the precipitation flux
     684! if ok_radocond_snow is true
     685IF ( ok_radocond_snow ) THEN
     686  IF ( .NOT. iftop ) THEN
     687    ! for the solid phase (crystals + snowflakes)
     688    ! we recalculate zradoice by summing
     689    ! the ice content calculated in the mesh
     690    ! + the contribution of the non-evaporated snowfall
     691    ! from the overlying layer
     692    DO i=1,klon
     693      IF (ziflprev(i) .NE. 0.0) THEN
     694        zradoice(i)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho_up(i)/zvelo_up(i)
     695      ELSE
     696        zradoice(i)=zoliqi(i)+zqpreci(i)
     697      ENDIF
     698      zradocond(i)=zradoliq(i)+zradoice(i)
     699    ENDDO
     700  ENDIF
     701  ! keep in memory air density and ice fallspeed velocity
     702  zrho_up(:) = zrho(:)
     703  zvelo_up(:) = zvelo(:)
     704ENDIF
     705
     706END SUBROUTINE histprecip_postcld
    307707
    308708
Note: See TracChangeset for help on using the changeset viewer.