source: codes/icosagcm/trunk/src/etat0.f90 @ 215

Last change on this file since 215 was 204, checked in by dubos, 10 years ago

Upgraded etat0_williamson to new interface (tested)
File size: 8.4 KB
Line 
1MODULE etat0_mod
2  USE icosa
3  PRIVATE
4
5    CHARACTER(len=255),SAVE :: etat0_type
6!$OMP THREADPRIVATE(etat0_type)
7
8    REAL(rstd) :: etat0_temp
9
10    PUBLIC :: etat0, etat0_type
11
12CONTAINS
13 
14  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
15    USE mpipara, ONLY : is_mpi_root
16    USE disvert_mod
17    ! New interface
18    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
19    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
20    ! Old interface
21    USE etat0_academic_mod, ONLY : etat0_academic=>etat0 
22    USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0
23    USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0
24    USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0 
25    USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 
26    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 
27    USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 
28    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0 
29
30    IMPLICIT NONE
31    TYPE(t_field),POINTER :: f_ps(:)
32    TYPE(t_field),POINTER :: f_mass(:)
33    TYPE(t_field),POINTER :: f_phis(:)
34    TYPE(t_field),POINTER :: f_theta_rhodz(:)
35    TYPE(t_field),POINTER :: f_u(:)
36    TYPE(t_field),POINTER :: f_q(:)
37   
38    REAL(rstd),POINTER :: ps(:), mass(:,:)
39    LOGICAL :: init_mass
40    INTEGER :: ind,i,j,ij,l
41
42    ! most etat0 routines set ps and not mass
43    ! in that case and if caldyn_eta == eta_lag
44    ! the initial distribution of mass is taken to be the same
45    ! as what the mass coordinate would dictate
46    ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE.
47    ! otherwise mass will be overwritten
48    init_mass = (caldyn_eta == eta_lag)
49
50    etat0_type='jablonowsky06'
51    CALL getin("etat0",etat0_type)
52   
53    SELECT CASE (TRIM(etat0_type))
54       !------------------- New interface ---------------------
55    CASE ('isothermal')
56       CALL getin_etat0_isothermal
57       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
58    CASE ('jablonowsky06')
59       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
60     CASE ('dcmip5')
61        CALL getin_etat0_dcmip5
62        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
63    CASE ('williamson91.6')
64       init_mass=.FALSE.
65       CALL getin_etat0_williamson
66       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
67       !------------------- Old interface --------------------
68    CASE ('academic')
69       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
70    CASE ('held_suarez')
71       PRINT *,"Held & Suarez (1994) test case"
72       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
73    CASE ('dcmip1')
74       CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
75    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
76       CALL etat0_dcmip2(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
77    CASE ('dcmip3')
78       CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
79     CASE ('dcmip4')
80        IF(nqtot<2) THEN
81           IF (is_mpi_root)  THEN
82              PRINT *, "nqtot must be at least 2 for test case DCMIP4"
83           END IF
84           STOP
85        END IF
86       CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
87     CASE ('readnf_start') 
88          print*,"readnf_start used"   
89       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
90        CASE ('readnf_hz') 
91          print*,"readnf_hz used"
92       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
93   CASE DEFAULT
94       PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
95            '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
96       STOP
97    END SELECT
98
99    IF(init_mass) THEN ! initialize mass distribution using ps
100!       !$OMP BARRIER
101       DO ind=1,ndomain
102          IF (.NOT. assigned_domain(ind)) CYCLE
103          CALL swap_dimensions(ind)
104          CALL swap_geometry(ind)
105          mass=f_mass(ind); ps=f_ps(ind)
106          CALL compute_rhodz(.TRUE., ps, mass)
107       END DO
108    END IF
109
110  END SUBROUTINE etat0
111
112  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
113    IMPLICIT NONE
114    TYPE(t_field),POINTER :: f_ps(:)
115    TYPE(t_field),POINTER :: f_mass(:)
116    TYPE(t_field),POINTER :: f_phis(:)
117    TYPE(t_field),POINTER :: f_theta_rhodz(:)
118    TYPE(t_field),POINTER :: f_u(:)
119    TYPE(t_field),POINTER :: f_q(:)
120 
121    REAL(rstd),POINTER :: ps(:)
122    REAL(rstd),POINTER :: mass(:,:)
123    REAL(rstd),POINTER :: phis(:)
124    REAL(rstd),POINTER :: theta_rhodz(:,:)
125    REAL(rstd),POINTER :: u(:,:)
126    REAL(rstd),POINTER :: q(:,:,:)
127    INTEGER :: ind
128
129    DO ind=1,ndomain
130      IF (.NOT. assigned_domain(ind)) CYCLE
131      CALL swap_dimensions(ind)
132      CALL swap_geometry(ind)
133      ps=f_ps(ind)
134      mass=f_mass(ind)
135      phis=f_phis(ind)
136      theta_rhodz=f_theta_rhodz(ind)
137      u=f_u(ind)
138      q=f_q(ind)
139      CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
140    ENDDO
141  END SUBROUTINE etat0_collocated
142
143  SUBROUTINE compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
144    USE theta2theta_rhodz_mod
145    USE wind_mod
146    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
147    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
148    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
149    IMPLICIT NONE
150    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
151    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
152    REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
153    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
154    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
155    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
156
157    REAL(rstd) :: lon_i(iim*jjm)
158    REAL(rstd) :: lat_i(iim*jjm)
159    REAL(rstd) :: temp_i(iim*jjm,llm)
160    REAL(rstd) :: ulon_i(iim*jjm,llm)
161    REAL(rstd) :: ulat_i(iim*jjm,llm)
162
163    REAL(rstd) :: lon_e(3*iim*jjm)
164    REAL(rstd) :: lat_e(3*iim*jjm)
165    REAL(rstd) :: ps_e(3*iim*jjm)
166    REAL(rstd) :: mass_e(3*iim*jjm,llm)
167    REAL(rstd) :: phis_e(3*iim*jjm)
168    REAL(rstd) :: temp_e(3*iim*jjm,llm)
169    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
170    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
171    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
172
173    INTEGER :: l,i,j,ij
174
175    DO l=1,llm
176       DO j=jj_begin-1,jj_end+1
177          DO i=ii_begin-1,ii_end+1
178             ij=(j-1)*iim+i
179             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon_i(ij),lat_i(ij))
180             CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon_e(ij+u_right),lat_e(ij+u_right))
181             CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon_e(ij+u_lup),lat_e(ij+u_lup))
182             CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon_e(ij+u_ldown),lat_e(ij+u_ldown))
183          END DO
184       END DO
185    END DO
186
187    SELECT CASE (TRIM(etat0_type))
188    CASE ('isothermal')
189       CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
190       CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
191    CASE('jablonowsky06')
192       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
193       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
194    CASE('dcmip5')
195       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
196       CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
197    CASE('williamson91.6')
198       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), theta_rhodz(:,1), ulon_i(:,1), ulat_i(:,1))
199       CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1))
200    END SELECT
201
202    SELECT CASE (TRIM(etat0_type))
203    CASE('williamson91.6')
204       ! Do nothing
205    CASE DEFAULT
206       CALL compute_temperature2theta_rhodz(ps,temp_i,theta_rhodz,0)   
207    END SELECT
208 
209    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
210
211  END SUBROUTINE compute_etat0_collocated
212
213!----------------------------- Resting isothermal state --------------------------------
214
215  SUBROUTINE getin_etat0_isothermal
216    etat0_temp=300
217    CALL getin("etat0_isothermal_temp",etat0_temp)
218  END SUBROUTINE getin_etat0_isothermal
219
220  SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
221    IMPLICIT NONE 
222    INTEGER, INTENT(IN)    :: ngrid
223    REAL(rstd),INTENT(OUT) :: phis(ngrid)
224    REAL(rstd),INTENT(OUT) :: ps(ngrid)
225    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
226    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
227    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
228    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
229    phis(:)=0
230    ps(:)=preff
231    temp(:,:)=etat0_temp
232    ulon(:,:)=0
233    ulat(:,:)=0
234    q(:,:,:)=0
235  END SUBROUTINE compute_etat0_isothermal
236
237END MODULE etat0_mod
Note: See TracBrowser for help on using the repository browser.