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

Last change on this file since 548 was 548, checked in by dubos, 7 years ago

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