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

Last change on this file since 876 was 876, checked in by jisesh, 5 years ago

devel: moved DYSL into compute_caldyn_slow_NH.F90 and compute_caldyn_Coriolis.F90
File size: 3.4 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
Note: See TracBrowser for help on using the repository browser.