MODULE etat0_venus_mod USE icosa IMPLICIT NONE PRIVATE TYPE(t_field),POINTER :: f_temp_eq( :) TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics REAL(rstd), SAVE :: kfrict !$OMP THREADPRIVATE(kfrict) PUBLIC :: etat0, init_physics, physics CONTAINS !-------------------------------- "Physics" ---------------------------------------- SUBROUTINE init_physics USE getin_mod IMPLICIT NONE REAL(rstd),POINTER :: temp(:,:) REAL(rstd) :: friction_time INTEGER :: ind friction_time=86400. !friction CALL getin('friction_time',friction_time) kfrict=1./friction_time CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use PRINT *, 'Initializing Temp_eq (venus)' CALL allocate_field(f_temp_eq,field_t,type_real,llm) DO ind=1,ndomain IF (.NOT. assigned_domain(ind)) CYCLE CALL swap_dimensions(ind) CALL swap_geometry(ind) temp=f_temp_eq(ind) CALL compute_temp_ref(temp, .TRUE.) ! FIXME With meridional gradient ENDDO END SUBROUTINE init_physics SUBROUTINE physics(f_ps,f_theta_rhodz,f_u) USE icosa USE theta2theta_rhodz_mod IMPLICIT NONE TYPE(t_field),POINTER :: f_theta_rhodz(:) TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_ps(:) REAL(rstd),POINTER :: temp(:,:) REAL(rstd),POINTER :: temp_eq(:,:) REAL(rstd),POINTER :: u(:,:) INTEGER :: ind CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp) DO ind=1,ndomain IF (.NOT. assigned_domain(ind)) CYCLE CALL swap_dimensions(ind) CALL swap_geometry(ind) u=f_u(ind) temp_eq=f_temp_eq(ind) temp=f_temp(ind) CALL compute_physics(temp_eq, temp, u) ENDDO CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) END SUBROUTINE physics SUBROUTINE compute_physics(temp_eq, temp, u) USE icosa USE theta2theta_rhodz_mod IMPLICIT NONE REAL(rstd),INTENT(IN) :: temp_eq(iim*jjm,llm) REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm) REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm) REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012 INTEGER :: i,j,l,ij 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 temp(ij,l) = temp(ij,l) - (temp(ij,l)-temp_eq(ij,l))*(dt*itau_physics/tauCLee) ENDDO ENDDO ENDDO u(:,1)=u(:,1)*(1.-dt*itau_physics*kfrict) END SUBROUTINE compute_physics !----------------------------- Initialize to T_eq -------------------------------------- SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) USE icosa USE theta2theta_rhodz_mod 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(:) TYPE(t_field),POINTER :: f_temp(:) REAL(rstd),POINTER :: temp(:,:) REAL(rstd),POINTER :: ps(:) REAL(rstd),POINTER :: phis(:) REAL(rstd),POINTER :: u(:,:) REAL(rstd),POINTER :: q(:,:,:) REAL(rstd) :: lat(iim*jjm) ! latitude REAL(rstd) :: pplay(iim*jjm, llm) ! pressure at full layers INTEGER :: ind CALL allocate_field(f_temp,field_t,type_real,llm) DO ind=1,ndomain IF (.NOT. assigned_domain(ind)) CYCLE CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) ps(:)=preff phis=f_phis(ind) phis(:)=0. u=f_u(ind) u(:,:)=0 q=f_q(ind) q(:,:,:)=1e2 temp=f_temp(ind) CALL compute_temp_ref(temp, .FALSE.) ! Without meridional gradient ENDDO CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) CALL deallocate_field(f_temp) END SUBROUTINE etat0 !------------------------- Compute reference temperature field ------------------------ SUBROUTINE compute_temp_ref(theta_eq, gradient) USE icosa USE disvert_mod USE omp_para USE math_const IMPLICIT NONE REAL(rstd), INTENT(OUT) :: theta_eq(iim*jjm,llm) LOGICAL, INTENT(IN) :: gradient REAL(rstd) :: clat(iim*jjm) ! latitude integer :: level REAL :: pressCLee(31), tempCLee(31), dt_epCLee(31), etaCLee(31) real(rstd) :: lon,lat, pplay, ztemp,zdt,fact logical, save :: firstcall integer :: i,j,ij, l,ll data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, & 4.273e-1, 3.373e-1, 2.610e-1,1.979e-1,1.472e-1, & 1.074e-1, 7.672e-2, 5.361e-2,3.657e-2,2.430e-2, & 1.569e-2, 9.814e-3, 5.929e-3,3.454e-3,1.934e-3, & 1.043e-3, 5.400e-4, 2.710e-4,1.324e-4,6.355e-5, & 3.070e-5, 1.525e-5, 7.950e-6,4.500e-6,2.925e-6, & 2.265e-6/ 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,:),lon,lat) clat(ij)=cos(lat) ENDDO ENDDO data tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, & 602.225, 574.542, 546.104, 517.339, 488.560, 459.932, & 431.741, 404.202, 377.555, 352.042, 327.887, 305.313, & 284.556, 265.697, 248.844, 233.771, 220.368, 208.247, & 197.127, 187.104, 178.489, 171.800, 167.598, 165.899, & 165.676/ data dt_epCLee/6.101 , 6.136 , 6.176 , 6.410 , 6.634 , 6.678 , & 6.719 , 6.762 , 7.167 , 7.524 , 9.840 ,14.948 , & 21.370 ,28.746 ,36.373 ,43.315 ,48.534 ,51.175 , & 50.757 ,47.342 ,41.536 ,34.295 ,26.758 ,19.807 , & 14.001 , 9.599 , 6.504 , 4.439 , 3.126 , 2.370 , & 2.000/ pressCLee = etaCLee*9.2e6 DO j=jj_begin-1,jj_end+1 DO i=ii_begin-1,ii_end+1 ij=(j-1)*iim+i DO l = 1, llm pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff ! look for largest level such that pressCLee(level) > pplay(ij,l)) ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) level = 1 DO ll = 1, 30 ! 30 data levels IF(pressCLee(ll) > pplay) THEN level = ll END IF END DO ! interpolate between level and level+1 ! interpolation is linear in log(pressure) fact = ( log10(pplay)-log10(pressCLee(level))) & /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact zdt = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact IF(gradient) THEN theta_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) ELSE theta_eq(ij,l) = ztemp END IF END DO END DO END DO END SUBROUTINE compute_temp_ref END MODULE etat0_venus_mod