source: codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin @ 837

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

devel : Cp(T) thermodynamics (TBC)
File size: 6.0 KB
Line 
1KERNEL(caldyn_wflux)
2  SEQUENCE_C0
3    BODY('llm-1,1,-1')
4      ! cumulate mass flux convergence from top to bottom
5      convm(CELL) = convm(CELL) + convm(UP(CELL))
6    END_BLOCK
7    EPILOGUE(1)
8      dmass_col(HIDX(CELL)) = convm(CELL)
9    END_BLOCK
10    BODY('2,llm')
11      ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
12      wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL)
13    END_BLOCK
14  END_BLOCK
15  ! make sure wflux is up to date
16  BARRIER
17END_BLOCK
18
19KERNEL(caldyn_dmass)
20  FORALL_CELLS()
21    ON_PRIMAL
22      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL))
23    END_BLOCK
24  END_BLOCK
25END_BLOCK
26
27KERNEL(caldyn_vert)
28
29  DO iq=1,nqdyn
30    FORALL_CELLS('2', 'llm')
31      ON_PRIMAL
32        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL) 
33      END_BLOCK
34    END_BLOCK
35    FORALL_CELLS('1', 'llm-1')
36      ON_PRIMAL
37        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL)) 
38      END_BLOCK
39    END_BLOCK
40  END DO
41
42  IF(caldyn_vert_variant == caldyn_vert_cons) THEN
43    ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
44    FORALL_CELLS('2','llm')
45      ON_EDGES
46        wwuu(EDGE) = .25*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)+u(DOWN(EDGE))) ! Fu
47      END_BLOCK
48    END_BLOCK
49    ! make sure wwuu is up to date
50    BARRIER
51
52    FORALL_CELLS()
53      ON_EDGES
54        dFu_deta = wwuu(UP(EDGE))-wwuu(EDGE)                                          ! d/deta (F*u)
55        dF_deta  = .5*(wflux(UP(CELL1))+wflux(UP(CELL2))-(wflux(CELL1)+wflux(CELL2))) ! d/deta(F)
56        du(EDGE) = du(EDGE) - (dFu_deta-u(EDGE)*dF_deta) / (.5*(rhodz(CELL1)+rhodz(CELL2)))  ! (F/m)du/deta
57      END_BLOCK
58    END_BLOCK
59  ELSE 
60    FORALL_CELLS('2','llm')
61      ON_EDGES
62        wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
63      END_BLOCK
64    END_BLOCK
65
66    ! make sure wwuu is up to date
67    BARRIER
68
69    FORALL_CELLS()
70      ON_EDGES
71        du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
72      END_BLOCK
73    END_BLOCK
74  END IF
75
76END_BLOCK
77
78KERNEL(gradient)
79  FORALL_CELLS_EXT()
80    ON_EDGES
81      grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
82    END_BLOCK
83  END_BLOCK
84END_BLOCK
85
86KERNEL(div)
87  FORALL_CELLS_EXT()
88    ON_PRIMAL
89      div_ij=0.
90      FORALL_EDGES
91        div_ij = div_ij + SIGN*LE_DE*u(EDGE)
92      END_BLOCK
93      divu(CELL) = div_ij / AI
94    END_BLOCK
95  END_BLOCK
96END_BLOCK
97
98KERNEL(theta)
99  IF(caldyn_eta==eta_mass) THEN ! Compute mass
100    ! FIXME : here mass_col is computed from rhodz
101    ! so that the DOFs are the same whatever caldyn_eta
102    ! in DYNAMICO mass_col is prognosed rather than rhodz
103    SEQUENCE_C1
104      PROLOGUE(0)
105        mass_col(HIDX(CELL))=0.
106      END_BLOCK
107      BODY('1,llm')
108          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)
109      END_BLOCK
110    END_BLOCK
111    FORALL_CELLS_EXT()
112      ON_PRIMAL
113        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
114        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
115        !        rhodz(CELL) = m/g
116        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
117      END_BLOCK
118    END_BLOCK
119  END IF
120  DO iq=1,nqdyn
121    FORALL_CELLS_EXT()
122      ON_PRIMAL
123        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
124      END_BLOCK
125    END_BLOCK
126  END DO
127END_BLOCK
128
129KERNEL(caldyn_fast)
130!
131  SELECT CASE(caldyn_thermo)
132  CASE(thermo_boussinesq)
133    FORALL_CELLS()
134      ON_PRIMAL
135        berni(CELL) = pk(CELL)
136        ! from now on pk contains the vertically-averaged geopotential
137        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
138      END_BLOCK
139    END_BLOCK
140  CASE(thermo_theta)
141    FORALL_CELLS()
142      ON_PRIMAL
143        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
144      END_BLOCK
145    END_BLOCK
146  CASE(thermo_entropy)
147    FORALL_CELLS()
148      ON_PRIMAL
149        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
150        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
151      END_BLOCK
152    END_BLOCK
153  CASE(thermo_variable_Cp)
154    ! thermodynamics with variable Cp
155    !           Cp(T) = Cp0 * (T/T0)^nu
156    ! =>            h = Cp(T).T/(nu+1)
157
158    FORALL_CELLS()
159      ON_PRIMAL
160        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
161        cp_ik = cpp*(pk(CELL)/Treff)**nu
162        berni(CELL) = berni(CELL) + pk(CELL)*(cp_ik/(nu+1.)-theta(CELL,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
163      END_BLOCK
164    END_BLOCK
165  CASE DEFAULT
166    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
167    STOP
168  END SELECT
169!
170  FORALL_CELLS()
171    ON_EDGES
172      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
173      du(EDGE) = du(EDGE) - SIGN*due
174      u(EDGE) = u(EDGE) + tau*du(EDGE)
175    END_BLOCK
176  END_BLOCK 
177!
178END_BLOCK
179
180KERNEL(pvort_only)
181  FORALL_CELLS_EXT()
182    ON_DUAL
183      etav = 0.d0
184      FORALL_EDGES
185         etav = etav + SIGN*u(EDGE)
186      END_BLOCK
187      hv=0.
188      FORALL_VERTICES
189        hv = hv + RIV2*rhodz(VERTEX)
190      END_BLOCK
191      qv(DUAL_CELL) = (etav + FV*AV )/(hv*AV)
192    END_BLOCK
193  END_BLOCK
194
195  FORALL_CELLS()
196    ON_EDGES
197      qu(EDGE)=0.5d0*(qv(VERTEX1)+qv(VERTEX2))
198    END_BLOCK
199  END_BLOCK
200END_BLOCK
201
202KERNEL(coriolis)
203!
204  DO iq=1,nqdyn
205    FORALL_CELLS_EXT()
206      ON_EDGES
207        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
208      END_BLOCK
209    END_BLOCK
210    FORALL_CELLS()
211      ON_PRIMAL
212        divF=0.
213        FORALL_EDGES
214          divF = divF + Ftheta(EDGE)*SIGN
215        END_BLOCK
216        dtheta_rhodz(CELL,iq) = -divF / AI
217      END_BLOCK
218    END_BLOCK
219  END DO ! iq
220!
221  FORALL_CELLS()
222    ON_PRIMAL
223      divF=0.
224      FORALL_EDGES
225        divF = divF + hflux(EDGE)*SIGN
226      END_BLOCK
227      convm(CELL) = -divF / AI
228    END_BLOCK
229  END_BLOCK
230!
231  FORALL_CELLS()
232    ON_EDGES
233      du_trisk=0.
234      FORALL_TRISK
235        du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
236      END_BLOCK
237      du(EDGE) = du(EDGE) + .5*du_trisk
238    END_BLOCK
239  END_BLOCK
240
241END_BLOCK
Note: See TracBrowser for help on using the repository browser.