c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c call sfcddep(deltim,sm(i,j),snd,z0,ustar(i,j),u10(i,j),v10(i,j) c + ,NTHP,NTHS,NTHV,PLM,TLM,depd,depw c + ,LM,rrcol_t(:,i,j),zint(i,j,:), scol_t(:,i,j)) subroutine sfcddep (deltim,sm,snd,z0,ustar,u10,v10 + ,NTHP,NTHS,NTHV,PLM,TLM,depd,depw . ,LMHK,RRCOL,DZINT,S) ! ****************************************************************** ! * * ! * surface layer * ! * responsible person: z.janjic * ! * * ! ****************************************************************** !----------------------------------------------------------------------- PARAMETER (KPSZ=8,MSIB=13,LZOB=9,NPS=KPSZ) dimension slm(KPSZ) c INCLUDE "CTLBLK.comm" !*********************************************************************** P A R A M E T E R &(SMALL=0.35, GLKBS=30.0,GLKBR=10.0,GRRS=GLKBR/GLKBS &,G=9.806,VISC=1.5E-5, TVISC=2.1E-5, QVISC=2.1E-5 &,CZIV=SMALL*GLKBS,USTFC=0.018/G &,DVISF1=5.00E-9 ! dust diffusn (smooth) &,DVISF2=3.00E-8 ! dust diffusn (transit) &,DVISF3=7.00E-8 ! dust diffusn (rough) &,SVISC1=DVISF1*2.1E-5 &,SVISC2=DVISF2*2.1E-5 &,SVISC3=DVISF3*2.1E-5 &,RSVISC1=1./SVISC1 &,RSVISC2=1./SVISC2 &,RSVISC3=1./SVISC3 &,DVISF=0.1 ! dust diffus. factor &,SVISC=DVISF*2.1E-5 ! partcle 'mlclr' vscsty &,RVISC=1./VISC,RTVISC=1./TVISC,RQVISC=1./QVISC c &,AA=0.00177 ! avrg. obstc. radius &,AAPAZI=0.00177 ! avrg. obstc. radius &,RSVISC=1./SVISC &,SQPR=0.84,SQSC=0.84,ZQRZT=SQSC/SQPR &,SQSCS=3.2,ZQRZS=SQSCS/SQPR &,USTRS=0.225,USTCS=0.7,USTARTR=0.15 &,FZU1=CZIV*VISC,FZT1=RVISC *TVISC*SQPR,FZQ1=RTVISC*QVISC*ZQRZT &,FZT2=CZIV*GRRS*TVISC*SQPR,FZQ2=RTVISC*QVISC*ZQRZT &,FZS1=RTVISC*SVISC*ZQRZS &,FZS2=RTVISC*SVISC*ZQRZS &,FZS11=RTVISC*SVISC1*ZQRZS &,FZS12=RTVISC*SVISC2*ZQRZS &,FZS13=RTVISC*SVISC3*ZQRZS &,FZS21=FZS11 &,FZS22=FZS12 &,FZS23=FZS13) DIMENSION PRADI(KPSZ),PDENS(KPSZ),GAMA(KPSZ),USTRTR(KPSZ) cbojan DIMENSION GPSS(KPSZ) !zender DIMENSION PSSMXX(KPSZ),PVTERM(KPSZ), WRAT(KPSZ),AA(KPSZ) DIMENSION WPR(LZOB),BCLY(LZOB),BSLT(LZOB),BSND(LZOB) DIMENSION CD(MSIB),CV(MSIB),ASM(MSIB),RSQCD0V(MSIB) DIMENSION BETA(KPSZ,LZOB),DELTA(KPSZ) DIMENSION VTERM(KPSZ) D I M E N S I O N & STKN(NPS),STKN2(NPS),STKNA(NPS),FSTK(NPS),GG(NPS),FBO(NPS) &,CEFF(NPS),VDIL(NPS),VDDP(NPS) &,SCHD(NPS) DIMENSION RR(NPS), RSM(NPS) DIMENSION RRCOL(LMHK),DZINT(LMHK+1),S(LMHK) dimension SSSS(200),sold(200),fprep(200) DIMENSION acw(NPS), acw1(NPS) DATA CD /.50,.50,.50,.40,.40,.40,.30,.35,.35,.35,0.,30.,0./ DATA CV /.16,.16,.16,.12,.12,.12,.08,.09,.09,.09,0.,08.,0./ DATA RSQCD0V/2.0,2.0,2.0,2.5,2.5,2.5,3.5,3.0,3.0,3.0,0.,3.5,0./ DATA ASM /.00,.00,.00,.00,.00,.00,1.4E-5,3.2E-5,3.2E-5,3.2E-5 &, .00,1.4E-5,.00/ !cboj DATA PRADI /0.15,0.25,0.47,0.80,1.36,2.29,3.93,7.24/ !microns ! vol rad DATA PRADI /0.175,0.225,0.375,0.725,1.5,3.0,6.0,9.0/ !microns ! vol rad DATA PDENS /2.50,2.50,2.50,2.50,2.65,2.65,2.65,2.65/ DATA ADENS /0.00122/ !part. density (kg/m^3) DATA GAMA /0.08,0.08,0.08,0.08,1.,1.,1.,1./ DATA PSSMXX /.0008,.0008,.0008,.0008 &, .00008,.00008,.00008,.00008/ DATA WRAT /5.E5,5.E5,5.E5,5.E5,5.E5,5.E5,5.E5,5.E5/ C BAGNOLD'S-TYPE FOR U* THRESHOLD SEE WHITE(1979) for A(PSZ) DATA AA /1.,0.9,0.8,0.8,0.7,0.6,0.5,0.4/ DATA WPR /2.5,6.8,11.5,2.5,10.,6.8,3.5,0.0,0.0/ DATA BCLY /.12,.34,.45,.12,.40,.34,.22,.00,.00/ DATA BSLT /.08,.56,.30,.18,.10,.36,.18,.00,.00/ DATA BSND /.80,.10,.25,.70,.50,.30,.60,.00,.00/ C----------------------------------------------------------------------- N = NTHP KK = N K = N LL = NTHS ALPHA = snd IVGTYP=NTHV slm(NTHP) = S(LMHK) c if ( SND.lt.1) return VTERM(K)=2*PDENS(K)*(PRADI(K))**2/(9*1.5E-5*9.8)/1E+6/10 c DO N=1,KPSZ USTRTR(N)=G*(2.*1E-6*PRADI(N))*(PDENS(N)-ADENS)/ADENS USTRTR(N)=AA(N)*SQRT(USTRTR(N)) c ENDDO C=========== C DEPOS.f c wetdepostion! c goto 21212 DO N=NTHP,NTHP ACW(N) = 0 DO L=1,LMHK DEZ=(DZINT(L)-DZINT(L+1)) IF (S(L).LT.0.) S(L)=0.0 fprep(l) = rrcol(l)/1000. ! kg/m2s to m/s if(fprep(l).gt.0.) then SOLD(L)=S(L) ! flx=SOLD(L)*WRAT(N)*FPREP(L) !! flux C*Wr*Vel ! dec = min(flx*deltim,SOLD(l)) ! DEC=SLM(N)*VDDP(N)/(DEZ/DELTIM) DEC=SOLD(L)*WRAT(N)*FPREP(L)/(DEZ/DELTIM) DEC=min(DEC,SOLD(L)) S(L) = SOLD(l)-dec ! ACW(N)=ACW(n)+dec ! deposition ACW(N)=ACW(n)+dec*DEZ ! print *,l,DEZ,n,acw(n),dec*1.e+9 endif enddo depw = ACW(n) ! DEPW(N)=DEPW(N)+ACW(N) ! DEPT(N)=DEPT(N)+ACW(N) enddo ! ACW(N)=SOLD(L)*WRAT(N)*FPREP(L) !! flux C*Wr*Vel ! DO N=NTHP,NTHP ! DO L=1,LMHK c print *,L,dzint(L) ! IF (S(L).LT.0.) S(L)=0.0 ! fprep(l) = rrcol(l)/1000. ! kg/m2s to m/s ! SOLD(L)=S(L) C ----------------- WET DEPOSITION ------------------------------------ ! IF (FPREP(L).GT.0.0) THEN ! IF(L.EQ.1)THEN ! ACW1(N)=0. ! ELSE ! ACW1(N)= ! & SOLD(L-1)*WRAT(N)*FPREP(L-1) ! & /(DZINT(L-1)-DZINT(L)) ! IF(ACW1(N).GT.SOLD(L-1)) ACW1(N)=SOLD(L-1) ! ENDIF ! dez = (DZINT(L-1)-DZINT(L)) ! alm=WRAT(N)*fprep(l)/dez ! Fc=alm*C*dz ! ACW(N)=SOLD(L)*WRAT(N)*FPREP(L) !! flux C*Wr*Vel ! dezl = DZINT(L)-DZINT(L+1) c if ( WRAT(N)*FPREP(L).gt.0.7*dezl) then c print *,"WRAT(N)",WRAT(N),FPREP(L),WRAT(N)*FPREP(L) c print *,"dzint ",l,dezl c endif cpazi ! carlos May2008 ! SSSS(N)=S(L)-ACW(N)/(DZINT(L)-DZINT(L+1)) c & +ACW1(N) ! IF(SSSS(N).LT.0)THEN ! SSSS(N)=0. c ACW(N)=(DZINT(I,J,L)-DZINT(I,J,L+1))*(S(I,J,L,N)+ACW1(N)) ! ACW(N)=ACW(N)+ACW1(N)*(DZINT(L)-DZINT(L+1)) ! ENDIF ! S(L)=max(0.,SSSS(N)) ! IF(L.EQ.LMHK) THEN ! IF(ACW(N).LT.0.)ACW(N)=0. ! DEPW=ACW(N) c DEPW(N)=DEPW(N)+ACW(N) c DEPT(N)=DEPT(N)+ACW(N) ! ENDIF ! ENDIF ! ENDDO !level ! ENDDO ! endof PARTICLEM cpazi 21212 continue slm(NTHP) = S(LMHK) C=========== C##### c IF(MOD(NTSD-NPHS/2,NPHS).EQ.0) T H E N c IF(MOD(NTSD-NPHS/2,NPHS).EQ.-100) T H E N C ----------------- DRY DEPOSITION ------------------------------------ C!!!!!!!!!!!!!!!!!!!!!!!!! BCONS=1.38*1.E-23 ! Boltzman constant RC=287.04 ! Gass constant PINO=3.14159265 ! Pi C DO N=NTHP,NTHP C CUNNINGHAN CORRECTION FACTOR (CCNN) FPT=0.065E-9 CCNN=1+FPT/(PRADI(N)*1E-9)*(1.257+0.4*EXP(-1.1*PRADI(N)*1E-9/FPT)) c LMHK=LM c pazi tlm plm APESTSS=PLM ! PDSL(I,J)*AETA(LMHK)+PT SCHD(N)=CCNN*VISC/(BCONS*TLM**2*RC/ &(6.*PINO*VISC*PRADI(N)*1E-9*APESTSS)) ENDDO C UVH=SQRT(U10*U10+V10*V10) UVH=AMAX1(UVH,0.01) CDR=(USTAR*USTAR)/(UVH*UVH) ! drag coeff. CDR=AMIN1(CDR,0.0005) C DO N=NTHP,NTHP C IF(SM.EQ.1. .OR. ! sea surface & IVGTYP.EQ.11 .OR. ! bare soil & IVGTYP.EQ.13 )THEN ! ice surface REYN=Z0*USTAR/VISC ! Re number STKN(N)=VTERM(N)*USTAR*USTAR/(G*VISC) ! Stokes number STKN2(N)=STKN(N)*STKN(N) FSTK(N)=STKN2(N)/(STKN2(N)+400.) C IF(REYN.LE.0.13) THEN ! smooth RSQCD0=13.5 SQCD0=1./RSQCD0 GG(N)=SCHD(N)**(-2./3.)+4.72*FSTK(N) ELSE IF(REYN.GT.0.13.AND.REYN.LE.2.) THEN ! transition RSQCD0=6.432*REYN**(-0.3634) SQCD0=1./RSQCD0 GG(N)=0.6667*REYN**(-0.2)*SCHD(N)**(-0.538*REYN**(-0.105)) & +2.225*REYN**(-0.3634)*FSTK(N) ELSE ! rough RSQCD0=5. SQCD0=1./RSQCD0 GG(N)=0.6849*REYN**(-0.25)/SQRT(SCHD(N))+1.75*FSTK(N) ENDIF VDSL=CDR*UVH*SQRT(CDR)/(SQCD0-SQRT(CDR)) !corrected error! VDIL(N)=SQRT(CDR)*USTAR*GG(N) !corrected error! cc print*,'corrected vdsl, vdil(n) ',vdsl*1000, vdil(n)*1000 VDSL=CDR*UVH*SQCD0/(SQCD0-SQRT(CDR)) VDIL(N)=SQCD0*USTAR*GG(N) cc print*,'OLD vdsl, vdil(n) ',vdsl*1000, vdil(n)*1000 cc print*,'------------------------' VDDP(N)=VDSL*VDIL(N)/(VDSL+VDIL(N)) cc print*,'vdsl= ', vdsl*1000 cc print*,'vdil= ', vdil(n)*1000 cc print*,'===========' ELSE ! vegetation RSQCD0VG=RSQCD0V(IVGTYP) SQCD0VG=1./RSQCD0VG RR(N)=(PRADI(N)*1.E-9)/AAPAZI STKNA(N)=VTERM(N)*USTAR/(G*AAPAZI) ! Stanton number (veg) FBO(N)=EXP(-SQRT(STKNA(N))) ! bluff-of factor RSM(N)=PRADI(N)/10.E-9 C CEFF(N)=CV(IVGTYP)*SCHD(N)**(-1.3) & +0.5*(RR(N)*RR(N))/(AAPAZI*AAPAZI) & +(STKNA(N)/(STKNA(N)+0.6))**3.2 & +ASM(IVGTYP)*RSM(N)*ALOG(1.+RSM(N)) GG(N)=SQRT(CEFF(N)/CD(IVGTYP)) C VDSL=CDR*UVH*SQRT(CDR)/(SQCD0VG-SQRT(CDR)) !corrected error! VDIL(N)=SQRT(CDR)*USTAR*GG(N) !corrected error! cc print*,'corrected vdsl, vdil(n) ',vdsl*1000, vdil(n)*1000 VDSL=CDR*UVH*SQCD0VG/(SQCD0VG-SQRT(CDR)) VDIL(N)=SQCD0VG*USTAR*GG(N) cc print*,'OLD vdsl, vdil(n) ',vdsl*1000, vdil(n)*1000 cc print*,'------------------------' VDDP(N)=VDSL*VDIL(N)*FBO(N)/(VDSL+VDIL(N)*FBO(N)) ENDIF C ENDDO C C!!!!!!!!!!!!!!!!!!!!!!!!! LMH = LMHK DEZ=DZINT(LMH)-DZINT(LMH+1) DO N=NTHP,NTHP DEC=SLM(N)*VDDP(N)/(DEZ/DELTIM) SLM(N)=SLM(N)-DEC IF(SLM(N).LT.0.)SLM(N)=0. c DEPD(N)=DEPD(N)+DEC*DEZ c DEPT=DEPT+DEC*DEZ DEPD=DEC*DEZ c print *,dez,IVGTYP,slm1,slm(n) S(LMHK) = slm(n) ENDDO C----------------------------------------------------------------------- C##### RETURN END