Changeset 78 for codes/icosagcm/trunk
- Timestamp:
- 08/02/12 23:14:56 (12 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0.f90
r67 r78 10 10 USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0 11 11 USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0 12 USE etat0_dcmip4 1_mod, ONLY : etat0_dcmip41=>etat012 USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 13 13 IMPLICIT NONE 14 14 TYPE(t_field),POINTER :: f_ps(:) … … 33 33 CASE ('dcmip3') 34 34 CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 35 CASE ('dcmip4 1')36 CALL etat0_dcmip4 1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)35 CASE ('dcmip4') 36 CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 37 37 CASE DEFAULT 38 38 PRINT*, 'Bad selector for variable etat0 <',etat0_type, & -
codes/icosagcm/trunk/src/etat0_dcmip41.f90
r67 r78 1 MODULE etat0_dcmip4 1_mod1 MODULE etat0_dcmip4_mod 2 2 USE icosa 3 3 PRIVATE … … 11 11 REAL(rstd),PARAMETER :: Gamma=0.005 12 12 REAL(rstd),PARAMETER :: up0=1 13 REAL(rstd) :: latw=2*Pi/9 14 REAL(rstd) :: pw=34000 15 REAL(rstd) :: q0=0.021 13 16 REAL(rstd) :: lonc 14 17 REAL(rstd) :: latc … … 42 45 u=f_u(ind) 43 46 q=f_q(ind) 44 CALL compute_etat0_dcmip4 1(ps, phis, theta_rhodz, u, q)47 CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) 45 48 ENDDO 46 49 47 50 END SUBROUTINE etat0 48 51 49 SUBROUTINE compute_etat0_dcmip4 1(ps, phis, theta_rhodz, u, q)52 SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) 50 53 USE icosa 51 54 USE disvert_mod … … 78 81 REAL(rstd) :: lonx,latx 79 82 REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r 80 lonc=Pi/9 81 latc=2*Pi/9 83 INTEGER :: testcase 84 85 lonc=Pi/9 86 latc=2*Pi/9 87 88 testcase=1 89 CALL getin("dcmip4_testcase",testcase) 90 82 91 83 92 DO l=1,llm … … 131 140 Tave=T0*eta(l)**(Rd*Gamma/g) 132 141 IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 133 DO j=jj_begin ,jj_end134 DO i=ii_begin ,ii_end142 DO j=jj_begin-1,jj_end+1 143 DO i=ii_begin-1,ii_end+1 135 144 ij=(j-1)*iim+i 136 145 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) … … 161 170 CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) 162 171 163 q(:,:,1)=theta(:,:) 164 165 166 DO l=1,llm 167 dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1) 168 IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa) & 172 173 IF (testcase==1) THEN 174 175 q(:,:,1)=theta(:,:) 176 177 DO l=1,llm 178 dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1) 179 IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa) & 169 180 + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1)) 170 DO j=jj_begin,jj_end171 DO i=ii_begin,ii_end172 ij=(j-1)*iim+i173 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)174 dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l)&175 + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l) &176 - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)&177 - 9./8. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l)) &181 DO j=jj_begin,jj_end 182 DO i=ii_begin,ii_end 183 ij=(j-1)*iim+i 184 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 185 dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & 186 + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l) & 187 - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)& 188 - 9./8. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l)) & 178 189 * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) 179 dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5 &180 *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) &181 + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat)))190 dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5 & 191 *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) & 192 + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat))) 182 193 183 duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l))194 duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l)) 184 195 185 K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)186 r=radius*acos(K)187 vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2) &188 + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) &196 K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 197 r=radius*acos(K) 198 vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2) & 199 + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) & 189 200 -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) 190 q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) 201 q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) 202 ENDDO 191 203 ENDDO 192 ENDDO 193 ENDDO 204 ENDDO 205 206 ELSE IF (testcase==2) THEN 207 208 DO l=1,llm 209 DO j=jj_begin,jj_end 210 DO i=ii_begin,ii_end 211 ij=(j-1)*iim+i 212 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 213 q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2) 214 ENDDO 215 ENDDO 216 ENDDO 217 218 ENDIF 194 219 195 220 196 END SUBROUTINE compute_etat0_dcmip4 1197 198 END MODULE etat0_dcmip4 1_mod221 END SUBROUTINE compute_etat0_dcmip4 222 223 END MODULE etat0_dcmip4_mod
Note: See TracChangeset
for help on using the changeset viewer.