source: codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90 @ 920

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

devel : separate module for compute_geopot
File size: 6.5 KB
Line 
1MODULE caldyn_kernels_base_mod
2  USE icosa
3  USE transfert_mod
4  USE disvert_mod
5  USE caldyn_vars_mod
6  USE omp_para
7  USE trace
8  IMPLICIT NONE
9  PRIVATE
10  SAVE
11
12  PUBLIC :: compute_caldyn_vert, compute_caldyn_vert_nh
13
14CONTAINS
15
16
17  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
18    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
19    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
20    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
21    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
22    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
23    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
24    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
25    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
26    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
27
28    ! temporary variable   
29    INTEGER :: i,j,ij,l,iq
30    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
31    INTEGER    :: ij_omp_begin, ij_omp_end
32
33    CALL trace_start("compute_caldyn_vert")
34
35!$OMP BARRIER
36
37    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
38
39    IF(dysl_caldyn_vert) THEN
40#define mass_bl(ij,l) bp(l)
41#define dmass_col(ij) dps(ij)
42#include "../kernels_hex/caldyn_wflux.k90"
43#include "../kernels_hex/caldyn_vert.k90"
44#undef mass_bl
45#undef dmass_col
46    ELSE
47
48!!! cumulate mass flux convergence from top to bottom
49    DO  l = llm-1, 1, -1
50       !DIR$ SIMD
51       DO ij=ij_omp_begin,ij_omp_end         
52          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
53       ENDDO
54    ENDDO
55    !  ENDIF
56
57    !$OMP BARRIER
58    ! FLUSH on convm
59    ! compute dmass_col
60    IF (is_omp_first_level) THEN
61       !DIR$ SIMD
62       DO ij=ij_begin,ij_end         
63          ! dps/dt = -int(div flux)dz
64          dps(ij) = convm(ij,1)
65       ENDDO
66    ENDIF
67
68!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
69    DO l=ll_beginp1,ll_end
70       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
71       !DIR$ SIMD
72       DO ij=ij_begin,ij_end         
73          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
74          ! => w>0 for upward transport
75          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
76       ENDDO
77    ENDDO
78
79    !--> flush wflux 
80    !$OMP BARRIER
81
82    DO iq=1,nqdyn
83       DO l=ll_begin,ll_endm1
84       !DIR$ SIMD
85          DO ij=ij_begin,ij_end         
86             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  -   0.5 * &
87                  ( wflux(ij,l+1) * (theta(ij,l,iq) + theta(ij,l+1,iq)))
88          END DO
89       END DO
90       DO l=ll_beginp1,ll_end
91          !DIR$ SIMD
92          DO ij=ij_begin,ij_end         
93             dtheta_rhodz(ij, l, iq) = dtheta_rhodz(ij, l, iq)  +   0.5 * &
94                  ( wflux(ij,l) * (theta(ij,l-1,iq) + theta(ij,l,iq) ) )
95          END DO
96       END DO
97    END DO
98
99    ! Compute vertical transport
100    DO l=ll_beginp1,ll_end
101       !DIR$ SIMD
102       DO ij=ij_begin,ij_end         
103          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))
104          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))
105          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))
106       ENDDO
107    ENDDO
108
109    !--> flush wwuu 
110    !$OMP BARRIER
111
112    ! Add vertical transport to du
113    DO l=ll_begin,ll_end
114       !DIR$ SIMD
115       DO ij=ij_begin,ij_end         
116          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))
117          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))
118          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))
119       ENDDO
120    ENDDO
121
122    END IF ! dysl
123
124    CALL trace_end("compute_caldyn_vert")
125
126  END SUBROUTINE compute_caldyn_vert
127
128  SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
129    REAL(rstd),INTENT(IN) :: mass(iim*jjm,llm)
130    REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1)
131    REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1)
132    REAL(rstd),INTENT(IN) :: wflux(iim*jjm,llm+1)
133    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
134    REAL(rstd),INTENT(INOUT) :: dPhi(iim*jjm,llm+1)
135    REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
136    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
137    ! local arrays
138    REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers
139    REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers
140    ! indices and temporary values
141    INTEGER    :: ij, l
142    REAL(rstd) :: wflux_ij, w_ij
143
144    CALL trace_start("compute_caldyn_vert_nh")
145
146    IF(dysl) THEN
147!$OMP BARRIER
148#include "../kernels_hex/caldyn_vert_NH.k90"
149!$OMP BARRIER
150    ELSE
151#define ETA_DOT(ij) eta_dot(ij,1)
152#define WCOV(ij) wcov(ij,1)
153
154    DO l=ll_begin,ll_end
155       ! compute the local arrays
156       !DIR$ SIMD
157       DO ij=ij_begin_ext,ij_end_ext
158          wflux_ij = .5*(wflux(ij,l)+wflux(ij,l+1))
159          w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l)
160          W_etadot(ij,l) = wflux_ij*w_ij
161          ETA_DOT(ij) = wflux_ij / mass(ij,l)
162          WCOV(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
163       ENDDO
164       ! add NH term to du
165      !DIR$ SIMD
166      DO ij=ij_begin,ij_end
167          du(ij+u_right,l) = du(ij+u_right,l) &
168               - .5*(WCOV(ij+t_right)+WCOV(ij)) &
169               *ne_right*(ETA_DOT(ij+t_right)-ETA_DOT(ij))
170          du(ij+u_lup,l) = du(ij+u_lup,l) &
171               - .5*(WCOV(ij+t_lup)+WCOV(ij)) &
172               *ne_lup*(ETA_DOT(ij+t_lup)-ETA_DOT(ij))
173          du(ij+u_ldown,l) = du(ij+u_ldown,l) &
174               - .5*(WCOV(ij+t_ldown)+WCOV(ij)) &
175               *ne_ldown*(ETA_DOT(ij+t_ldown)-ETA_DOT(ij))
176       END DO
177    ENDDO
178    ! add NH terms to dW, dPhi
179    ! FIXME : TODO top and bottom
180    DO l=ll_beginp1,ll_end ! inner interfaces only
181       !DIR$ SIMD
182       DO ij=ij_begin,ij_end
183          dPhi(ij,l) = dPhi(ij,l) - wflux(ij,l) &
184               * (geopot(ij,l+1)-geopot(ij,l-1))/(mass(ij,l-1)+mass(ij,l))
185       END DO
186    END DO
187    DO l=ll_begin,ll_end
188       !DIR$ SIMD
189       DO ij=ij_begin,ij_end
190          dW(ij,l+1) = dW(ij,l+1) + W_etadot(ij,l) ! update inner+top interfaces
191          dW(ij,l)   = dW(ij,l)   - W_etadot(ij,l) ! update bottom+inner interfaces
192       END DO
193    END DO
194
195#undef ETA_DOT
196#undef WCOV
197
198    END IF ! dysl
199    CALL trace_end("compute_caldyn_vert_nh")
200
201  END SUBROUTINE compute_caldyn_vert_NH
202
203END MODULE caldyn_kernels_base_mod
Note: See TracBrowser for help on using the repository browser.