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

Last change on this file since 615 was 615, checked in by dubos, 6 years ago

devel/unstructured : import Python layer from HEAT
File size: 4.8 KB
Line 
1KERNEL(caldyn_vert)
2  SEQUENCE
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
16  ! make sure wflux is up to date
17  BARRIER
18 
19  FORALL_CELLS('2','llm')
20    ON_EDGES
21      wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
22    END_BLOCK
23  END_BLOCK
24
25  ! make sure wwuu is up to date
26  BARRIER
27
28  FORALL_CELLS()
29    ON_PRIMAL
30      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL)) ! FIXME : we shoud prognose mass_col
31    END_BLOCK
32    ON_EDGES
33      du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
34    END_BLOCK
35  END_BLOCK
36
37  DO iq=1,nqdyn
38    FORALL_CELLS('2', 'llm')
39      ON_PRIMAL
40        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL) 
41      END_BLOCK
42    END_BLOCK
43    FORALL_CELLS('1', 'llm-1')
44      ON_PRIMAL
45        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL)) 
46      END_BLOCK
47    END_BLOCK
48  END DO
49END_BLOCK
50
51KERNEL(gradient)
52  FORALL_CELLS_EXT()
53    ON_EDGES
54      grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
55    END_BLOCK
56  END_BLOCK
57END_BLOCK
58
59KERNEL(div)
60  FORALL_CELLS_EXT()
61    ON_PRIMAL
62      div_ij=0.
63      FORALL_EDGES
64        div_ij = div_ij + SIGN*LE_DE*u(EDGE)
65      END_BLOCK
66      divu(CELL) = div_ij / AI
67    END_BLOCK
68  END_BLOCK
69END_BLOCK
70
71KERNEL(theta)
72  IF(caldyn_eta==eta_mass) THEN ! Compute mass
73    ! FIXME : here mass_col is computed from rhodz
74    ! so that the DOFs are the same whatever caldyn_eta
75    ! in DYNAMICO mass_col is prognosed rather than rhodz
76    SEQUENCE_EXT
77      PROLOGUE(0)
78        mass_col(HIDX(CELL))=0.
79      END_BLOCK
80      BODY('1,llm')
81          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)
82      END_BLOCK
83    END_BLOCK
84    FORALL_CELLS_EXT()
85      ON_PRIMAL
86        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
87        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
88        !        rhodz(CELL) = m/g
89        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
90      END_BLOCK
91    END_BLOCK
92  END IF
93  DO iq=1,nqdyn
94    FORALL_CELLS_EXT()
95      ON_PRIMAL
96        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
97      END_BLOCK
98    END_BLOCK
99  END DO
100END_BLOCK
101
102KERNEL(caldyn_fast)
103!
104  SELECT CASE(caldyn_thermo)
105  CASE(thermo_boussinesq)
106    FORALL_CELLS()
107      ON_PRIMAL
108        berni(CELL) = pk(CELL)
109        ! from now on pk contains the vertically-averaged geopotential
110        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
111      END_BLOCK
112    END_BLOCK
113  CASE(thermo_theta)
114    FORALL_CELLS()
115      ON_PRIMAL
116        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
117      END_BLOCK
118    END_BLOCK
119  CASE(thermo_entropy)
120    FORALL_CELLS()
121      ON_PRIMAL
122        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
123        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
124      END_BLOCK
125    END_BLOCK
126  CASE DEFAULT
127    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
128    STOP
129  END SELECT
130!
131  FORALL_CELLS()
132    ON_EDGES
133      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
134      du(EDGE) = du(EDGE) - SIGN*due
135      u(EDGE) = u(EDGE) + tau*du(EDGE)
136    END_BLOCK
137  END_BLOCK 
138!
139END_BLOCK
140
141KERNEL(pvort_only)
142  FORALL_CELLS_EXT()
143    ON_DUAL
144      etav = 0.d0
145      FORALL_EDGES
146         etav = etav + SIGN*u(EDGE)
147      END_BLOCK
148      hv=0.
149      FORALL_VERTICES
150        hv = hv + RIV2*rhodz(VERTEX)
151      END_BLOCK
152      qv(DUAL_CELL) = (etav + FV*AV )/(hv*AV)
153    END_BLOCK
154  END_BLOCK
155
156  FORALL_CELLS()
157    ON_EDGES
158      qu(EDGE)=0.5d0*(qv(VERTEX1)+qv(VERTEX2))
159    END_BLOCK
160  END_BLOCK
161END_BLOCK
162
163KERNEL(coriolis)
164!
165  DO iq=1,nqdyn
166    FORALL_CELLS_EXT()
167      ON_EDGES
168        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
169      END_BLOCK
170    END_BLOCK
171    FORALL_CELLS()
172      ON_PRIMAL
173        divF=0.
174        FORALL_EDGES
175          divF = divF + Ftheta(EDGE)*SIGN
176        END_BLOCK
177        dtheta_rhodz(CELL,iq) = -divF / AI
178      END_BLOCK
179    END_BLOCK
180  END DO ! iq
181!
182  FORALL_CELLS()
183    ON_PRIMAL
184      divF=0.
185      FORALL_EDGES
186        divF = divF + hflux(EDGE)*SIGN
187      END_BLOCK
188      convm(CELL) = -divF / AI
189    END_BLOCK
190  END_BLOCK
191!
192  FORALL_CELLS()
193    ON_EDGES
194      du_trisk=0.
195      FORALL_TRISK
196        du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
197      END_BLOCK
198      du(EDGE) = du(EDGE) + .5*du_trisk
199    END_BLOCK
200  END_BLOCK
201
202END_BLOCK
Note: See TracBrowser for help on using the repository browser.