| 1 | module coare_cp_mod |
|---|
| 2 | implicit none |
|---|
| 3 | private |
|---|
| 4 | public psit_30, psiuo, coare_cp |
|---|
| 5 | |
|---|
| 6 | contains |
|---|
| 7 | |
|---|
| 8 | real function psit_30(zet) |
|---|
| 9 | implicit none |
|---|
| 10 | real, intent(in) :: zet |
|---|
| 11 | |
|---|
| 12 | real :: x, psik, psic, f, c |
|---|
| 13 | |
|---|
| 14 | if(zet<0) then |
|---|
| 15 | x=(1.-(15.*zet))**.5 |
|---|
| 16 | psik=2.*log((1.+x)/2) |
|---|
| 17 | x=(1.-(34.15*zet))**.3333 |
|---|
| 18 | psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) |
|---|
| 19 | f=zet*zet/(1.+zet*zet) |
|---|
| 20 | psit_30=(1.-f)*psik+f*psic |
|---|
| 21 | else |
|---|
| 22 | c=min(50.,.35*zet) |
|---|
| 23 | psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) |
|---|
| 24 | endif |
|---|
| 25 | |
|---|
| 26 | end function psit_30 |
|---|
| 27 | |
|---|
| 28 | real function psiuo(zet) |
|---|
| 29 | implicit none |
|---|
| 30 | real, intent(in) :: zet |
|---|
| 31 | |
|---|
| 32 | real :: x, psik, psic, f, c |
|---|
| 33 | |
|---|
| 34 | if (zet<0) then |
|---|
| 35 | x=(1.-15.*zet)**.25 |
|---|
| 36 | psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) |
|---|
| 37 | x=(1.-10.15*zet)**.3333 |
|---|
| 38 | psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) |
|---|
| 39 | f=zet*zet/(1+zet*zet) |
|---|
| 40 | psiuo=(1-f)*psik+f*psic |
|---|
| 41 | else |
|---|
| 42 | c=min(50.,.35*zet) |
|---|
| 43 | psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) |
|---|
| 44 | endif |
|---|
| 45 | |
|---|
| 46 | end function psiuo |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | subroutine coare_cp(du,dt,dq,& |
|---|
| 50 | t,q,& |
|---|
| 51 | zu,zt,zq,& |
|---|
| 52 | p,zi,& |
|---|
| 53 | nits,& |
|---|
| 54 | coeffs,rugosm,rugosh) |
|---|
| 55 | ! zo,tau,hsb,hlb,& |
|---|
| 56 | ! var) |
|---|
| 57 | !version with shortened iteration modified Rt and Rq |
|---|
| 58 | !uses wave information wave period in s and wave ht in m |
|---|
| 59 | !no wave, standard coare 2.6 charnock: jwave=0 |
|---|
| 60 | !Oost et al. zo=50/2/pi L (u*/c)**4.5 if jwave=1 |
|---|
| 61 | !taylor and yelland zo=1200 h*(L/h)**4.5 jwave=2 |
|---|
| 62 | ! |
|---|
| 63 | USE MODD_CSTS, ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, & |
|---|
| 64 | XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, & |
|---|
| 65 | XP00 |
|---|
| 66 | |
|---|
| 67 | IMPLICIT NONE |
|---|
| 68 | |
|---|
| 69 | real, intent(in) :: du,dt,dq,t,q |
|---|
| 70 | real, intent(in) :: zu,zt,zq,p,zi |
|---|
| 71 | integer, intent(in) :: nits |
|---|
| 72 | ! real, dimension (nits), intent(out) :: zo,tau,hsb,hlb |
|---|
| 73 | ! real, dimension(nits,3), intent(out) :: var |
|---|
| 74 | real, dimension(3), intent(out) :: coeffs |
|---|
| 75 | real, intent(out) :: rugosm |
|---|
| 76 | real, intent(out) :: rugosh |
|---|
| 77 | |
|---|
| 78 | real, parameter :: beta=1.2, von=.4, fdg = 1. ,& |
|---|
| 79 | tdk = 273.16, pi = 3.141593, grav = 9.82,& |
|---|
| 80 | rgas = 287.1 |
|---|
| 81 | |
|---|
| 82 | integer, dimension(3) :: shape_input |
|---|
| 83 | |
|---|
| 84 | real bf, cc, rhoa, visa,& |
|---|
| 85 | u10, ut, uts, ut0, ug,& |
|---|
| 86 | cd10, ch10, ct, ct10,& |
|---|
| 87 | charn, ribu,& |
|---|
| 88 | l, l10, & |
|---|
| 89 | ribcu,rr,& |
|---|
| 90 | zet,zetu, zom10, zoh10, zot, zoq, zot10,& |
|---|
| 91 | cd, ch, le, cpa, cpv,& |
|---|
| 92 | usr,qsr,tsr,& |
|---|
| 93 | zom, t_c |
|---|
| 94 | |
|---|
| 95 | !, ZTWAVE, ZHWAVE, ZCWAVE, ZLWAVE |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | |
|---|
| 100 | real old_usr, old_tsr, old_qsr,tmp |
|---|
| 101 | |
|---|
| 102 | real, external :: grv |
|---|
| 103 | integer i,j,k |
|---|
| 104 | !---------------- Rajout pour prendre en compte différent Z0 --------------------------------! |
|---|
| 105 | ! INTEGER :: NGRVWAVES ! Pour le choix du z0 |
|---|
| 106 | ! NGRVWAVES = 2 |
|---|
| 107 | !--------------------------------------------------------------------------------------------- |
|---|
| 108 | !-------------------- Attention Modif réalisée pas SURRR ------------------------------------- |
|---|
| 109 | !--------------------------------------------------------------------------------------------- |
|---|
| 110 | |
|---|
| 111 | Ribcu=-zu/zi/.004/Beta**3 |
|---|
| 112 | cpa=1004.67 |
|---|
| 113 | |
|---|
| 114 | t_c = t - 273.15 |
|---|
| 115 | |
|---|
| 116 | |
|---|
| 117 | Le=(2.501-.00237*(t_c-dt))*1e6 |
|---|
| 118 | |
|---|
| 119 | |
|---|
| 120 | visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | cpv=cpa*(1+0.84*Q) |
|---|
| 124 | rhoa=P/(Rgas*t*(1+0.61*Q)) |
|---|
| 125 | |
|---|
| 126 | |
|---|
| 127 | |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | |
|---|
| 131 | ug=.2 |
|---|
| 132 | |
|---|
| 133 | ut=sqrt(du*du+ug*ug) |
|---|
| 134 | ut=MAX(ut , 0.1 * MIN(10.,zu) ) |
|---|
| 135 | |
|---|
| 136 | u10=ut*log(10/1e-4)/log(zu/1e-4) |
|---|
| 137 | usr=.035*u10 ! turbulent friction velocity (m/s), including gustiness |
|---|
| 138 | zom10=0.011*usr*usr/grav+0.11*visa/usr ! roughness length for u (smith 88) |
|---|
| 139 | Cd10=(von/log(10/zom10))**2 |
|---|
| 140 | ! zoh10=0.40*visa/usr+1.4e-5 |
|---|
| 141 | ! Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca |
|---|
| 142 | Ch10=0.00115 |
|---|
| 143 | Ct10=Ch10/sqrt(Cd10) |
|---|
| 144 | zot10=10/exp(von/Ct10) ! roughness length for t |
|---|
| 145 | Cd=(von/log(zu/zom10))**2 |
|---|
| 146 | Ct=von/log(zt/zot10) |
|---|
| 147 | CC=von*Ct/Cd |
|---|
| 148 | |
|---|
| 149 | ut0 = ut |
|---|
| 150 | |
|---|
| 151 | |
|---|
| 152 | ut = ut0 |
|---|
| 153 | Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 |
|---|
| 154 | |
|---|
| 155 | if (Ribu < 0) then |
|---|
| 156 | zetu=CC*Ribu/(1+Ribu/Ribcu) |
|---|
| 157 | else |
|---|
| 158 | zetu=CC*Ribu*(1+27/9*Ribu/CC) |
|---|
| 159 | endif |
|---|
| 160 | |
|---|
| 161 | L10=zu/zetu |
|---|
| 162 | |
|---|
| 163 | ! if (zetu .GT. 50) then |
|---|
| 164 | ! nits=1 |
|---|
| 165 | ! endif |
|---|
| 166 | |
|---|
| 167 | usr=ut*von/(log(zu/zom10)-psiuo(zetu)) |
|---|
| 168 | tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) |
|---|
| 169 | qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) |
|---|
| 170 | ! tkt=.001 |
|---|
| 171 | |
|---|
| 172 | ! charnock constant - lin par morceau - constant |
|---|
| 173 | if ( ut <= 10. ) then |
|---|
| 174 | charn=0.011 |
|---|
| 175 | else |
|---|
| 176 | if (ut > 18) then |
|---|
| 177 | charn=0.018 |
|---|
| 178 | else |
|---|
| 179 | charn=0.011+(ut-10)/(18-10)*(0.018-0.011) |
|---|
| 180 | endif |
|---|
| 181 | endif |
|---|
| 182 | |
|---|
| 183 | ! ZHWAVE = 0.018*ut*ut*(1.+0.015*ut) |
|---|
| 184 | ! ZTWAVE = 0.729*ut |
|---|
| 185 | ! ZCWAVE = XG*ZTWAVE/(2.*pi) |
|---|
| 186 | ! ZLWAVE = ZTWAVE*ZCWAVE |
|---|
| 187 | |
|---|
| 188 | |
|---|
| 189 | !*************** bulk loop ************ |
|---|
| 190 | do i=1, nits |
|---|
| 191 | |
|---|
| 192 | zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q) |
|---|
| 193 | |
|---|
| 194 | ! IF (NGRVWAVES==0) THEN |
|---|
| 195 | ! zom = charn*usr*usr/XG + 0.11*visa/usr !Smith 1988 |
|---|
| 196 | ! ELSE IF (NGRVWAVES==1) THEN |
|---|
| 197 | ! zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 & |
|---|
| 198 | ! + 0.11*visa/usr !Oost et al. 2002 |
|---|
| 199 | ! ELSE IF (NGRVWAVES==2) THEN |
|---|
| 200 | ! zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 & |
|---|
| 201 | ! + 0.11*visa/usr !Taulor and Yelland 2001 |
|---|
| 202 | ! ENDIF |
|---|
| 203 | |
|---|
| 204 | zom=charn*usr*usr/grav+0.11*visa/usr |
|---|
| 205 | |
|---|
| 206 | rr=zom*usr/visa |
|---|
| 207 | L=zu/zet |
|---|
| 208 | zoq=min(1.15e-4,5.5e-5/rr**.6) ! a modifier |
|---|
| 209 | zot=zoq ! a modifier |
|---|
| 210 | ! zot=0.40*visa/usr+1.4e-5 |
|---|
| 211 | |
|---|
| 212 | old_usr = usr |
|---|
| 213 | old_tsr = tsr |
|---|
| 214 | old_qsr = qsr |
|---|
| 215 | |
|---|
| 216 | usr=ut*von/(log(zu/zom)-psiuo(zu/L)) |
|---|
| 217 | tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) |
|---|
| 218 | qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) |
|---|
| 219 | |
|---|
| 220 | Bf=-grav/t*usr*(tsr+.61*t*qsr) |
|---|
| 221 | |
|---|
| 222 | if (Bf > 0) then |
|---|
| 223 | ug=Beta*(Bf*zi)**.333 |
|---|
| 224 | else |
|---|
| 225 | ug=.2 |
|---|
| 226 | endif |
|---|
| 227 | |
|---|
| 228 | ut=sqrt(du*du+ug*ug) |
|---|
| 229 | ut=MAX(ut , 0.1 * MIN(10.,zu) ) |
|---|
| 230 | |
|---|
| 231 | |
|---|
| 232 | enddo !bulk iter loop |
|---|
| 233 | |
|---|
| 234 | ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress |
|---|
| 235 | ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr |
|---|
| 236 | ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr |
|---|
| 237 | |
|---|
| 238 | tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du |
|---|
| 239 | |
|---|
| 240 | coeffs = [tmp**2,& |
|---|
| 241 | von*fdg/(log(zt/zot)-psit_30(zt/L)) * tmp,& |
|---|
| 242 | von*fdg/(log(zq/zoq)-psit_30(zq/L)) * tmp] |
|---|
| 243 | |
|---|
| 244 | rugosm = zom |
|---|
| 245 | rugosh = zot |
|---|
| 246 | |
|---|
| 247 | end subroutine coare_cp |
|---|
| 248 | |
|---|
| 249 | end module coare_cp_mod |
|---|
