source: codes/icosagcm/devel/src/initial/etat0_jablonowsky06.f90 @ 834

Last change on this file since 834 was 834, checked in by dubos, 5 years ago

devel : moved thermodynamics from caldyn_vars to earth_const
File size: 2.3 KB
Line 
1MODULE etat0_jablonowsky06_mod
2  USE icosa, ignore_Rd=>Rd ! avoid name conflict with Rd from earth_const                                                                                                                                    PRIVATE
3  REAL(rstd),PARAMETER :: eta0=0.252
4  REAL(rstd),PARAMETER :: etat=0.2
5  REAL(rstd),PARAMETER :: ps0=1e5
6  REAL(rstd),PARAMETER :: u0=35
7  REAL(rstd),PARAMETER :: T0=288
8  REAL(rstd),PARAMETER :: DeltaT=4.8e5
9  REAL(rstd),PARAMETER :: Rd=287
10  REAL(rstd),PARAMETER :: Gamma=0.005
11  REAL(rstd),PARAMETER :: up0=1
12
13  PUBLIC compute_etat0
14
15CONTAINS
16 
17  SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, temp, ulon, ulat)
18    USE disvert_mod
19    IMPLICIT NONE 
20    INTEGER, INTENT(IN) :: ngrid
21    REAL(rstd),INTENT(IN) :: lat(ngrid)
22    REAL(rstd),INTENT(IN) :: lon(ngrid)
23    REAL(rstd),INTENT(OUT) :: phis(ngrid)
24    REAL(rstd),INTENT(OUT) :: ps(ngrid)
25    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
26    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
27    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
28   
29    INTEGER :: l,ij
30    REAL(rstd) :: eta(llm)
31    REAL(rstd) :: etav(llm)
32    REAL(rstd) :: etas, etavs, Tave, phis_ave, r2
33   
34    DO l=1,llm
35       eta(l)  = 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
36       etav(l) = (eta(l)-eta0)*Pi/2
37    ENDDO
38    etas  = ap(1)+bp(1)
39    etavs = (etas-eta0)*Pi/2
40       
41    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
42    DO ij=1,ngrid
43       ps(ij)=ps0
44       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  &
45            +(8./5*cos(lat(ij))**3 * (sin(lat(ij))**2 + 2./3) - Pi/4)*radius*Omega )
46    ENDDO
47
48    DO l=1,llm
49       Tave=T0*eta(l)**(Rd*Gamma/g)
50       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
51       DO ij=1,ngrid
52          r2 = arc(Pi/9.,2.*Pi/9., lon(ij),lat(ij))**2
53          temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5  &
54               * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
55               + (8./5*cos(lat(ij))**3*(sin(lat(ij))**2+2./3)-Pi/4)*radius*Omega)
56          ulon(ij,l) = u0*cos(etav(l))**1.5*sin(2*lat(ij))**2 + up0*exp(-r2/0.01)
57          ulat(ij,l) = 0
58       ENDDO
59    ENDDO
60   
61   
62  END SUBROUTINE compute_etat0
63 
64END MODULE etat0_jablonowsky06_mod
Note: See TracBrowser for help on using the repository browser.