source: codes/icosagcm/trunk/src/etat0_heldsz.f90 @ 156

Last change on this file since 156 was 149, checked in by sdubey, 11 years ago
Added few new routines to read NC files and compute diagnostics to r145.
Few routines of dry physics including radiation module, surface process and convective adjustment in new routine phyparam.f90. dynetat to read start files for dynamics. check_conserve routine to compute conservation of quatities like mass, energy etc.etat0_heldsz.f90 for held-suarez test case initial conditions. new Key time_style=lmd or dcmip to use day_step, ndays like in LMDZ
File size: 7.1 KB
Line 
1    MODULE etat0_heldsz_mod
2         USE icosa
3      IMPLICIT NONE
4  REAL(rstd),ALLOCATABLE::knewt_t(:),kfrict(:)
5  REAL(rstd)::knewt_g
6  TYPE(t_field),POINTER :: f_tetarappel(:)
7  TYPE(t_field),POINTER :: f_clat(:)
8
9CONTAINS
10 
11  SUBROUTINE test_etat0_heldsz
12  USE icosa
13  USE kinetic_mod
14  IMPLICIT NONE
15    TYPE(t_field),POINTER :: f_ps(:)
16    TYPE(t_field),POINTER :: f_phis(:)
17    TYPE(t_field),POINTER :: f_theta_rhodz(:)
18    TYPE(t_field),POINTER :: f_u(:)
19    TYPE(t_field),POINTER :: f_q(:)
20    TYPE(t_field),POINTER :: f_Ki(:)
21 
22    REAL(rstd),POINTER :: Ki(:,:)
23    INTEGER :: ind
24   
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   
32    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
33    CALL kinetic(f_u,f_Ki)
34
35    CALL writefield('ps',f_ps)
36    CALL writefield('theta',f_theta_rhodz)
37  END SUBROUTINE test_etat0_heldsz
38   
39   
40  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
41  USE icosa
42  IMPLICIT NONE
43    TYPE(t_field),POINTER :: f_ps(:)
44    TYPE(t_field),POINTER :: f_phis(:)
45    TYPE(t_field),POINTER :: f_theta_rhodz(:)
46    TYPE(t_field),POINTER :: f_u(:)
47    TYPE(t_field),POINTER :: f_q(:)
48 
49    REAL(rstd),POINTER :: ps(:)
50    REAL(rstd),POINTER :: phis(:)
51    REAL(rstd),POINTER :: theta_rhodz(:,:)
52    REAL(rstd),POINTER :: u(:,:)
53    REAL(rstd),POINTER :: q(:,:,:)
54    INTEGER :: ind
55    REAL(rstd),POINTER::clat(:) 
56    REAL(rstd),POINTER::tetarappel(:,:) 
57   
58    CALL allocate_field(f_tetarappel,field_t,type_real,llm)
59    CALL allocate_field(f_clat,field_t,type_real)
60    ALLOCATE(knewt_t(llm)); ALLOCATE( kfrict(llm)) 
61
62    DO ind=1,ndomain
63      CALL swap_dimensions(ind)
64      CALL swap_geometry(ind)
65      ps=f_ps(ind)
66      phis=f_phis(ind)
67      theta_rhodz=f_theta_rhodz(ind)
68      tetarappel=f_tetarappel(ind) 
69      u=f_u(ind)
70      q=f_q(ind)
71      q=1e2
72      clat=f_clat(ind) 
73      CALL compute_etat0_heldsz(ps, phis, theta_rhodz, u,clat,tetarappel)
74    ENDDO
75  END SUBROUTINE etat0
76
77  SUBROUTINE compute_etat0_heldsz(ps, phis, theta_rhodz, u,clat,tetarappel)
78  USE icosa
79  USE disvert_mod
80  USE pression_mod
81  USE exner_mod
82  USE geopotential_mod
83  USE theta2theta_rhodz_mod
84  IMPLICIT NONE 
85  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
86  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
87  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
88  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
89  REAL(rstd),INTENT(OUT) :: clat(iim*jjm)
90  REAL(rstd),INTENT(OUT) :: tetarappel(iim*jjm,llm) 
91 
92  INTEGER :: i,j,l,ij
93  REAL(rstd) :: r
94  REAL(rstd) :: theta(iim*jjm,llm)
95  REAL(rstd) :: zsig
96  INTEGER    :: lsup
97  REAL(rstd) :: ddsin
98  REAL(rstd) :: lon,lat
99  REAL(rstd) :: p(iim*jjm,llm+1)
100  REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm)
101  REAL(rstd) :: delta
102  REAL(rstd) :: pks(iim*jjm),pk(iim*jjm,llm)
103  REAL(rstd) :: phi(iim*jjm,llm)
104  REAL(rstd) :: x 
105  REAL(rstd) :: fact(3*iim*jjm)
106  REAL(rstd) :: ut(3*iim*jjm,llm)
107
108  REAL(rstd) :: teta0,ttp,delt_y,delt_z,eps
109  REAL(rstd) :: k_f,k_c_a,k_c_s
110  REAL(rstd) :: zz,ran1
111  REAL(rstd) :: tetastrat,tetajl(iim*jjm,llm) 
112  REAL(rstd) :: slat(iim*jjm) 
113!-------------choces of parametes and get it 
114     k_f=1.                !friction
115     CALL getin('k_j',k_f)
116     k_f=1./(daysec*k_f)
117     k_c_s=4.  !cooling surface
118     CALL getin('k_c_s',k_c_s)
119     k_c_s=1./(daysec*k_c_s)
120     k_c_a=40. !cooling free atm
121     CALL getin('k_c_a',k_c_a)
122     k_c_a=1./(daysec*k_c_a)
123     ! Constants for Teta equilibrium profile
124     teta0=315.     ! mean Teta (S.H. 315K)
125     CALL getin('teta0',teta0)
126     ttp=200.       ! Tropopause temperature (S.H. 200K)
127     CALL getin('ttp',ttp)
128     eps=0.         ! Deviation to N-S symmetry(~0-20K)
129     CALL getin('eps',eps)
130     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
131     CALL getin('delt_y',delt_y)
132     delt_z=10.     ! Vertical Gradient (S.H. 10K)
133     CALL getin('delt_z',delt_z)
134!-----------------------------------------------------------
135    knewt_g=k_c_a 
136    DO l=1,llm
137       zsig=ap(l)/preff+bp(l)
138       knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 
139       kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 
140    ENDDO
141         DO j=jj_begin-1,jj_end+1
142           DO i=ii_begin-1,ii_end+1
143             ij=(j-1)*iim+i
144             CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
145             clat(ij)=cos(lat) 
146             slat(ij)=sin(lat) 
147            ENDDO
148         ENDDO
149
150    DO l=1,llm
151       zsig=ap(l)/preff+bp(l)
152       tetastrat=ttp*zsig**(-kappa)
153         DO j=jj_begin-1,jj_end+1
154           DO i=ii_begin-1,ii_end+1
155             ij=(j-1)*iim+i
156             ddsin=slat(ij) 
157             tetajl(ij,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
158                -delt_z*(1.-ddsin*ddsin)*log(zsig)
159             tetajl(ij,l)=MAX(tetajl(ij,l),tetastrat) 
160             tetarappel(ij,l)=tetajl(ij,l) 
161            ENDDO
162          ENDDO
163        ENDDO 
164     
165     DO j=jj_begin-1,jj_end+1
166       DO i=ii_begin-1,ii_end+1
167         ij=(j-1)*iim+i
168         ps(ij)=100000.0
169         phis(ij)=0.0
170       ENDDO
171     ENDDO
172     
173
174    CALL compute_pression(ps,p,1)     
175    CALL compute_exner(ps,p,pks,pk,1) 
176        theta(:,:)=tetarappel(:,:) 
177    CALL compute_geopotential(phis,pks,pk,theta,phi,1)
178
179        u=0.0 !!wind 0
180!============================================================
181     DO l=1,llm
182      DO j=jj_begin-1,jj_end+1
183        DO i=ii_begin-1,ii_end+1
184          ij=(j-1)*iim+i
185          CALL RANDOM_NUMBER(r); r = 0.0 
186          theta(ij,l)=theta(ij,l)*(1.+0.0005*r)
187       ENDDO
188      ENDDO
189     ENDDO
190    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
191
192  END SUBROUTINE compute_etat0_heldsz
193
194
195        SUBROUTINE held_saurez(f_ps,f_theta_rhodz,f_u) 
196        USE icosa
197            IMPLICIT NONE
198    TYPE(t_field),POINTER :: f_theta_rhodz(:)
199    TYPE(t_field),POINTER :: f_u(:)
200    TYPE(t_field),POINTER :: f_ps(:)
201    REAL(rstd),POINTER :: theta_rhodz(:,:)
202    REAL(rstd),POINTER :: u(:,:)
203    REAL(rstd),POINTER :: ps(:)
204    REAL(rstd),POINTER :: tetarappel(:,:)
205    REAL(rstd),POINTER :: clat(:)
206    INTEGER::ind
207
208    DO ind=1,ndomain
209      CALL swap_dimensions(ind)
210      CALL swap_geometry(ind)
211      theta_rhodz=f_theta_rhodz(ind)
212      u=f_u(ind)
213      ps=f_ps(ind) 
214      tetarappel=f_tetarappel(ind) 
215      clat=f_clat(ind) 
216      CALL compute_heldsz(ps,theta_rhodz,u,clat,tetarappel) 
217    ENDDO
218        END SUBROUTINE held_saurez 
219
220        SUBROUTINE compute_heldsz(ps,theta_rhodz,u,clat,tetarappel) 
221        USE icosa
222        USE theta2theta_rhodz_mod
223           IMPLICIT NONE
224   REAL(rstd),INTENT(IN)::ps(iim*jjm) 
225   REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
226   REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
227   REAL(rstd)::theta(iim*jjm,llm) 
228   REAL(rstd),INTENT(IN)::tetarappel(iim*jjm,llm) 
229   REAL(rstd),INTENT(IN):: clat(iim*jjm) 
230   INTEGER :: i,j,l,ij
231
232      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
233    DO l=1,llm
234      DO j=jj_begin-1,jj_end+1
235        DO i=ii_begin-1,ii_end+1
236          ij=(j-1)*iim+i
237        theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-tetarappel(ij,l))* &
238        (knewt_g+knewt_t(l)*clat(ij)**4 )
239        ENDDO
240      ENDDO
241    ENDDO 
242       CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
243
244        Do l=1,llm
245        u(:,l)=u(:,l)*(1.-dt*kfrict(l))
246        END DO
247
248          END SUBROUTINE compute_heldsz
249
250END MODULE etat0_heldsz_mod
Note: See TracBrowser for help on using the repository browser.