source: codes/icosagcm/trunk/src/initial/etat0_academic.f90

Last change on this file was 1046, checked in by ymipsl, 4 years ago

Introduce modification from A. Durocher github to make held&suarez testcase working on GPU

YM & AD
File size: 8.7 KB
Line 
1MODULE etat0_academic_mod
2  USE icosa
3  IMPLICIT NONE
4
5  PRIVATE
6 
7  PUBLIC :: etat0
8 
9CONTAINS
10 
11  SUBROUTINE test_etat0_academic
12  USE icosa
13  USE kinetic_mod
14    TYPE(t_field),POINTER,SAVE :: f_ps(:)
15    TYPE(t_field),POINTER,SAVE :: f_phis(:)
16    TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
17    TYPE(t_field),POINTER,SAVE :: f_u(:)
18    TYPE(t_field),POINTER,SAVE :: f_q(:)
19    TYPE(t_field),POINTER,SAVE :: f_Ki(:)
20    TYPE(t_field),POINTER,SAVE :: f_temp(:)
21       
22    CALL allocate_field(f_ps,field_t,type_real)
23    CALL allocate_field(f_phis,field_t,type_real)
24    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
25    CALL allocate_field(f_u,field_u,type_real,llm)
26    CALL allocate_field(f_Ki,field_t,type_real,llm)
27    CALL allocate_field(f_temp,field_t,type_real)
28   
29    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
30
31    CALL kinetic(f_u,f_Ki)
32
33    CALL writefield('ps',f_ps)
34    CALL writefield('theta',f_theta_rhodz)
35
36    CALL writefield('Ki',f_Ki)
37   
38  END SUBROUTINE test_etat0_academic
39   
40  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
41  USE icosa
42    TYPE(t_field),POINTER :: f_ps(:)
43    TYPE(t_field),POINTER :: f_phis(:)
44    TYPE(t_field),POINTER :: f_theta_rhodz(:)
45    TYPE(t_field),POINTER :: f_u(:)
46    TYPE(t_field),POINTER :: f_q(:)
47 
48    REAL(rstd),POINTER :: ps(:)
49    REAL(rstd),POINTER :: phis(:)
50    REAL(rstd),POINTER :: theta_rhodz(:,:)
51    REAL(rstd),POINTER :: u(:,:)
52    INTEGER :: ind
53
54    PRINT *, 'etat0_academic needs an upgrade for 4D theta'
55    PRINT *, 'STOP in etat0_academic.f90/etat0'
56    STOP
57
58    DO ind=1,ndomain
59      IF (.NOT. assigned_domain(ind)) CYCLE
60      CALL swap_dimensions(ind)
61      CALL swap_geometry(ind)
62      ps=f_ps(ind)
63      phis=f_phis(ind)
64      theta_rhodz=f_theta_rhodz(ind)
65      u=f_u(ind)
66      CALL compute_etat0_academic(ps, phis, theta_rhodz, u)
67    ENDDO
68
69  END SUBROUTINE etat0
70 
71  SUBROUTINE compute_etat0_academic(ps, phis, theta_rhodz, u)
72  USE icosa
73  USE disvert_mod
74  USE pression_mod
75  USE exner_mod
76  USE geopotential_mod
77  USE theta2theta_rhodz_mod
78  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
79  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
80  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
81  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
82 
83  INTEGER :: i,j,l,ij
84  REAL(rstd) :: r
85  REAL(rstd) :: theta(iim*jjm,llm)
86  REAL(rstd) :: zsig
87  INTEGER    :: lsup
88  REAL(rstd) :: ddsin
89  REAL(rstd) :: thetarappel
90  REAL(rstd) :: lat
91  REAL(rstd) :: p(iim*jjm,llm+1)
92  REAL(rstd) :: phi(iim*jjm,llm)
93  REAL(rstd) :: x 
94  REAL(rstd) :: fact(3*iim*jjm)
95  REAL(rstd) :: ut(3*iim*jjm,llm)
96   
97    DO l=1,llm
98       zsig=ap(l)/preff+bp(l)
99       IF (zsig.gt.0.3) THEN
100         lsup=l
101         thetarappel=1./8.*(-log(zsig)-.5)
102         DO j=jj_begin-1,jj_end+1
103           DO i=ii_begin-1,ii_end+1
104             ij=(j-1)*iim+i
105             ddsin=sin(lat_i(ij))-sin(pi/20.)
106             theta(ij,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+thetarappel)
107            ENDDO
108         ENDDO
109       ELSE
110!   Choix isotherme au-dessus de 300 mbar
111         DO j=jj_begin-1,jj_end+1
112           DO i=ii_begin-1,ii_end+1
113             ij=(j-1)*iim+i
114             theta(ij,l)=theta(ij,lsup)*(0.3/zsig)**kappa
115           ENDDO
116         ENDDO
117       ENDIF
118     ENDDO
119     
120     
121     DO j=jj_begin-1,jj_end+1
122       DO i=ii_begin-1,ii_end+1
123         ij=(j-1)*iim+i
124         ps(ij)=1.e5
125         phis(ij)=0
126       ENDDO
127     ENDDO
128     
129
130!$OMP BARRIER
131    CALL compute_pression(ps,p,1)     
132!$OMP BARRIER
133    CALL compute_geopotential(phis,ps,theta,phi,1)
134
135
136     DO j=jj_begin-1,jj_end+1
137       DO i=ii_begin-1,ii_end+1
138         ij=(j-1)*iim+i
139
140         lat=lat_e(ij+u_right)
141         IF (abs(sin(lat))<1.e-4) lat=1e-4
142         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
143         fact(ij+u_right)=(1.-x)/(2.*omega*sin(lat))
144
145         lat=lat_e(ij+u_lup)
146         IF (abs(sin(lat))<1.e-4) lat=1e-4
147         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
148         fact(ij+u_lup)=(1.-x)/(2.*omega*sin(lat))
149
150         lat=lat_e(ij+u_ldown)
151         IF (abs(sin(lat))<1.e-4) lat=1e-4
152         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x
153         fact(ij+u_ldown)=(1.-x)/(2.*omega*sin(lat))
154       ENDDO
155     ENDDO
156
157 ! gradient ==> tangential component     
158     DO l=1,llm   
159       DO j=jj_begin-1,jj_end+1
160         DO i=ii_begin-1,ii_end+1
161           ij=(j-1)*iim+i
162           ut(ij+u_right,l)=1/de(ij+u_right)*(ne(ij,right)*phi(ij,l)+ ne(ij+t_right,left)*phi(ij+t_right,l) )
163           ut(ij+u_lup,l)=1/de(ij+u_lup)*(ne(ij,lup)*phi(ij,l)+ ne(ij+t_lup,rdown)*phi(ij+t_lup,l) )
164           ut(ij+u_ldown,l)=1/de(ij+u_ldown)*(ne(ij,ldown)*phi(ij,l)+ ne(ij+t_ldown,rup)*phi(ij+t_ldown,l) )
165         ENDDO
166       ENDDO
167     ENDDO
168     
169! attention au signe ??
170! reconstruction of perpendicular component             
171     DO l=1,llm   
172       DO j=jj_begin,jj_end
173         DO i=ii_begin,ii_end
174           ij=(j-1)*iim+i
175           u(ij+u_right,l) =  -fact(ij+u_right)/de(ij+u_right) *                                        & 
176                              ( wee(ij+u_right,1,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                        &
177                                wee(ij+u_right,2,1)*le(ij+u_lup)*ut(ij+u_lup,l)+                        &
178                                wee(ij+u_right,3,1)*le(ij+u_left)*ut(ij+u_left,l)+                      &
179                                wee(ij+u_right,4,1)*le(ij+u_ldown)*ut(ij+u_ldown,l)+                    &
180                                wee(ij+u_right,5,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                    & 
181                                wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ut(ij+t_right+u_ldown,l)+    &
182                                wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ut(ij+t_right+u_rdown,l)+    &
183                                wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ut(ij+t_right+u_right,l)+    &
184                                wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ut(ij+t_right+u_rup,l)+        &
185                                wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ut(ij+t_right+u_lup,l) )
186                               
187           u(ij+u_lup,l) = -fact(ij+u_lup)/de(ij+u_lup) *                                        & 
188                       ( wee(ij+u_lup,1,1)*le(ij+u_left)*ut(ij+u_left,l)+                        &
189                         wee(ij+u_lup,2,1)*le(ij+u_ldown)*ut(ij+u_ldown,l)+                      &
190                         wee(ij+u_lup,3,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                      &
191                         wee(ij+u_lup,4,1)*le(ij+u_right)*ut(ij+u_right,l)+                      &
192                         wee(ij+u_lup,5,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                          & 
193                         wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ut(ij+t_lup+u_right,l)+          &
194                         wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ut(ij+t_lup+u_rup,l)+              &
195                         wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ut(ij+t_lup+u_lup,l)+              &
196                         wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ut(ij+t_lup+u_left,l)+            &
197                         wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ut(ij+t_lup+u_ldown,l) )
198
199
200           u(ij+u_ldown,l) = -fact(ij+u_ldown)/de(ij+u_ldown) *                                  & 
201                       ( wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ut(ij+u_rdown,l)+                    &
202                         wee(ij+u_ldown,2,1)*le(ij+u_right)*ut(ij+u_right,l)+                    &
203                         wee(ij+u_ldown,3,1)*le(ij+u_rup)*ut(ij+u_rup,l)+                        &
204                         wee(ij+u_ldown,4,1)*le(ij+u_lup)*ut(ij+u_lup,l)+                        &
205                         wee(ij+u_ldown,5,1)*le(ij+u_left)*ut(ij+u_left,l)+                      & 
206                         wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ut(ij+t_ldown+u_lup,l)+        &
207                         wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ut(ij+t_ldown+u_left,l)+      &
208                         wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ut(ij+t_ldown+u_ldown,l)+    &
209                         wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ut(ij+t_ldown+u_rdown,l)+    &
210                         wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ut(ij+t_ldown+u_right,l) )
211                               
212!           u(ij+u_right,l)=fact(ij+u_right)*ut(ij+u_right,l)
213!           u(ij+u_lup,l)=fact(ij+u_lup)*ut(ij+u_lup,l)
214!           u(ij+u_ldown,l)=fact(ij+u_ldown)*ut(ij+u_ldown,l)
215!            ps(ij)=phi(ij,1)
216         ENDDO
217       ENDDO
218     ENDDO
219   
220     DO l=1,llm
221      DO j=jj_begin,jj_end
222        DO i=ii_begin,ii_end
223          ij=(j-1)*iim+i
224          CALL RANDOM_NUMBER(r) ;r=0
225          theta(ij,l)=theta(ij,l)*(1-0.0005*r)
226       ENDDO
227     ENDDO
228    ENDDO
229
230    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1, ondevice=.false.)
231       
232
233  END SUBROUTINE compute_etat0_academic
234 
235END MODULE etat0_academic_mod
Note: See TracBrowser for help on using the repository browser.