MODULE etat0_dcmip4_mod USE icosa PRIVATE REAL(rstd),PARAMETER :: eta0=0.252 REAL(rstd),PARAMETER :: etat=0.2 REAL(rstd),PARAMETER :: ps0=1e5 REAL(rstd),PARAMETER :: u0=35 REAL(rstd),PARAMETER :: T0=288 REAL(rstd),PARAMETER :: DeltaT=4.8e5 REAL(rstd),PARAMETER :: Rd=287 REAL(rstd),PARAMETER :: Gamma=0.005 REAL(rstd),PARAMETER :: up0=1 REAL(rstd) :: latw=2*Pi/9 REAL(rstd) :: pw=34000 REAL(rstd) :: q0=0.021 REAL(rstd) :: lonc REAL(rstd) :: latc PUBLIC etat0 CONTAINS SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) USE icosa IMPLICIT NONE TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_phis(:) TYPE(t_field),POINTER :: f_theta_rhodz(:) TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_q(:) REAL(rstd),POINTER :: ps(:) REAL(rstd),POINTER :: phis(:) REAL(rstd),POINTER :: theta_rhodz(:,:) REAL(rstd),POINTER :: u(:,:) REAL(rstd),POINTER :: q(:,:,:) INTEGER :: ind DO ind=1,ndomain CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) phis=f_phis(ind) theta_rhodz=f_theta_rhodz(ind) u=f_u(ind) q=f_q(ind) CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) ENDDO END SUBROUTINE etat0 SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) USE icosa USE disvert_mod USE pression_mod USE exner_mod USE geopotential_mod USE theta2theta_rhodz_mod IMPLICIT NONE REAL(rstd),INTENT(OUT) :: ps(iim*jjm) REAL(rstd),INTENT(OUT) :: phis(iim*jjm) REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) INTEGER :: i,j,l,ij REAL(rstd) :: theta(iim*jjm,llm) REAL(rstd) :: Y(iim*jjm,llm) REAL(rstd) :: vort REAL(rstd) :: eta(llm) REAL(rstd) :: etav(llm) REAL(rstd) :: etas, etavs REAL(rstd) :: lon,lat REAL(rstd) :: ulon(3) REAL(rstd) :: ep(3), norm_ep REAL(rstd) :: Tave, T REAL(rstd) :: phis_ave REAL(rstd) :: V0(3) REAL(rstd) :: r2 REAL(rstd) :: utot REAL(rstd) :: lonx,latx REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r INTEGER :: testcase lonc=Pi/9 latc=2*Pi/9 testcase=1 CALL getin("dcmip4_testcase",testcase) DO l=1,llm eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) ) etav(l)=(eta(l)-eta0)*Pi/2 ENDDO etas=ap(1)/preff+bp(1) etavs=(etas-eta0)*Pi/2 DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i ps(ij)=ps0 ENDDO ENDDO CALL lonlat2xyz(lonc,latc,V0) u(:,:)=1e10 DO l=1,llm DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) r=radius*acos(K) utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:)) CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) r=radius*acos(K) utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:)) CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) r=radius*acos(K) utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:)) ENDDO ENDDO ENDDO DO l=1,llm Tave=T0*eta(l)**(Rd*Gamma/g) IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) theta(ij,l)=T*eta(l)**(-kappa) ENDDO ENDDO ENDDO phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5 & +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) ! phis(ij)=phis_ave+u0*cos(etavs)**1.5 ENDDO ENDDO CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) IF (testcase==1) THEN q(:,:,1)=theta(:,:) DO l=1,llm dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1) IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa) & + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1)) DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l) & - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)& - 9./8. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l)) & * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5 & *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) & + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat))) duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l)) K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) r=radius*acos(K) vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2) & + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) & -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) ENDDO ENDDO ENDDO ELSE IF (testcase==2) THEN DO l=1,llm DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE compute_etat0_dcmip4 END MODULE etat0_dcmip4_mod