source: codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90

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

devel : Cp(T) thermodynamics (TBC)
File size: 4.7 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- compute_geopot ----------------------------------
3   SELECT CASE(caldyn_thermo)
4   CASE(thermo_boussinesq)
5      ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure)
6      !$OMP DO SCHEDULE(STATIC)
7      DO ij=1,primal_num
8         pk(llm,ij) = ptop + .5*g* theta(llm,ij,1)*rhodz(llm,ij)
9         DO l = llm-1,1,-1
10            pk(l,ij) = pk(l+1,ij) + (.5*g)*( theta(l,ij,1)*rhodz(l,ij) + theta(l+1,ij,1)*rhodz(l+1,ij) )
11         END DO
12         IF(caldyn_eta == eta_lag) THEN
13            ps(ij) = pk(1,ij) + .5*g* theta(1,ij,1)*rhodz(1,ij)
14         END IF
15      END DO
16      !$OMP END DO
17      ! now pk contains the Lagrange multiplier (pressure)
18      ! specific volume 1 = dphi/g/rhodz
19      !$OMP DO SCHEDULE(STATIC)
20      DO ij=1,primal_num
21         DO l = 1,llm
22            geopot(l+1,ij) = geopot(l,ij) + g*rhodz(l,ij)
23         END DO
24      END DO
25      !$OMP END DO
26   CASE(thermo_theta)
27      !$OMP DO SCHEDULE(STATIC)
28      DO ij=1,primal_num
29         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
30         DO l = llm-1,1,-1
31            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
32         END DO
33         IF(caldyn_eta == eta_lag) THEN
34            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
35         END IF
36      END DO
37      !$OMP END DO
38      !$OMP DO SCHEDULE(STATIC)
39      DO ij=1,primal_num
40         DO l = 1,llm
41            p_ik = pk(l,ij)
42            exner_ik = cpp * (p_ik/preff) ** kappa
43            gv = (g*kappa)*theta(l,ij,1)*exner_ik/p_ik
44            pk(l,ij) = exner_ik
45            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
46         END DO
47      END DO
48      !$OMP END DO
49   CASE(thermo_entropy)
50      !$OMP DO SCHEDULE(STATIC)
51      DO ij=1,primal_num
52         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
53         DO l = llm-1,1,-1
54            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
55         END DO
56         IF(caldyn_eta == eta_lag) THEN
57            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
58         END IF
59      END DO
60      !$OMP END DO
61      !$OMP DO SCHEDULE(STATIC)
62      DO ij=1,primal_num
63         DO l = 1,llm
64            p_ik = pk(l,ij)
65            temp_ik = Treff*exp((theta(l,ij,1) + Rd*log(p_ik/preff))/cpp)
66            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
67            pk(l,ij) = temp_ik
68            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
69         END DO
70      END DO
71      !$OMP END DO
72   CASE(thermo_variable_Cp)
73      ! thermodynamics with variable Cp
74      ! Cp.dT = dh = Tds + vdp
75      ! pv = RT
76      ! => ds = (dh+v.dp)/T = Cp.dT/T - R dp/p
77      ! Cp(T) = Cp0 * (T/T0)^nu
78      ! => s(p,T) = Cp(T)/nu - R log(p/preff)
79      ! h = Cp(T).T/(nu+1)
80      !$OMP DO SCHEDULE(STATIC)
81      DO ij=1,primal_num
82         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
83         DO l = llm-1,1,-1
84            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
85         END DO
86         IF(caldyn_eta == eta_lag) THEN
87            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
88         END IF
89      END DO
90      !$OMP END DO
91      !$OMP DO SCHEDULE(STATIC)
92      DO ij=1,primal_num
93         DO l = 1,llm
94            p_ik = pk(l,ij)
95            Cp_ik = nu*( theta(l,ij,1) + Rd*log(p_ik/preff) )
96            temp_ik = Treff* (Cp_ik/cpp)**(1./nu)
97            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
98            pk(l,ij) = temp_ik
99            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
100         END DO
101      END DO
102      !$OMP END DO
103   CASE(thermo_moist)
104      !$OMP DO SCHEDULE(STATIC)
105      DO ij=1,primal_num
106         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)*(1.+theta(llm,ij,2))
107         DO l = llm-1,1,-1
108            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij)*(1.+theta(l,ij,2)) + rhodz(l+1,ij)*(1.+theta(l+1,ij,2)) )
109         END DO
110         IF(caldyn_eta == eta_lag) THEN
111            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)*(1.+theta(1,ij,2))
112         END IF
113      END DO
114      !$OMP END DO
115      !$OMP DO SCHEDULE(STATIC)
116      DO ij=1,primal_num
117         DO l = 1,llm
118            p_ik = pk(l,ij)
119            qv = theta(l,ij,2) ! water vaper mixing ratio = mv/md
120            Rmix = Rd+qv*Rv
121            chi = ( theta(l,ij,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
122            temp_ik = Treff*exp(chi)
123            ! specific volume v = R*T/p
124            ! R = (Rd + qv.Rv)/(1+qv)
125            gv = g*Rmix*temp_ik/(p_ik*(1+qv))
126            pk(l,ij) = temp_ik
127            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
128         END DO
129      END DO
130      !$OMP END DO
131   END SELECT
132   !---------------------------- compute_geopot ----------------------------------
133   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.