source: codes/icosagcm/devel/src/kernels_hex/caldyn_fast.k90 @ 837

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

devel : Cp(T) thermodynamics (TBC)
File size: 2.5 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- caldyn_fast ----------------------------------
3   !
4   SELECT CASE(caldyn_thermo)
5   CASE(thermo_boussinesq)
6      DO l = ll_begin, ll_end
7         !DIR$ SIMD
8         DO ij=ij_begin, ij_end
9            berni(ij,l) = pk(ij,l)
10            ! from now on pk contains the vertically-averaged geopotential
11            pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
12         END DO
13      END DO
14   CASE(thermo_theta)
15      DO l = ll_begin, ll_end
16         !DIR$ SIMD
17         DO ij=ij_begin, ij_end
18            berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
19         END DO
20      END DO
21   CASE(thermo_entropy)
22      DO l = ll_begin, ll_end
23         !DIR$ SIMD
24         DO ij=ij_begin, ij_end
25            berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
26            berni(ij,l) = berni(ij,l) + pk(ij,l)*(cpp-theta(ij,l,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
27         END DO
28      END DO
29   CASE(thermo_variable_Cp)
30      ! thermodynamics with variable Cp
31      ! Cp(T) = Cp0 * (T/T0)^nu
32      ! => h = Cp(T).T/(nu+1)
33      DO l = ll_begin, ll_end
34         !DIR$ SIMD
35         DO ij=ij_begin, ij_end
36            berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
37            cp_ik = cpp*(pk(ij,l)/Treff)**nu
38            berni(ij,l) = berni(ij,l) + pk(ij,l)*(cp_ik/(nu+1.)-theta(ij,l,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
39         END DO
40      END DO
41   CASE DEFAULT
42      PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo ! FIXME
43      STOP
44   END SELECT
45   !
46   DO l = ll_begin, ll_end
47      !DIR$ SIMD
48      DO ij=ij_begin, ij_end
49         due = .5*(theta(ij,l,1)+theta(ij+t_right,l,1))*(pk(ij+t_right,l)-pk(ij,l)) + berni(ij+t_right,l)-berni(ij,l)
50         du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due
51         u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
52         due = .5*(theta(ij,l,1)+theta(ij+t_lup,l,1))*(pk(ij+t_lup,l)-pk(ij,l)) + berni(ij+t_lup,l)-berni(ij,l)
53         du(ij+u_lup,l) = du(ij+u_lup,l) - ne_lup*due
54         u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l)
55         due = .5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))*(pk(ij+t_ldown,l)-pk(ij,l)) + berni(ij+t_ldown,l)-berni(ij,l)
56         du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due
57         u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
58      END DO
59   END DO
60   !
61   !---------------------------- caldyn_fast ----------------------------------
62   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.