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

Last change on this file since 385 was 382, checked in by ymipsl, 8 years ago

Add etat0 for DCMIP2016 baroclinic wave testcase.

YM
File size: 11.8 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_geopot,f_w, 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_bubble_mod, ONLY : getin_etat0_bubble=>getin_etat0
24    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
25    USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
26    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : getin_etat0_dcmip2016_baroclinic_wave=>getin_etat0
27    ! Ad hoc interfaces
28    USE etat0_academic_mod, ONLY : etat0_academic=>etat0 
29    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 
30    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0 
31    USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0 
32
33    IMPLICIT NONE
34    TYPE(t_field),POINTER :: f_ps(:)
35    TYPE(t_field),POINTER :: f_mass(:)
36    TYPE(t_field),POINTER :: f_phis(:)
37    TYPE(t_field),POINTER :: f_theta_rhodz(:)
38    TYPE(t_field),POINTER :: f_u(:)
39    TYPE(t_field),POINTER :: f_geopot(:)
40    TYPE(t_field),POINTER :: f_w(:)
41    TYPE(t_field),POINTER :: f_q(:)
42   
43    REAL(rstd),POINTER :: ps(:), mass(:,:)
44    LOGICAL :: autoinit_mass, autoinit_geopot, collocated
45    INTEGER :: ind,i,j,ij,l
46
47    ! most etat0 routines set ps and not mass
48    ! in that case and if caldyn_eta == eta_lag
49    ! the initial distribution of mass is taken to be the same
50    ! as what the mass coordinate would dictate
51    ! however if etat0_XXX defines mass then the flag autoinit_mass must be set to .FALSE.
52    ! otherwise mass will be overwritten
53    autoinit_mass = (caldyn_eta == eta_lag)
54
55    etat0_type='jablonowsky06'
56    CALL getin("etat0",etat0_type)
57   
58    !------------------- Generic interface ---------------------
59    collocated=.TRUE.
60    SELECT CASE (TRIM(etat0_type))
61    CASE ('isothermal')
62       CALL getin_etat0_isothermal
63    CASE ('temperature_profile')
64       CALL getin_etat0_temperature
65    CASE ('jablonowsky06')
66    CASE ('dcmip1')
67        CALL getin_etat0_dcmip1
68    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
69       CALL getin_etat0_dcmip2
70    CASE ('dcmip3')
71    CASE ('dcmip4')
72        CALL getin_etat0_dcmip4
73    CASE ('dcmip5')
74        CALL getin_etat0_dcmip5
75    CASE ('bubble')
76        CALL getin_etat0_bubble
77    CASE ('williamson91.6')
78       autoinit_mass=.FALSE.
79       CALL getin_etat0_williamson
80    CASE ('dcmip2016_baroclinic_wave')
81        CALL getin_etat0_dcmip2016_baroclinic_wave
82    CASE DEFAULT
83       collocated=.FALSE.
84       autoinit_mass = .FALSE.
85    END SELECT
86
87    !------------------- Ad hoc interfaces --------------------
88    SELECT CASE (TRIM(etat0_type))
89    CASE ('start_file')
90       CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
91    CASE ('academic')
92       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
93    CASE ('held_suarez')
94       PRINT *,"Held & Suarez (1994) test case"
95       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
96    CASE ('venus')
97       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
98       PRINT *, "Venus (Lebonnois et al., 2012) test case"
99   CASE DEFAULT
100      IF(collocated) THEN
101         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
102      ELSE
103         PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
104              '> options are <isothermal>, <temperature_profile>, <held_suarez>, &
105              &<bubble>, <venus>, <jablonowsky06>, <academic>, <dcmip[1-4]> '
106         STOP
107      END IF
108    END SELECT
109
110!       !$OMP BARRIER
111    IF(autoinit_mass) THEN
112       DO ind=1,ndomain
113          IF (.NOT. assigned_domain(ind)) CYCLE
114          CALL swap_dimensions(ind)
115          CALL swap_geometry(ind)
116          mass=f_mass(ind); ps=f_ps(ind)
117          CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
118       END DO
119    END IF
120 
121  END SUBROUTINE etat0
122
123  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
124    USE theta2theta_rhodz_mod
125    IMPLICIT NONE
126    TYPE(t_field),POINTER :: f_ps(:)
127    TYPE(t_field),POINTER :: f_mass(:)
128    TYPE(t_field),POINTER :: f_phis(:)
129    TYPE(t_field),POINTER :: f_theta_rhodz(:)
130    TYPE(t_field),POINTER :: f_u(:)
131    TYPE(t_field),POINTER :: f_geopot(:)
132    TYPE(t_field),POINTER :: f_W(:)
133    TYPE(t_field),POINTER :: f_q(:)
134 
135    TYPE(t_field),POINTER,SAVE :: f_temp(:)
136    REAL(rstd),POINTER :: ps(:)
137    REAL(rstd),POINTER :: mass(:,:)
138    REAL(rstd),POINTER :: phis(:)
139    REAL(rstd),POINTER :: theta_rhodz(:,:)
140    REAL(rstd),POINTER :: temp(:,:)
141    REAL(rstd),POINTER :: u(:,:)
142    REAL(rstd),POINTER :: geopot(:,:)
143    REAL(rstd),POINTER :: W(:,:)
144    REAL(rstd),POINTER :: q(:,:,:)
145    INTEGER :: ind
146
147    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp')
148
149    DO ind=1,ndomain
150      IF (.NOT. assigned_domain(ind)) CYCLE
151      CALL swap_dimensions(ind)
152      CALL swap_geometry(ind)
153      ps=f_ps(ind)
154      mass=f_mass(ind)
155      phis=f_phis(ind)
156      theta_rhodz=f_theta_rhodz(ind)
157      temp=f_temp(ind)
158      u=f_u(ind)
159      geopot=f_geopot(ind)
160      w=f_w(ind)
161      q=f_q(ind)
162
163      IF( TRIM(etat0_type)=='williamson91.6' ) THEN
164       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, geopot, W, q)
165      ELSE
166       CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
167      ENDIF
168    ENDDO
169   
170    IF( TRIM(etat0_type)/='williamson91.6' ) CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
171   
172    CALL deallocate_field(f_temp)
173   
174  END SUBROUTINE etat0_collocated
175
176  SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)
177    USE wind_mod
178    USE disvert_mod
179    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
180    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
181    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0
182    USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
183    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
184    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
185    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0 
186    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
187    USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0
188    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0
189    IMPLICIT NONE
190    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
191    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
192    REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
193    REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)
194    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
195    REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1)
196    REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1)
197    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
198
199    REAL(rstd) :: ulon_i(iim*jjm,llm)
200    REAL(rstd) :: ulat_i(iim*jjm,llm)
201
202    REAL(rstd) :: ps_e(3*iim*jjm)
203    REAL(rstd) :: mass_e(3*iim*jjm,llm)
204    REAL(rstd) :: phis_e(3*iim*jjm)
205    REAL(rstd) :: temp_e(3*iim*jjm,llm)
206    REAL(rstd) :: geopot_e(3*iim*jjm,llm+1)
207    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
208    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
209    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
210
211    INTEGER :: l,i,j,ij
212    REAL :: p_ik, v_ik, mass_ik
213    LOGICAL :: autoinit_mass, autoinit_NH
214
215    ! For NH geopotential and vertical momentum must be initialized.
216    ! Unless autoinit_NH is set to .FALSE. , they will be initialized
217    ! to hydrostatic geopotential and zero
218    autoinit_mass = .TRUE.
219    autoinit_NH = .NOT. hydrostatic
220    w(:,:) = 0
221
222    !$OMP BARRIER
223
224    SELECT CASE (TRIM(etat0_type))
225    CASE ('isothermal')
226       CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
227       CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
228    CASE ('temperature_profile')
229       CALL compute_etat0_temperature(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
230       CALL compute_etat0_temperature(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
231    CASE('jablonowsky06')
232       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
233       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
234    CASE('dcmip1')
235       CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
236       CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
237    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
238       CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
239       CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)     
240    CASE('dcmip3')
241       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
242       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
243       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot
244    CASE('dcmip4')
245       CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
246       CALL compute_dcmip4(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
247    CASE('dcmip5')
248       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
249       CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
250    CASE('bubble')
251       CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
252       CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
253!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot
254    CASE('williamson91.6')
255       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1))
256       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))
257       autoinit_mass = .FALSE. ! do not overwrite mass
258    CASE('dcmip2016_baroclinic_wave')
259       CALL compute_dcmip2016_baroclinic_wave(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
260       CALL compute_dcmip2016_baroclinic_wave(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
261    END SELECT
262
263    IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
264    IF(autoinit_NH) THEN
265       geopot(:,1) = phis(:) ! surface geopotential
266       DO l = 1, llm
267          DO ij=1,iim*jjm
268             ! hybrid pressure coordinate
269             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
270             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
271             ! v=R.T/p, R=kappa*cpp
272             v_ik = kappa*cpp*temp_i(ij,l)/p_ik
273             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g
274          END DO
275       END DO
276    END IF
277
278    !$OMP BARRIER
279
280    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
281
282  END SUBROUTINE compute_etat0_collocated
283
284!----------------------------- Resting isothermal state --------------------------------
285
286  SUBROUTINE getin_etat0_isothermal
287    etat0_temp=300
288    CALL getin("etat0_isothermal_temp",etat0_temp)
289  END SUBROUTINE getin_etat0_isothermal
290
291  SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
292    IMPLICIT NONE 
293    INTEGER, INTENT(IN)    :: ngrid
294    REAL(rstd),INTENT(OUT) :: phis(ngrid)
295    REAL(rstd),INTENT(OUT) :: ps(ngrid)
296    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
297    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
298    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
299    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
300    phis(:)=0
301    ps(:)=preff
302    temp(:,:)=etat0_temp
303    ulon(:,:)=0
304    ulat(:,:)=0
305    q(:,:,:)=0
306  END SUBROUTINE compute_etat0_isothermal
307
308END MODULE etat0_mod
Note: See TracBrowser for help on using the repository browser.