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

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

OpenMP fixes for DCMIP
File size: 9.2 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_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0
22    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
23    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
24    USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
25    ! Ad hoc interfaces
26    USE etat0_academic_mod, ONLY : etat0_academic=>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 ('dcmip4')
68        CALL getin_etat0_dcmip4
69    CASE ('dcmip5')
70        CALL getin_etat0_dcmip5
71    CASE ('williamson91.6')
72       init_mass=.FALSE.
73       CALL getin_etat0_williamson
74    CASE DEFAULT
75       collocated=.FALSE.
76    END SELECT
77
78    !------------------- Ad hoc interfaces --------------------
79    SELECT CASE (TRIM(etat0_type))
80    CASE ('start_file')
81       CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
82    CASE ('academic')
83       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
84    CASE ('held_suarez')
85       PRINT *,"Held & Suarez (1994) test case"
86       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
87    CASE ('venus')
88       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
89       PRINT *, "Venus (Lebonnois et al., 2012) test case"
90   CASE DEFAULT
91      IF(collocated) THEN
92         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
93      ELSE
94         PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
95              '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
96         STOP
97      END IF
98    END SELECT
99
100    IF(init_mass) THEN ! initialize mass distribution using ps
101!       !$OMP BARRIER
102       DO ind=1,ndomain
103          IF (.NOT. assigned_domain(ind)) CYCLE
104          CALL swap_dimensions(ind)
105          CALL swap_geometry(ind)
106          mass=f_mass(ind); ps=f_ps(ind)
107          CALL compute_rhodz(.TRUE., ps, mass)
108       END DO
109    END IF
110
111  END SUBROUTINE etat0
112
113  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
114    USE theta2theta_rhodz_mod
115    IMPLICIT NONE
116    TYPE(t_field),POINTER :: f_ps(:)
117    TYPE(t_field),POINTER :: f_mass(:)
118    TYPE(t_field),POINTER :: f_phis(:)
119    TYPE(t_field),POINTER :: f_theta_rhodz(:)
120    TYPE(t_field),POINTER :: f_u(:)
121    TYPE(t_field),POINTER :: f_q(:)
122 
123    TYPE(t_field),POINTER,SAVE :: f_temp(:)
124    REAL(rstd),POINTER :: ps(:)
125    REAL(rstd),POINTER :: mass(:,:)
126    REAL(rstd),POINTER :: phis(:)
127    REAL(rstd),POINTER :: theta_rhodz(:,:)
128    REAL(rstd),POINTER :: temp(:,:)
129    REAL(rstd),POINTER :: u(:,:)
130    REAL(rstd),POINTER :: q(:,:,:)
131    INTEGER :: ind
132
133    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp')
134
135    DO ind=1,ndomain
136      IF (.NOT. assigned_domain(ind)) CYCLE
137      CALL swap_dimensions(ind)
138      CALL swap_geometry(ind)
139      ps=f_ps(ind)
140      mass=f_mass(ind)
141      phis=f_phis(ind)
142      theta_rhodz=f_theta_rhodz(ind)
143      temp=f_temp(ind)
144      u=f_u(ind)
145      q=f_q(ind)
146
147      IF( TRIM(etat0_type)=='williamson91.6' ) THEN
148       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
149      ELSE
150       CALL compute_etat0_collocated(ps,mass, phis, temp, u, q)
151      ENDIF
152    ENDDO
153   
154    IF( TRIM(etat0_type)/='williamson91.6' ) CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
155   
156    CALL deallocate_field(f_temp)
157   
158  END SUBROUTINE etat0_collocated
159
160  SUBROUTINE compute_etat0_collocated(ps,mass, phis, temp_i, u, q)
161    USE wind_mod
162    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
163    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
164    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0
165    USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
166    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
167    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
168    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
169    USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0
170    IMPLICIT NONE
171    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
172    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
173    REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
174    REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)
175    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
176    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
177
178    REAL(rstd) :: ulon_i(iim*jjm,llm)
179    REAL(rstd) :: ulat_i(iim*jjm,llm)
180
181    REAL(rstd) :: ps_e(3*iim*jjm)
182    REAL(rstd) :: mass_e(3*iim*jjm,llm)
183    REAL(rstd) :: phis_e(3*iim*jjm)
184    REAL(rstd) :: temp_e(3*iim*jjm,llm)
185    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
186    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
187    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
188
189    INTEGER :: l,i,j,ij
190
191    !$OMP BARRIER
192
193    SELECT CASE (TRIM(etat0_type))
194    CASE ('isothermal')
195       CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
196       CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
197    CASE ('temperature_profile')
198       CALL compute_etat0_temperature(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
199       CALL compute_etat0_temperature(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
200    CASE('jablonowsky06')
201       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
202       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
203    CASE('dcmip1')
204       CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
205       CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
206    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
207       CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
208       CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)     
209    CASE('dcmip3')
210       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
211       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
212    CASE('dcmip4')
213       CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
214       CALL compute_dcmip4(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    !$OMP BARRIER
224
225    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
226
227  END SUBROUTINE compute_etat0_collocated
228
229!----------------------------- Resting isothermal state --------------------------------
230
231  SUBROUTINE getin_etat0_isothermal
232    etat0_temp=300
233    CALL getin("etat0_isothermal_temp",etat0_temp)
234  END SUBROUTINE getin_etat0_isothermal
235
236  SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
237    IMPLICIT NONE 
238    INTEGER, INTENT(IN)    :: ngrid
239    REAL(rstd),INTENT(OUT) :: phis(ngrid)
240    REAL(rstd),INTENT(OUT) :: ps(ngrid)
241    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
242    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
243    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
244    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
245    phis(:)=0
246    ps(:)=preff
247    temp(:,:)=etat0_temp
248    ulon(:,:)=0
249    ulat(:,:)=0
250    q(:,:,:)=0
251  END SUBROUTINE compute_etat0_isothermal
252
253END MODULE etat0_mod
Note: See TracBrowser for help on using the repository browser.