5 SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6 d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, &
7 pfrac_impa, pfrac_nucl, pfrac_1nucl, &
8 frac_impa, frac_nucl,
beta, &
9 prfl, psfl, rhcl, zqta, fraca, &
10 ztv, zpspsk, ztla, zthl, iflag_cld_th, &
54 REAL ztfondue, qsl, qsi
56 logical lognormale(
klon)
82 INTEGER iflag_ice_thermo
102 REAL zqs(
klon), zdqs(
klon), zdelta, zcor, zcvm5
104 LOGICAL convergence(
klon)
107 INTEGER n_i(
klon), iter
110 REAL zrfl(
klon), zrfln(
klon), zqev, zqevt
111 REAL zifl(
klon), zifln(
klon), zqev0,zqevi, zqevti
118 INTEGER exposant_glace_old
121 REAL zchau ,zfroi ,zfice(
klon),zneb(
klon)
122 REAL zmelt, zpluie, zice, zcondold
144 REAL zprec_cond(
klon)
149 REAL zmair, zcpair, zcpeau
151 REAL zlh_solid(
klon), zm_solid
162 fallvc(zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
163 fallvs(zzz) = 3.29/2.0 * ((zzz)**0.16) *
ffallv_lsc
165 DATA appel1er /.
true./
169 ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
175 WRITE(
lunout,*)
'fisrtilp, ninter:', ninter
176 WRITE(
lunout,*)
'fisrtilp, evap_prec:', evap_prec
177 WRITE(
lunout,*)
'fisrtilp, cpartiel:', cpartiel
178 IF (abs(dtime/
REAL(ninter)-360.0).GT.0.001) then
179 WRITE(
lunout,*)
'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
180 WRITE(
lunout,*)
'Je prefere un sous-intervalle de 6 minutes'
215 t_glace_min_old = rtt - 15.0
219 exposant_glace_old = 2
222 exposant_glace_old = 6
300 zmair=(paprs(i,k)-paprs(i,k+1))/
rg
301 zcpair=rcpd*(1.0+rvtmp2*zq(i))
303 zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
304 + zmair*zcpair*zt(i) ) &
305 / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
323 IF (zrfl(i)+zifl(i).GT.0.)
THEN
326 zdelta=max(0.,sign(1.,rtt-zt(i)))
327 zqs(i)= r2es*foeew(zt(i),zdelta)/pplay(i,k)
328 zqs(i)=min(0.5,zqs(i))
329 zcor=1./(1.-retv*zqs(i))
332 IF (zt(i) .LT. t_coup)
THEN
333 zqs(i) = qsats(zt(i)) / pplay(i,k)
335 zqs(i) = qsatl(zt(i)) / pplay(i,k)
341 IF (.NOT. ice_thermo)
THEN
345 IF (zrfl(i)+zifl(i).GT.0.)
THEN
347 zqev = max(0.0, (zqs(i)-zq(i))*zneb(i) )
348 zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * sqrt(zrfl(i)) &
349 * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/
rg
350 zqevt = max(0.0,min(zqevt,zrfl(i))) &
351 *
rg*dtime/(paprs(i,k)-paprs(i,k+1))
352 zqev = min(zqev, zqevt)
353 zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
361 IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
363 zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
364 * (
rg/(paprs(i,k)-paprs(i,k+1)))*dtime
365 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
366 * (
rg/(paprs(i,k)-paprs(i,k+1)))*dtime &
367 * rlvtt/rcpd/(1.0+rvtmp2*zq(i))
378 IF (zrfl(i)+zifl(i).GT.0.)
THEN
386 zqev0 = max(0.0, (zqs(i)-zq(i))*zneb(i) )
395 qsl= r2es*foeew(zt(i),0.)/pplay(i,k)
397 zcor= 1./(1.-retv*qsl)
400 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*sqrt(zrfl(i)) &
401 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/
rg
402 zqevt = max(0.0,min(zqevt,zrfl(i))) &
403 *
rg*dtime/(paprs(i,k)-paprs(i,k+1))
405 qsi= r2es*foeew(zt(i),1.)/pplay(i,k)
407 zcor= 1./(1.-retv*qsi)
410 zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*sqrt(zifl(i)) &
411 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/
rg
412 zqevti = max(0.0,min(zqevti,zifl(i))) &
413 *
rg*dtime/(paprs(i,k)-paprs(i,k+1))
420 IF (zqevt+zqevti.GT.zqev0)
THEN
421 zqev=zqev0*zqevt/(zqevt+zqevti)
422 zqevi=zqev0*zqevti/(zqevt+zqevti)
425 IF (zqevt+zqevti.GT.0.)
THEN
426 zqev=min(zqev0*zqevt/(zqevt+zqevti),zqevt)
427 zqevi=min(zqev0*zqevti/(zqevt+zqevti),zqevti)
434 zrfln(i) = max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
436 zifln(i) = max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
445 zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
446 * (
rg/(paprs(i,k)-paprs(i,k+1)))*dtime
447 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
448 * (
rg/(paprs(i,k)-paprs(i,k+1)))*dtime &
449 * rlvtt/rcpd/(1.0+rvtmp2*zq(i)) &
450 + (zifln(i)-zifl(i)) &
451 * (
rg/(paprs(i,k)-paprs(i,k+1)))*dtime &
452 * rlstt/rcpd/(1.0+rvtmp2*zq(i))
458 zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
459 zmelt = min(max(zmelt,0.),1.)
460 zrfl(i)=zrfl(i)+zmelt*zifl(i)
461 zifl(i)=zifl(i)*(1.-zmelt)
463 zt(i)=zt(i)-zifl(i)*zmelt*(
rg*dtime)/(paprs(i,k)-paprs(i,k+1)) &
464 *rlmlt/rcpd/(1.0+rvtmp2*zq(i))
481 zdelta = max(0.,sign(1.,rtt-zt(i)))
482 zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
483 zcvm5 = zcvm5 /rcpd/(1.0+rvtmp2*zq(i))
484 zqs(i) = r2es*foeew(zt(i),zdelta)/pplay(i,k)
485 zqs(i) = min(0.5,zqs(i))
486 zcor = 1./(1.-retv*zqs(i))
488 zdqs(i) = foede(zt(i),zdelta,zcvm5,zqs(i),zcor)
492 IF (zt(i).LT.t_coup)
THEN
493 zqs(i) = qsats(zt(i))/pplay(i,k)
494 zdqs(i) = dqsats(zt(i),zqs(i))
496 zqs(i) = qsatl(zt(i))/pplay(i,k)
497 zdqs(i) = dqsatl(zt(i),zqs(i))
526 if (iflag_pdf.eq.0)
then
529 zdelq = min(ratqs(i,k),0.99) * zq(i)
530 rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
531 zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
538 if(zq(i).lt.1.e-15)
then
544 if (iflag_cld_th>=5)
then
548 qcloud,ctot,zpspsk,paprs,ztla,zthl, &
558 if (iflag_cld_th <= 4)
then
560 elseif (iflag_cld_th >= 6)
then
562 lognormale = fraca(:,k) < 1e-10
580 if (iflag_fisrtilp_qsat.ge.0)
then
581 do iter=1,iflag_fisrtilp_qsat+1
585 convergence(i)=abs(dt(i)).gt.ddt0
586 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i))
then
587 tbef(i)=tbef(i)+dt(i)
588 if (.not.ice_thermo)
then
589 zdelta = max(0.,sign(1.,rtt-tbef(i)))
592 zdelta = max(0.,sign(1.,t_glace_min_old-tbef(i)))
597 zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
598 zcvm5 = zcvm5 /rcpd/(1.0+rvtmp2*zq(i))
599 zqs(i) = r2es*foeew(tbef(i),zdelta)/pplay(i,k)
600 zqs(i) = min(0.5,zqs(i))
601 zcor = 1./(1.-retv*zqs(i))
603 zdqs(i) = foede(tbef(i),zdelta,zcvm5,zqs(i),zcor)
604 zpdf_sig(i)=ratqs(i,k)*zq(i)
605 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
606 zpdf_delta(i)=log(zq(i)/zqs(i))
607 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
608 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
609 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
610 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
611 zpdf_e1(i)=1.-erf(zpdf_e1(i))
612 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
613 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
614 zpdf_e2(i)=1.-erf(zpdf_e2(i))
616 if (zpdf_e1(i).lt.1.e-10)
then
620 rneb(i,k)=0.5*zpdf_e1(i)
621 zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
627 if (.not. ice_thermo)
then
630 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i))
then
632 qlbef(i)=max(0.,zqn(i)-zqs(i))
633 num=-tbef(i)+zt(i)+rneb(i,k)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))*qlbef(i)
634 denom=1.+rneb(i,k)*zdqs(i)
654 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i))
then
657 zfice(i) = 1.0 - (tbef(i)-t_glace_min_old) / (rtt-t_glace_min_old)
658 zfice(i) = min(max(zfice(i),0.0),1.0)
659 zfice(i) = zfice(i)**exposant_glace_old
660 dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - rtt)
667 if ((zfice(i).eq.0).or.(zfice(i).eq.1))
then
671 if (zfice(i).lt.1)
then
677 qlbef(i)=max(0.,zqn(i)-zqs(i))
678 num=-tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*rlvtt+zfice(i)*rlstt)/rcpd/(1.0+rvtmp2*zq(i))*qlbef(i)
679 denom=1.+rneb(i,k)*((1-zfice(i))*rlvtt+zfice(i)*rlstt)/cste*zdqs(i) &
680 -(rlstt-rlvtt)/rcpd/(1.0+rvtmp2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
704 IF (rneb(i,k) .LE. 0.0)
THEN
708 rhcl(i,k)=zq(i)/zqs(i)
709 ELSE IF (rneb(i,k) .GE. 1.0)
THEN
712 zcond(i) = max(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
715 zcond(i) = max(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
716 rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
724 IF (rneb(i,k) .LE. 0.0)
THEN
728 rhcl(i,k)=zq(i)/zqs(i)
729 ELSE IF (rneb(i,k) .GE. 1.0)
THEN
732 zcond(i) = max(0.0,zqn(i)-zqs(i))
735 zcond(i) = max(0.0,zqn(i)-zqs(i))*rneb(i,k)
736 rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
761 IF (zq(i).GT.zqs(i))
THEN
766 zcond(i) = max(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
771 zq(i) = zq(i) - zcond(i)
775 IF (.NOT. ice_thermo)
THEN
776 if (iflag_fisrtilp_qsat.lt.1)
then
778 zt(i) = zt(i) + zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*zq(i))
780 else if (iflag_fisrtilp_qsat.gt.0)
then
782 zt(i) = zt(i) + zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*(zq(i)+zcond(i)))
789 if (iflag_fisrtilp_qsat.lt.1)
then
795 zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (rtt-t_glace_min_old)
796 zfice(i) = min(max(zfice(i),0.0),1.0)
797 zfice(i) = zfice(i)**exposant_glace_old
799 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*zq(i)) &
800 +zfice(i)*zcond(i) * rlstt/rcpd/(1.0+rvtmp2*zq(i))
808 zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (rtt-t_glace_min_old)
809 zfice(i) = min(max(zfice(i),0.0),1.0)
810 zfice(i) = zfice(i)**exposant_glace_old
812 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * rlvtt/rcpd/(1.0+rvtmp2*(zq(i)+zcond(i))) &
813 +zfice(i)*zcond(i) * rlstt/rcpd/(1.0+rvtmp2*(zq(i)+zcond(i)))
823 IF (rneb(i,k).GT.0.0)
THEN
825 zrho(i) = pplay(i,k) / zt(i) / rd
826 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*
rg)
830 IF (.NOT. ice_thermo)
THEN
833 IF (rneb(i,k).GT.0.0)
THEN
834 zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
835 zfice(i) = min(max(zfice(i),0.0),1.0)
836 zfice(i) = zfice(i)**exposant_glace_old
853 IF (rneb(i,k).GT.0.0)
THEN
854 zneb(i) = max(rneb(i,k), seuil_neb)
858 radliq(i,k) = zoliq(i)/
REAL(ninter+1)
864 IF (rneb(i,k).GT.0.0)
THEN
865 zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
869 IF (zneb(i).EQ.seuil_neb)
THEN
874 if (ptconv(i,k))
then
877 zfroi = dtime/
REAL(ninter)/zdz(i)*zoliq(i) &
878 *fallvc(zrhol(i)) * zfice(i)
882 zfroi = dtime/
REAL(ninter)/zdz(i)*zoliq(i) &
883 *fallvs(zrhol(i)) * zfice(i)
885 zchau = zct *dtime/
REAL(ninter) * zoliq(i) &
886 *(1.0-exp(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i))
888 IF (.NOT. ice_thermo)
THEN
891 zpluie = min(max(zchau,0.0),zoliq(i)*(1.-zfice(i)))
892 zice = min(max(zfroi,0.0),zoliq(i)*zfice(i))
896 ztot = max(ztot ,0.0)
898 ztot = min(ztot,zoliq(i))
902 zoliqp(i) = max(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0)
903 zoliqi(i) = max(zoliq(i)*zfice(i)-1.*zice , 0.0)
904 zoliq(i) = max(zoliq(i)-ztot , 0.0)
906 radliq(i,k) = radliq(i,k) + zoliq(i)/
REAL(ninter+1)
911 IF (.NOT. ice_thermo)
THEN
913 IF (rneb(i,k).GT.0.0)
THEN
915 zrfl(i) = zrfl(i)+ max(zcond(i)-zoliq(i),0.0) &
916 * (paprs(i,k)-paprs(i,k+1))/(
rg*dtime)
921 IF (rneb(i,k).GT.0.0)
THEN
923 if (.not.ice_thermo)
then
927 d_ql(i,k) = (1-zfice(i))*zoliq(i)
928 d_qi(i,k) = zfice(i)*zoliq(i)
931 zrfl(i) = zrfl(i)+ max(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
932 *(paprs(i,k)-paprs(i,k+1))/(
rg*dtime)
933 zifl(i) = zifl(i)+ max(zcond(i)*zfice(i)-zoliqi(i),0.0) &
934 *(paprs(i,k)-paprs(i,k+1))/(
rg*dtime)
959 IF (.NOT. ice_thermo)
THEN
961 IF (zt(i).LT.rtt)
THEN
986 d_q(i,k) = zq(i) - q(i,k)
987 d_t(i,k) = zt(i) - t(i,k)
994 if(zcond(i).gt.zoliq(i)+1.e-10)
then
995 beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
999 zprec_cond(i) = max(zcond(i)-zoliq(i),0.0) &
1000 * (paprs(i,k)-paprs(i,k+1))/
rg
1001 IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.)
THEN
1004 if (t(i,k) .GE. t_glace_min_old)
THEN
1005 zalpha_tr = a_tr_sca(3)
1007 zalpha_tr = a_tr_sca(4)
1011 zalpha_tr = a_tr_sca(3)
1013 zalpha_tr = a_tr_sca(4)
1016 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
1017 pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1018 frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1021 zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
1022 pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1030 IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.)
THEN
1032 if (t(i,kk) .GE. t_glace_min_old)
THEN
1033 zalpha_tr = a_tr_sca(1)
1035 zalpha_tr = a_tr_sca(2)
1039 zalpha_tr = a_tr_sca(1)
1041 zalpha_tr = a_tr_sca(2)
1044 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
1045 pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1046 frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1060 IF (.NOT.ice_thermo)
THEN
1062 IF ((t(i,1)+d_t(i,1)) .LT. rtt)
THEN
1065 snow(i) = zrfl(i)+zifl(i)
1067 zlh_solid(i) = rlstt-rlvtt
1085 IF (.not.ice_thermo)
THEN
1088 zcpair=rcpd*(1.0+rvtmp2*(q(i,k)+d_q(i,k)))
1089 zmair=(paprs(i,k)-paprs(i,k+1))/
rg
1090 zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1091 d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
1097 if (ncoreczq>0)
then
1098 WRITE(
lunout,*)
'WARNING : ZQ dans fisrtilp ',ncoreczq,
' val < 1.e-15.'
subroutine fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, beta, prfl, psfl, rhcl, zqta, fraca, ztv, zpspsk, ztla, zthl, iflag_cld_th, iflag_ice_thermo)
subroutine icefrac_lsc(np, temp, sig, icefrac)
!$Id t_glace_min REAL exposant_glace REAL rei_max REAL coefw_cld_cv REAL tmax_fonte_cv INTEGER iflag_cld_cv common nuagecom && t_glace_min
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real beta
!$Id t_glace_min REAL exposant_glace REAL rei_max REAL coefw_cld_cv REAL tmax_fonte_cv INTEGER iflag_cld_cv common nuagecom exposant_glace
!$Id t_glace_min REAL exposant_glace REAL rei_max REAL coefw_cld_cv REAL tmax_fonte_cv INTEGER iflag_t_glace
!$Id cld_lc_con REAL cld_tau_con REAL ffallv_lsc
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
!$Id cld_lc_con REAL cld_tau_lsc
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
subroutine cloudth(ngrid, klev, ind2, ztv, po, zqta, fraca, qcloud, ctot, zpspsk, paprs, ztla, zthl, ratqs, zqs, t)
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout