module coare_cp_mod implicit none private public psit_30, psiuo, coare_cp contains real function psit_30(zet) implicit none real, intent(in) :: zet real :: x, psik, psic, f, c if(zet<0) then x=(1.-(15.*zet))**.5 psik=2.*log((1.+x)/2) x=(1.-(34.15*zet))**.3333 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) f=zet*zet/(1.+zet*zet) psit_30=(1.-f)*psik+f*psic else c=min(50.,.35*zet) psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) endif end function psit_30 real function psiuo(zet) implicit none real, intent(in) :: zet real :: x, psik, psic, f, c if (zet<0) then x=(1.-15.*zet)**.25 psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) x=(1.-10.15*zet)**.3333 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) f=zet*zet/(1+zet*zet) psiuo=(1-f)*psik+f*psic else c=min(50.,.35*zet) psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) endif end function psiuo subroutine coare_cp(du,dt,dq,& t,q,& zu,zt,zq,& p,zi,& nits,& coeffs,rugosm,rugosh) ! zo,tau,hsb,hlb,& ! var) !version with shortened iteration modified Rt and Rq !uses wave information wave period in s and wave ht in m !no wave, standard coare 2.6 charnock: jwave=0 !Oost et al. zo=50/2/pi L (u*/c)**4.5 if jwave=1 !taylor and yelland zo=1200 h*(L/h)**4.5 jwave=2 ! USE MODD_CSTS, ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, & XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, & XP00 IMPLICIT NONE real, intent(in) :: du,dt,dq,t,q real, intent(in) :: zu,zt,zq,p,zi integer, intent(in) :: nits ! real, dimension (nits), intent(out) :: zo,tau,hsb,hlb ! real, dimension(nits,3), intent(out) :: var real, dimension(3), intent(out) :: coeffs real, intent(out) :: rugosm real, intent(out) :: rugosh real, parameter :: beta=1.2, von=.4, fdg = 1. ,& tdk = 273.16, pi = 3.141593, grav = 9.82,& rgas = 287.1 integer, dimension(3) :: shape_input real bf, cc, rhoa, visa,& u10, ut, uts, ut0, ug,& cd10, ch10, ct, ct10,& charn, ribu,& l, l10, & ribcu,rr,& zet,zetu, zom10, zoh10, zot, zoq, zot10,& cd, ch, le, cpa, cpv,& usr,qsr,tsr,& zom, t_c !, ZTWAVE, ZHWAVE, ZCWAVE, ZLWAVE real old_usr, old_tsr, old_qsr,tmp real, external :: grv integer i,j,k !---------------- Rajout pour prendre en compte différent Z0 --------------------------------! ! INTEGER :: NGRVWAVES ! Pour le choix du z0 ! NGRVWAVES = 2 !--------------------------------------------------------------------------------------------- !-------------------- Attention Modif réalisée pas SURRR ------------------------------------- !--------------------------------------------------------------------------------------------- Ribcu=-zu/zi/.004/Beta**3 cpa=1004.67 t_c = t - 273.15 Le=(2.501-.00237*(t_c-dt))*1e6 visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) cpv=cpa*(1+0.84*Q) rhoa=P/(Rgas*t*(1+0.61*Q)) ug=.2 ut=sqrt(du*du+ug*ug) ut=MAX(ut , 0.1 * MIN(10.,zu) ) u10=ut*log(10/1e-4)/log(zu/1e-4) usr=.035*u10 ! turbulent friction velocity (m/s), including gustiness zom10=0.011*usr*usr/grav+0.11*visa/usr ! roughness length for u (smith 88) Cd10=(von/log(10/zom10))**2 ! zoh10=0.40*visa/usr+1.4e-5 ! Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca Ch10=0.00115 Ct10=Ch10/sqrt(Cd10) zot10=10/exp(von/Ct10) ! roughness length for t Cd=(von/log(zu/zom10))**2 Ct=von/log(zt/zot10) CC=von*Ct/Cd ut0 = ut ut = ut0 Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 if (Ribu < 0) then zetu=CC*Ribu/(1+Ribu/Ribcu) else zetu=CC*Ribu*(1+27/9*Ribu/CC) endif L10=zu/zetu ! if (zetu .GT. 50) then ! nits=1 ! endif usr=ut*von/(log(zu/zom10)-psiuo(zetu)) tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) ! tkt=.001 ! charnock constant - lin par morceau - constant if ( ut <= 10. ) then charn=0.011 else if (ut > 18) then charn=0.018 else charn=0.011+(ut-10)/(18-10)*(0.018-0.011) endif endif ! ZHWAVE = 0.018*ut*ut*(1.+0.015*ut) ! ZTWAVE = 0.729*ut ! ZCWAVE = XG*ZTWAVE/(2.*pi) ! ZLWAVE = ZTWAVE*ZCWAVE !*************** bulk loop ************ do i=1, nits zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q) ! IF (NGRVWAVES==0) THEN ! zom = charn*usr*usr/XG + 0.11*visa/usr !Smith 1988 ! ELSE IF (NGRVWAVES==1) THEN ! zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 & ! + 0.11*visa/usr !Oost et al. 2002 ! ELSE IF (NGRVWAVES==2) THEN ! zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 & ! + 0.11*visa/usr !Taulor and Yelland 2001 ! ENDIF zom=charn*usr*usr/grav+0.11*visa/usr rr=zom*usr/visa L=zu/zet zoq=min(1.15e-4,5.5e-5/rr**.6) ! a modifier zot=zoq ! a modifier ! zot=0.40*visa/usr+1.4e-5 old_usr = usr old_tsr = tsr old_qsr = qsr usr=ut*von/(log(zu/zom)-psiuo(zu/L)) tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) Bf=-grav/t*usr*(tsr+.61*t*qsr) if (Bf > 0) then ug=Beta*(Bf*zi)**.333 else ug=.2 endif ut=sqrt(du*du+ug*ug) ut=MAX(ut , 0.1 * MIN(10.,zu) ) enddo !bulk iter loop ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du coeffs = [tmp**2,& von*fdg/(log(zt/zot)-psit_30(zt/L)) * tmp,& von*fdg/(log(zq/zoq)-psit_30(zq/L)) * tmp] rugosm = zom rugosh = zot end subroutine coare_cp end module coare_cp_mod