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

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

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