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

Last change on this file since 344 was 344, checked in by dubos, 9 years ago

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