MODULE etat0_jablonowsky06_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 PUBLIC test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06, compute_etat0_new CONTAINS SUBROUTINE test_etat0_jablonowsky06 USE icosa USE kinetic_mod USE pression_mod USE exner_mod USE geopotential_mod USE vorticity_mod IMPLICIT NONE TYPE(t_field),POINTER,SAVE :: f_ps(:) TYPE(t_field),POINTER,SAVE :: f_phis(:) TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:) TYPE(t_field),POINTER,SAVE :: f_u(:) TYPE(t_field),POINTER,SAVE :: f_Ki(:) TYPE(t_field),POINTER,SAVE :: f_temp(:) TYPE(t_field),POINTER,SAVE :: f_p(:) TYPE(t_field),POINTER,SAVE :: f_pks(:) TYPE(t_field),POINTER,SAVE :: f_pk(:) TYPE(t_field),POINTER,SAVE :: f_phi(:) TYPE(t_field),POINTER,SAVE :: f_vort(:) TYPE(t_field),POINTER,SAVE :: f_q(:) REAL(rstd),POINTER :: Ki(:,:) REAL(rstd),POINTER :: temp(:) INTEGER :: ind CALL allocate_field(f_ps,field_t,type_real) CALL allocate_field(f_phis,field_t,type_real) CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) CALL allocate_field(f_u,field_u,type_real,llm) CALL allocate_field(f_p,field_t,type_real,llm+1) CALL allocate_field(f_Ki,field_t,type_real,llm) CALL allocate_field(f_pks,field_t,type_real) CALL allocate_field(f_pk,field_t,type_real,llm) CALL allocate_field(f_phi,field_t,type_real,llm) CALL allocate_field(f_temp,field_t,type_real) CALL allocate_field(f_vort,field_z,type_real,llm) CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) CALL kinetic(f_u,f_Ki) CALL vorticity(f_u,f_vort) CALL pression(f_ps,f_p) CALL exner(f_ps,f_p,f_pks,f_pk) CALL geopotential(f_phis,f_pks,f_pk,f_theta_rhodz,f_phi) CALL writefield('ps',f_ps) CALL writefield('phis',f_phis) CALL writefield('theta',f_theta_rhodz) CALL writefield('f_phi',f_phi) CALL writefield('Ki',f_Ki) CALL writefield('vort',f_vort) END SUBROUTINE test_etat0_jablonowsky06 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 IF (.NOT. assigned_domain(ind)) CYCLE 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) q=0 CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) ENDDO END SUBROUTINE etat0 SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 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) INTEGER :: i,j,l,ij REAL(rstd) :: theta(iim*jjm,llm) 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 DO l=1,llm eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) ) etav(l)=(eta(l)-eta0)*Pi/2 ENDDO etas=ap(1)+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(Pi/9,2*Pi/9,V0) u(:,:)=1e10 DO l=1,llm DO j=jj_begin,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep) r2=(asin(sqrt(sum(ep*ep))))**2 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) ulon(1) = -sin(lon) * utot ulon(2) = cos(lon) * utot ulon(3) = 0 * utot CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep) norm_ep=sqrt(sum(ep(:)**2)) IF (norm_ep>1e-30) THEN ep=-ep/norm_ep*ne(ij,right) u(ij+u_right,l)=sum(ep(:)*ulon(:)) ENDIF CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep) r2=(asin(sqrt(sum(ep*ep))))**2 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) ulon(1) = -sin(lon) * utot ulon(2) = cos(lon) * utot ulon(3) = 0 * utot CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep) norm_ep=sqrt(sum(ep(:)**2)) IF (norm_ep>1e-30) THEN ep=-ep/norm_ep*ne(ij,lup) u(ij+u_lup,l)=sum(ep(:)*ulon(:)) ENDIF CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep) r2=(asin(sqrt(sum(ep*ep))))**2 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) ulon(1) = -sin(lon) * utot ulon(2) = cos(lon) * utot ulon(3) = 0 * utot CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep) norm_ep=sqrt(sum(ep(:)**2)) IF (norm_ep>1e-30) THEN ep=-ep/norm_ep*ne(ij,ldown) u(ij+u_ldown,l)=sum(ep(:)*ulon(:)) ENDIF 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,jj_end DO i=ii_begin,ii_end ij=(j-1)*iim+i CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 & * ( (-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) 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 ) ENDDO ENDDO CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) END SUBROUTINE compute_etat0_jablonowsky06 SUBROUTINE compute_etat0_new(ngrid,lat,lon, phis, ps, temp, ulon, ulat) USE disvert_mod IMPLICIT NONE INTEGER, INTENT(IN) :: ngrid REAL(rstd),INTENT(IN) :: lat(ngrid) REAL(rstd),INTENT(IN) :: lon(ngrid) REAL(rstd),INTENT(OUT) :: phis(ngrid) REAL(rstd),INTENT(OUT) :: ps(ngrid) REAL(rstd),INTENT(OUT) :: temp(ngrid,llm) REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) INTEGER :: l,ij REAL(rstd) :: eta(llm) REAL(rstd) :: etav(llm) REAL(rstd) :: etas, etavs, Tave, phis_ave, r2 DO l=1,llm eta(l) = 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) ) etav(l) = (eta(l)-eta0)*Pi/2 ENDDO etas = ap(1)+bp(1) etavs = (etas-eta0)*Pi/2 phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) DO ij=1,ngrid ps(ij)=ps0 phis(ij) = phis_ave + u0*cos(etavs)**1.5*( (-2*sin(lat(ij))**6 * (cos(lat(ij))**2+1./3) + 10./63 )*u0*cos(etavs)**1.5 & +(8./5*cos(lat(ij))**3 * (sin(lat(ij))**2 + 2./3) - Pi/4)*radius*Omega ) ENDDO DO l=1,llm Tave=T0*eta(l)**(Rd*Gamma/g) IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 DO ij=1,ngrid r2 = arc(Pi/9,2*Pi/9, lon(ij),lat(ij))**2 temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 & * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & + (8./5*cos(lat(ij))**3*(sin(lat(ij))**2+2./3)-Pi/4)*radius*Omega) ulon(ij,l) = u0*cos(etav(l))**1.5*sin(2*lat(ij))**2 + up0*exp(-r2/0.01) ulat(ij,l) = 0 ENDDO ENDDO END SUBROUTINE compute_etat0_new END MODULE etat0_jablonowsky06_mod