source: codes/icosagcm/trunk/src/kernels/compute_geopot.k90 @ 580

Last change on this file since 580 was 580, checked in by dubos, 7 years ago

trunk : upgrading to devel
File size: 3.9 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      DO ij=ij_omp_begin_ext,ij_omp_end_ext
7         pk(ij,llm) = ptop + .5*g* theta(ij,llm,1)*rhodz(ij,llm)
8      END DO
9      DO l = llm-1,1,-1
10         DO ij=ij_omp_begin_ext,ij_omp_end_ext
11            pk(ij,l) = pk(ij,l+1) + (.5*g)*( theta(ij,l,1)*rhodz(ij,l) + theta(ij,l+1,1)*rhodz(ij,l+1) )
12         END DO
13      END DO
14      IF(caldyn_eta == eta_lag) THEN
15         DO ij=ij_omp_begin_ext,ij_omp_end_ext
16            ps(ij) = pk(ij,1) + .5*g* theta(ij,1,1)*rhodz(ij,1)
17         END DO
18      END IF
19      ! now pk contains the Lagrange multiplier (pressure)
20      ! specific volume 1 = dphi/g/rhodz
21      DO l = 1,llm
22         DO ij=ij_omp_begin_ext,ij_omp_end_ext
23            geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
24         END DO
25      END DO
26   CASE(thermo_theta)
27      DO ij=ij_omp_begin_ext,ij_omp_end_ext
28         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)
29      END DO
30      DO l = llm-1,1,-1
31         DO ij=ij_omp_begin_ext,ij_omp_end_ext
32            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l) + rhodz(ij,l+1) )
33         END DO
34      END DO
35      IF(caldyn_eta == eta_lag) THEN
36         DO ij=ij_omp_begin_ext,ij_omp_end_ext
37            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)
38         END DO
39      END IF
40      DO l = 1,llm
41         DO ij=ij_omp_begin_ext,ij_omp_end_ext
42            p_ik = pk(ij,l)
43            exner_ik = cpp * (p_ik/preff) ** kappa
44            gv = (g*kappa)*theta(ij,l,1)*exner_ik/p_ik
45            pk(ij,l) = exner_ik
46            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
47         END DO
48      END DO
49   CASE(thermo_entropy)
50      DO ij=ij_omp_begin_ext,ij_omp_end_ext
51         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)
52      END DO
53      DO l = llm-1,1,-1
54         DO ij=ij_omp_begin_ext,ij_omp_end_ext
55            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l) + rhodz(ij,l+1) )
56         END DO
57      END DO
58      IF(caldyn_eta == eta_lag) THEN
59         DO ij=ij_omp_begin_ext,ij_omp_end_ext
60            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)
61         END DO
62      END IF
63      DO l = 1,llm
64         DO ij=ij_omp_begin_ext,ij_omp_end_ext
65            p_ik = pk(ij,l)
66            temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
67            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
68            pk(ij,l) = temp_ik
69            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
70         END DO
71      END DO
72   CASE(thermo_moist)
73      DO ij=ij_omp_begin_ext,ij_omp_end_ext
74         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)*(1.+theta(ij,llm,2))
75      END DO
76      DO l = llm-1,1,-1
77         DO ij=ij_omp_begin_ext,ij_omp_end_ext
78            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l)*(1.+theta(ij,l,2)) + rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) )
79         END DO
80      END DO
81      IF(caldyn_eta == eta_lag) THEN
82         DO ij=ij_omp_begin_ext,ij_omp_end_ext
83            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)*(1.+theta(ij,1,2))
84         END DO
85      END IF
86      DO l = 1,llm
87         DO ij=ij_omp_begin_ext,ij_omp_end_ext
88            p_ik = pk(ij,l)
89            qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
90            Rmix = Rd+qv*Rv
91            chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
92            temp_ik = Treff*exp(chi)
93            ! specific volume v = R*T/p
94            ! R = (Rd + qv.Rv)/(1+qv)
95            gv = g*Rmix*temp_ik/(p_ik*(1+qv))
96            pk(ij,l) = temp_ik
97            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
98         END DO
99      END DO
100   END SELECT
101   !---------------------------- compute_geopot ----------------------------------
102   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.