source: codes/icosagcm/devel/src/dynamics/compute_caldyn_vert.F90 @ 977

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

devel : separate modules for caldyn_vert and caldyn_vert_NH
File size: 7.0 KB
Line 
1MODULE compute_caldyn_vert_mod
2  USE prec, ONLY : rstd
3  USE caldyn_vars_mod
4  USE grid_param
5  USE disvert_mod
6  USE omp_para
7  USE trace
8  IMPLICIT NONE
9  PRIVATE
10  SAVE
11
12  PUBLIC :: compute_caldyn_vert_manual, &
13       compute_caldyn_vert_hex
14CONTAINS
15
16#ifdef BEGIN_DYSL
17
18KERNEL(caldyn_wflux)
19  SEQUENCE_C0
20    BODY('llm-1,1,-1')
21      ! cumulate mass flux convergence from top to bottom
22      convm(CELL) = convm(CELL) + convm(UP(CELL))
23    END_BLOCK
24    EPILOGUE(1)
25      dmass_col(HIDX(CELL)) = convm(CELL)
26    END_BLOCK
27    BODY('2,llm')
28      ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
29      wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL) 
30    END_BLOCK
31  END_BLOCK
32  ! make sure wflux is up to date
33  BARRIER
34END_BLOCK
35
36KERNEL(caldyn_dmass)
37  FORALL_CELLS()
38    ON_PRIMAL
39      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL))
40    END_BLOCK
41  END_BLOCK
42END_BLOCK
43
44KERNEL(caldyn_vert)
45  DO iq=1,nqdyn
46    FORALL_CELLS('2', 'llm')
47      ON_PRIMAL
48        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL) 
49      END_BLOCK
50    END_BLOCK
51    FORALL_CELLS('1', 'llm-1')
52      ON_PRIMAL
53        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL)) 
54      END_BLOCK
55    END_BLOCK
56  END DO
57
58  IF(caldyn_vert_variant == caldyn_vert_cons) THEN
59    ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
60    FORALL_CELLS('2','llm')
61      ON_EDGES
62        wwuu(EDGE) = .25*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)+u(DOWN(EDGE))) ! Fu
63      END_BLOCK
64    END_BLOCK
65    ! make sure wwuu is up to date
66    BARRIER
67
68    FORALL_CELLS()
69      ON_EDGES
70        dFu_deta = wwuu(UP(EDGE))-wwuu(EDGE)                                          ! d/deta (F*u)
71        dF_deta  = .5*(wflux(UP(CELL1))+wflux(UP(CELL2))-(wflux(CELL1)+wflux(CELL2))) ! d/deta(F)
72        du(EDGE) = du(EDGE) - (dFu_deta-u(EDGE)*dF_deta) / (.5*(rhodz(CELL1)+rhodz(CELL2)))  ! (F/m)du/deta
73      END_BLOCK
74    END_BLOCK
75  ELSE 
76    FORALL_CELLS('2','llm')
77      ON_EDGES
78        wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
79      END_BLOCK
80    END_BLOCK
81
82    ! make sure wwuu is up to date
83    BARRIER
84
85    FORALL_CELLS()
86      ON_EDGES
87        du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
88      END_BLOCK
89    END_BLOCK
90  END IF
91END_BLOCK
92
93#endif END_DYSL
94
95  SUBROUTINE compute_caldyn_vert_hex(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
96    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
97    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
98    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
99    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
100    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
101    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
102    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
103    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
104    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
105
106    ! temporary variable   
107    INTEGER :: i,j,ij,l,iq
108    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
109    INTEGER    :: ij_omp_begin, ij_omp_end
110
111    CALL trace_start("compute_caldyn_vert")
112
113!$OMP BARRIER
114
115    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
116
117#define mass_bl(ij,l) bp(l)
118#define dmass_col(ij) dps(ij)
119#include "../kernels_hex/caldyn_wflux.k90"
120#include "../kernels_hex/caldyn_vert.k90"
121#undef mass_bl
122#undef dmass_col
123
124    CALL trace_end("compute_caldyn_vert")
125
126  END SUBROUTINE compute_caldyn_vert_hex
127
128  SUBROUTINE compute_caldyn_vert_manual(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
129    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
130    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
131    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
132    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
133    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
134    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
135    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
136    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
137    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
138
139    ! temporary variable   
140    INTEGER :: i,j,ij,l,iq
141    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
142    INTEGER    :: ij_omp_begin, ij_omp_end
143
144    CALL trace_start("compute_caldyn_vert")
145
146!$OMP BARRIER
147
148    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
149
150!!! cumulate mass flux convergence from top to bottom
151    DO  l = llm-1, 1, -1
152       !DIR$ SIMD
153       DO ij=ij_omp_begin,ij_omp_end         
154          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
155       ENDDO
156    ENDDO
157    !  ENDIF
158
159    !$OMP BARRIER
160    ! FLUSH on convm
161    ! compute dmass_col
162    IF (is_omp_first_level) THEN
163       !DIR$ SIMD
164       DO ij=ij_begin,ij_end         
165          ! dps/dt = -int(div flux)dz
166          dps(ij) = convm(ij,1)
167       ENDDO
168    ENDIF
169!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
170    DO l=ll_beginp1,ll_end
171       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
172       !DIR$ SIMD
173       DO ij=ij_begin,ij_end         
174          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
175          ! => w>0 for upward transport
176          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
177       ENDDO
178    ENDDO
179
180    !--> flush wflux 
181    !$OMP BARRIER
182
183    DO iq=1,nqdyn
184       DO l=ll_begin,ll_endm1
185       !DIR$ SIMD
186          DO ij=ij_begin,ij_end         
187             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  -   0.5 * &
188                  ( wflux(ij,l+1) * (theta(ij,l,iq) + theta(ij,l+1,iq)))
189          END DO
190       END DO
191       DO l=ll_beginp1,ll_end
192          !DIR$ SIMD
193          DO ij=ij_begin,ij_end         
194             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  +   0.5 * &
195                  ( wflux(ij,l) * (theta(ij,l-1,iq) + theta(ij,l,iq) ) )
196          END DO
197
198       END DO
199    END DO
200
201    ! Compute vertical transport
202    DO l=ll_beginp1,ll_end
203       !DIR$ SIMD
204       DO ij=ij_begin,ij_end         
205          wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
206          wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
207          wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
208       ENDDO
209    ENDDO
210
211    !--> flush wwuu 
212    !$OMP BARRIER
213
214    ! Add vertical transport to du
215    DO l=ll_begin,ll_end
216       !DIR$ SIMD
217       DO ij=ij_begin,ij_end         
218          du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
219          du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
220          du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
221       ENDDO
222    ENDDO
223
224    CALL trace_end("compute_caldyn_vert")
225
226  END SUBROUTINE compute_caldyn_vert_manual
227
228END MODULE compute_caldyn_vert_mod
Note: See TracBrowser for help on using the repository browser.