source: codes/icosagcm/devel/src/diagnostics/kinetic.f90 @ 728

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

devel : diagnose consistent kinetic energy using new req_z1_scal (r711)
File size: 5.1 KB
Line 
1MODULE kinetic_mod
2IMPLICIT NONE
3PRIVATE
4
5PUBLIC :: kinetic, kinetic_new, gradient
6
7CONTAINS
8 
9  SUBROUTINE kinetic(f_ue,f_Ki)
10  USE icosa
11  IMPLICIT NONE
12    TYPE(t_field), POINTER :: f_ue(:)
13    TYPE(t_field), POINTER :: f_Ki(:)
14 
15    REAL(rstd), POINTER :: ue(:,:)
16    REAL(rstd), POINTER :: Ki(:,:)
17    INTEGER :: ind
18 
19    CALL transfert_request(f_ue,req_e1_vect)
20
21    DO ind=1,ndomain
22      IF (.NOT. assigned_domain(ind)) CYCLE
23      CALL swap_dimensions(ind)
24      CALL swap_geometry(ind)
25      ue=f_ue(ind)
26      Ki=f_Ki(ind)
27      CALL compute_kinetic_trisk(ue, Ki)
28    ENDDO 
29  END SUBROUTINE kinetic
30 
31  SUBROUTINE kinetic_new(f_ue,f_Kv,f_Ki)
32    USE icosa
33    IMPLICIT NONE
34    TYPE(t_field), POINTER :: f_ue(:)
35    TYPE(t_field), POINTER :: f_Kv(:)
36    TYPE(t_field), POINTER :: f_Ki(:)
37   
38    REAL(rstd), POINTER :: ue(:,:)
39    REAL(rstd), POINTER :: Kv(:,:)
40    REAL(rstd), POINTER :: Ki(:,:)
41    INTEGER :: ind
42   
43    CALL transfert_request(f_ue,req_e1_vect)
44
45    DO ind=1,ndomain
46      IF (.NOT. assigned_domain(ind)) CYCLE
47      CALL swap_dimensions(ind)
48      CALL swap_geometry(ind)
49      ue=f_ue(ind)
50      Kv=f_Kv(ind)
51      CALL compute_kv(ue, Kv)
52    ENDDO
53   
54    CALL transfert_request(f_Kv,req_z1_scal)
55
56    DO ind=1,ndomain
57       IF (.NOT. assigned_domain(ind)) CYCLE
58       CALL swap_dimensions(ind)
59       CALL swap_geometry(ind)
60       Kv=f_Kv(ind)
61       Ki=f_Ki(ind)
62       CALL compute_Ki_from_Kv(Kv, Ki)
63    ENDDO
64  END SUBROUTINE kinetic_new
65 
66  SUBROUTINE compute_kinetic_trisk(ue, Ki)
67  USE icosa
68  USE omp_para
69  IMPLICIT NONE
70    REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
71    REAL(rstd),INTENT(OUT) :: Ki(iim*jjm,llm)
72    INTEGER :: i,j,ij,l
73   
74    DO l=ll_begin,ll_end
75      DO j=jj_begin,jj_end
76        DO i=ii_begin,ii_end
77          ij=(j-1)*iim+i
78
79          Ki(ij,l)=1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*ue(ij+u_right,l)**2 +  &
80                               le(ij+u_rup)*de(ij+u_rup)*ue(ij+u_rup,l)**2 +        &
81                               le(ij+u_lup)*de(ij+u_lup)*ue(ij+u_lup,l)**2 +        &
82                               le(ij+u_left)*de(ij+u_left)*ue(ij+u_left,l)**2 +     &
83                               le(ij+u_ldown)*de(ij+u_ldown)*ue(ij+u_ldown,l)**2 +  &
84                               le(ij+u_rdown)*de(ij+u_rdown)*ue(ij+u_rdown,l)**2 ) 
85       
86        ENDDO
87      ENDDO
88    ENDDO
89  END SUBROUTINE compute_kinetic_trisk
90 
91  SUBROUTINE compute_kv(ue, Kv)
92  USE icosa
93  USE omp_para
94  IMPLICIT NONE
95    REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
96    REAL(rstd),INTENT(OUT) :: Kv(2*iim*jjm,llm)
97    INTEGER :: ij,l, u_up, u_down
98
99    u_up    = t_lup + u_right
100    u_down  = t_rdown + u_left
101   
102    DO l=ll_begin,ll_end
103      Kv(:,l)=0.
104      DO ij=ij_begin,ij_end
105          Kv(ij+z_up,l) = (radius**2/Av(ij+z_up))*(                         &
106                               S1(ij,vup)*ue(ij+u_rup,l)**2 +        &
107                               S2(ij,vup)*ue(ij+u_lup,l)**2 +        &
108                               S2(ij+t_lup,vrdown)*ue(ij+u_up,l)**2)
109
110          Kv(ij+z_down,l) = (radius**2/Av(ij+z_down))*(                         &
111                               S1(ij,vdown)*ue(ij+u_ldown,l)**2 +       &
112                               S2(ij,vdown)*ue(ij+u_rdown,l)**2 +       &
113                               S2(ij+t_rdown,vlup)*ue(ij+u_down,l)**2 )
114       ENDDO
115    ENDDO
116  END SUBROUTINE compute_kv
117 
118  SUBROUTINE compute_Ki_from_Kv(Kv, Ki)
119    USE icosa
120    USE omp_para
121    IMPLICIT NONE
122    REAL(rstd),INTENT(OUT):: Ki(iim*jjm,llm)
123    REAL(rstd), INTENT(IN) :: Kv(2*iim*jjm,llm)
124    INTEGER :: ij,l
125   
126    DO l=ll_begin,ll_end
127       DO ij=ij_begin,ij_end
128          Ki(ij,l) = Riv(ij,vup)*Kv(ij+z_up,l) + &
129               Riv(ij,vlup)  *Kv(ij+z_lup,l) + &
130               Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
131               Riv(ij,vdown) *Kv(ij+z_down,l) + &
132               Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
133               Riv(ij,vrup)  *Kv(ij+z_rup,l)
134       END DO
135    END DO
136  END SUBROUTINE compute_Ki_from_Kv
137 
138  SUBROUTINE gradient(f_berni, f_du)
139  USE icosa
140  IMPLICIT NONE
141    TYPE(t_field), POINTER :: f_berni(:)
142    TYPE(t_field), POINTER :: f_du(:)
143 
144    REAL(rstd), POINTER :: du(:,:)
145    REAL(rstd), POINTER :: berni(:,:)
146    INTEGER :: ind
147 
148    CALL transfert_request(f_du,req_e1_vect)
149    CALL transfert_request(f_du,req_e1_vect)
150
151    DO ind=1,ndomain
152      IF (.NOT. assigned_domain(ind)) CYCLE
153      CALL swap_dimensions(ind)
154      CALL swap_geometry(ind)
155      berni=f_berni(ind)
156      du=f_du(ind)
157      CALL compute_grad(berni, du)
158    ENDDO
159 
160  END SUBROUTINE gradient
161 
162  SUBROUTINE compute_grad(berni, du)
163  USE icosa
164  USE omp_para
165  IMPLICIT NONE
166    REAL(rstd),INTENT(IN) :: berni(iim*jjm,llm)
167    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
168    INTEGER :: i,j,ij,l
169   
170    DO l=ll_begin,ll_end
171      DO j=jj_begin,jj_end
172        DO i=ii_begin,ii_end
173          ij=(j-1)*iim+i
174            du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l))/de(ij+u_right)
175            du(ij+u_lup,l)   = ne_lup*(berni(ij,l)-berni(ij+t_lup,l))/de(ij+u_right)
176            du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l))/de(ij+u_right)
177         ENDDO
178      ENDDO
179    ENDDO
180
181  END SUBROUTINE compute_grad
182
183
184END MODULE kinetic_mod
Note: See TracBrowser for help on using the repository browser.