source: codes/icosagcm/trunk/src/initial/etat0.f90

Last change on this file was 1017, checked in by ymipsl, 4 years ago

Add possibilty to set a white noise at restart. The amplitude of white noise is fixed by the folowing parameter that can be used in run.def :

  • etat0_ps_white_noise
  • etat0_theta_rhodz_white_noise
  • etat0_u_white_noise
  • etat0_q_white_noise

The value set the relative amplitude of the random perturbation.
ex :

etat0_ps_white_noise=1e-15 -> ps(ij)=ps(ij)*(1+1e-15*ran) where ran inside [0,1[

default perturbation is 0

YM
File size: 18.9 KB
Line 
1MODULE etat0_mod
2  USE icosa
3  USE omp_para
4  IMPLICIT NONE         
5  PRIVATE
6
7    CHARACTER(len=255),SAVE :: etat0_type
8!$OMP THREADPRIVATE(etat0_type)
9
10    REAL(rstd) :: etat0_temp
11
12    PUBLIC :: etat0, init_etat0, etat0_type
13
14! Important notes for OpenMP
15! When etat0 is called, vertical OpenMP parallelism is deactivated.
16! Therefore only the omp_level_master thread must work, i.e. :
17!   !$OMP BARRIER
18!    DO ind=1,ndomain
19!      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
20!      ...
21!    END DO
22!   !$OMP BARRIER
23! There MUST be NO OMP BARRIER inside the DO-LOOP or any routine it calls.
24
25CONTAINS
26 
27  SUBROUTINE init_etat0
28    USE etat0_database_mod, ONLY: init_etat0_database => init_etat0 
29    USE etat0_start_file_mod, ONLY: init_etat0_start_file => init_etat0 
30    USE etat0_heldsz_mod, ONLY: init_etat0_held_suarez => init_etat0 
31   
32    CALL getin("etat0",etat0_type)
33
34    SELECT CASE (TRIM(etat0_type))
35      CASE ('isothermal')
36      CASE ('temperature_profile')
37      CASE ('jablonowsky06')
38      CASE ('dcmip5')
39      CASE ('williamson91.6')
40      CASE ('start_file')
41        CALL init_etat0_start_file
42      CASE ('database')
43        CALL init_etat0_database
44      CASE ('academic')
45      CASE ('held_suarez')
46         CALL init_etat0_held_suarez
47      CASE ('venus')
48      CASE ('dcmip1')
49      CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
50      CASE ('dcmip3')
51      CASE ('dcmip4')
52      CASE ('dcmip2016_baroclinic_wave')
53      CASE ('dcmip2016_cyclone')
54      CASE ('dcmip2016_supercell')
55      CASE ('bubble')
56      CASE DEFAULT
57         PRINT*, 'Bad selector for variable etat0 <',TRIM(etat0_type),'>'// &
58            ' options are  <isothermal>, <temperature_profile>, <jablonowsky06>, <dcmip5>, <williamson91.6>,'& 
59                         //' <start_file>, <database>, <academic>, <held_suarez>, <venus>, <dcmip1>,'         &
60                         //' <dcmip2_mountain,dcmip2_schaer_noshear,dcmip2_schaer_shear>, <dcmip3>, <dcmip4>,'&
61                         //' <dcmip2016_baroclinic_wave>, <dcmip2016_cyclone>, <dcmip2016_supercell>', 'bubble'
62         STOP
63    END SELECT
64
65  END SUBROUTINE init_etat0
66
67  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q)
68    USE disvert_mod
69    ! Generic interface
70    USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0
71    USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0
72    USE etat0_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0
73    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
74    USE etat0_bubble_mod, ONLY : getin_etat0_bubble=>getin_etat0
75    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
76    USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
77    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : getin_etat0_dcmip2016_baroclinic_wave=>getin_etat0
78    USE etat0_dcmip2016_cyclone_mod, ONLY : getin_etat0_dcmip2016_cyclone=>getin_etat0
79    USE etat0_dcmip2016_supercell_mod, ONLY : getin_etat0_dcmip2016_supercell=>getin_etat0
80    ! Ad hoc interfaces
81    USE etat0_academic_mod, ONLY : etat0_academic=>etat0
82    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0
83    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0
84    USE etat0_database_mod, ONLY : etat0_database=>etat0
85    USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0 
86
87    TYPE(t_field),POINTER :: f_ps(:)
88    TYPE(t_field),POINTER :: f_mass(:)
89    TYPE(t_field),POINTER :: f_phis(:)
90    TYPE(t_field),POINTER :: f_theta_rhodz(:)
91    TYPE(t_field),POINTER :: f_u(:)
92    TYPE(t_field),POINTER :: f_geopot(:)
93    TYPE(t_field),POINTER :: f_w(:)
94    TYPE(t_field),POINTER :: f_q(:)
95   
96    REAL(rstd),POINTER :: ps(:), mass(:,:)
97    REAL    :: etat0_ps_white_noise
98    REAL    :: etat0_theta_rhodz_white_noise
99    REAL    :: etat0_u_white_noise
100    REAL    :: etat0_q_white_noise
101    LOGICAL :: autoinit_mass, collocated
102    INTEGER :: ind
103
104    ! most etat0 routines set ps and not mass
105    ! in that case and if caldyn_eta == eta_lag
106    ! the initial distribution of mass is taken to be the same
107    ! as what the mass coordinate would dictate
108    ! however if etat0_XXX defines mass then the flag autoinit_mass must be set to .FALSE.
109    ! otherwise mass will be overwritten
110    autoinit_mass = (caldyn_eta == eta_lag)
111
112    etat0_type='jablonowsky06'
113    CALL getin("etat0",etat0_type)
114   
115    !------------------- Generic interface ---------------------
116    collocated=.TRUE.
117    SELECT CASE (TRIM(etat0_type))
118    CASE ('isothermal')
119       CALL getin_etat0_isothermal
120    CASE ('temperature_profile')
121       CALL getin_etat0_temperature
122    CASE ('jablonowsky06')
123    CASE ('dcmip1')
124        CALL getin_etat0_dcmip1
125    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
126       CALL getin_etat0_dcmip2
127    CASE ('dcmip3')
128    CASE ('dcmip4')
129        CALL getin_etat0_dcmip4
130    CASE ('dcmip5')
131        CALL getin_etat0_dcmip5
132    CASE ('bubble')
133        CALL getin_etat0_bubble
134    CASE ('williamson91.6')
135       autoinit_mass=.FALSE.
136       CALL getin_etat0_williamson
137    CASE ('dcmip2016_baroclinic_wave')
138        CALL getin_etat0_dcmip2016_baroclinic_wave
139    CASE ('dcmip2016_cyclone')
140        CALL getin_etat0_dcmip2016_cyclone
141    CASE ('dcmip2016_supercell')
142        CALL getin_etat0_dcmip2016_supercell
143    CASE DEFAULT
144       collocated=.FALSE.
145       autoinit_mass = .FALSE.
146    END SELECT
147
148    !------------------- Ad hoc interfaces --------------------
149    SELECT CASE (TRIM(etat0_type))
150     CASE ('database')
151        CALL etat0_database(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
152    CASE ('start_file')
153       CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
154    CASE ('academic')
155       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
156    CASE ('held_suarez')
157       PRINT *,"Held & Suarez (1994) test case"
158       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
159    CASE ('venus')
160       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
161       PRINT *, "Venus (Lebonnois et al., 2012) test case"
162   CASE DEFAULT
163      IF(collocated) THEN
164         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
165      ELSE
166         PRINT*, 'Bad selector for variable etat0 <',TRIM(etat0_type),'>'// &
167            ' options are  <isothermal>, <temperature_profile>, <jablonowsky06>, <dcmip5>, <williamson91.6>,'& 
168                         //' <start_file>, <database>, <academic>, <held_suarez>, <venus>, <dcmip1>,'         &
169                         //' <dcmip2_mountain,dcmip2_schaer_noshear,dcmip2_schaer_shear>, <dcmip3>, <dcmip4>,'&
170                         //' <dcmip2016_baroclinic_wave>, <dcmip2016_cyclone>, <dcmip2016_supercell>'
171         STOP
172      END IF
173    END SELECT
174
175    IF(autoinit_mass) THEN
176       DO ind=1,ndomain
177          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
178          CALL swap_dimensions(ind)
179          CALL swap_geometry(ind)
180          mass=f_mass(ind); ps=f_ps(ind)
181          CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
182       END DO
183    END IF
184   
185    etat0_ps_white_noise=0.
186    CALL getin("etat0_ps_white_noise",etat0_ps_white_noise)
187    CALL add_white_noise(f_ps, etat0_ps_white_noise)
188
189    etat0_theta_rhodz_white_noise=0.
190    CALL getin("etat0_theta_rhodz_white_noise",etat0_theta_rhodz_white_noise)
191    CALL add_white_noise(f_theta_rhodz, etat0_theta_rhodz_white_noise)
192 
193    etat0_u_white_noise=0.
194    CALL getin("etat0_u_white_noise",etat0_u_white_noise)
195    CALL add_white_noise(f_u, etat0_u_white_noise)
196
197    etat0_q_white_noise=0.
198    CALL getin("etat0_q_white_noise",etat0_q_white_noise)
199    CALL add_white_noise(f_q, etat0_q_white_noise)
200   
201  END SUBROUTINE etat0
202
203
204  SUBROUTINE add_white_noise(field, factor)
205  USE icosa
206  IMPLICIT NONE
207    TYPE(t_field),POINTER :: field(:)   ! INOUT
208    REAL,INTENT(IN) :: factor   
209   
210    INTEGER,ALLOCATABLE :: seed(:)
211    REAL,ALLOCATABLE :: random2d(:)
212    REAL,ALLOCATABLE :: random3d(:,:)
213    REAL,ALLOCATABLE :: random4d(:,:,:)
214    REAL,POINTER :: field2d(:)
215    REAL,POINTER :: field3d(:,:)
216    REAL,POINTER :: field4d(:,:,:)
217    INTEGER :: ind
218    INTEGER :: m
219   
220    CALL RANDOM_SEED(SIZE=m)
221    ALLOCATE(seed(m))
222   
223    DO ind=1,ndomain
224      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
225      seed=domloc_glo_ind(ind) ! to be reproducible
226      CALL RANDOM_SEED(PUT=seed)
227      IF (field(ind)%ndim==2) THEN
228        field2d=field(ind)
229        ALLOCATE(random2d(size(field2d,1)))
230        CALL RANDOM_NUMBER(random2d)
231        field2d=field2d*(1.+random2d*factor)
232        DEALLOCATE(random2d)
233      ELSE IF (field(ind)%ndim==3) THEN
234        field3d=field(ind)
235        ALLOCATE(random3d(size(field3d,1),size(field3d,2)))
236        CALL RANDOM_NUMBER(random3d)
237        field3d=field3d*(1.+random3d*factor)
238        DEALLOCATE(random3d)
239      ELSE IF (field(ind)%ndim==4) THEN
240        field4d=field(ind)
241        ALLOCATE(random4d(size(field4d,1),size(field4d,2),size(field4d,3)))
242        CALL RANDOM_NUMBER(random4d)
243        field4d=field4d*(1.+random4d*factor)
244        DEALLOCATE(random4d)
245      ENDIF
246    ENDDO
247   
248  END SUBROUTINE
249   
250   
251  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
252    USE theta2theta_rhodz_mod
253    TYPE(t_field),POINTER :: f_ps(:)
254    TYPE(t_field),POINTER :: f_mass(:)
255    TYPE(t_field),POINTER :: f_phis(:)
256    TYPE(t_field),POINTER :: f_theta_rhodz(:)
257    TYPE(t_field),POINTER :: f_u(:)
258    TYPE(t_field),POINTER :: f_geopot(:)
259    TYPE(t_field),POINTER :: f_W(:)
260    TYPE(t_field),POINTER :: f_q(:)
261 
262    TYPE(t_field),POINTER,SAVE :: f_temp(:)
263    REAL(rstd),POINTER :: ps(:)
264    REAL(rstd),POINTER :: mass(:,:)
265    REAL(rstd),POINTER :: phis(:)
266    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
267    REAL(rstd),POINTER :: temp(:,:)
268    REAL(rstd),POINTER :: u(:,:)
269    REAL(rstd),POINTER :: geopot(:,:)
270    REAL(rstd),POINTER :: W(:,:)
271    REAL(rstd),POINTER :: q(:,:,:)
272    INTEGER :: ind
273
274    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp')
275
276    DO ind=1,ndomain
277      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
278!      IF (.NOT. assigned_domain(ind)) CYCLE
279      CALL swap_dimensions(ind)
280      CALL swap_geometry(ind)
281      ps=f_ps(ind)
282      mass=f_mass(ind)
283      phis=f_phis(ind)
284      theta_rhodz=f_theta_rhodz(ind)
285      temp=f_temp(ind)
286      u=f_u(ind)
287      geopot=f_geopot(ind)
288      w=f_w(ind)
289      q=f_q(ind)
290
291      IF( TRIM(etat0_type)=='williamson91.6' ) THEN
292         CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q)
293      ELSE
294         CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
295      ENDIF
296
297      IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1)
298   
299    ENDDO
300   
301    CALL deallocate_field(f_temp)
302   
303  END SUBROUTINE etat0_collocated
304
305  SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset)
306    USE icosa
307    USE pression_mod
308    USE exner_mod
309    USE omp_para
310    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
311    REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
312    REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot)
313    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
314    INTEGER,INTENT(IN) :: offset
315
316    REAL(rstd) :: p(iim*jjm,llm+1)
317    REAL(rstd) :: cppd,Rd, mass, p_ij, chi,nu, entropy, theta
318    INTEGER :: i,j,ij,l
319
320    cppd=cpp
321    Rd=kappa*cppd
322
323    CALL compute_pression(ps,p,offset)
324    ! flush p
325    DO    l    = ll_begin, ll_end
326       DO j=jj_begin-offset,jj_end+offset
327          DO i=ii_begin-offset,ii_end+offset
328             ij=(j-1)*iim+i
329             mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass
330             p_ij = .5*(p(ij,l)+p(ij,l+1))  ! pressure at full level
331             SELECT CASE(caldyn_thermo)
332             CASE(thermo_theta)
333                theta = temp(ij,l)*(p_ij/preff)**(-kappa) 
334                theta_rhodz(ij,l) = mass * theta
335             CASE(thermo_entropy)
336                nu = log(p_ij/preff)
337                chi = log(temp(ij,l)/Treff)
338                entropy = cppd*chi-Rd*nu
339                theta_rhodz(ij,l) = mass * entropy
340!             CASE(thermo_moist)
341!                q_ij=q(ij,l,1)
342!                r_ij=1.-q_ij
343!                mass=mass*(1-q_ij) ! dry mass
344!                nu = log(p_ij/preff)
345!                chi = log(temp(ij,l)/Treff)
346!                entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu)
347!                theta_rhodz(ij,l) = mass * entropy               
348                CASE DEFAULT
349                   STOP
350             END SELECT
351          ENDDO
352       ENDDO
353    ENDDO
354  END SUBROUTINE compute_temperature2entropy
355
356  SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)
357    USE wind_mod
358    USE disvert_mod
359    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
360    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
361    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0
362    USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
363    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
364    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
365    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0 
366    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
367    USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0
368    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0
369    USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0
370    USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0
371    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
372    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
373    REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
374    REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)
375    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
376    REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1)
377    REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1)
378    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
379
380    REAL(rstd) :: ulon_i(iim*jjm,llm)
381    REAL(rstd) :: ulat_i(iim*jjm,llm)
382
383    REAL(rstd) :: ps_e(3*iim*jjm)
384    REAL(rstd) :: mass_e(3*iim*jjm,llm)
385    REAL(rstd) :: phis_e(3*iim*jjm)
386    REAL(rstd) :: temp_e(3*iim*jjm,llm)
387    REAL(rstd) :: geopot_e(3*iim*jjm,llm+1)
388    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
389    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
390    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
391
392    INTEGER :: l,ij
393    REAL :: p_ik, v_ik, mass_ik
394    LOGICAL :: autoinit_mass, autoinit_NH
395
396    ! For NH geopotential and vertical momentum must be initialized.
397    ! Unless autoinit_NH is set to .FALSE. , they will be initialized
398    ! to hydrostatic geopotential and zero
399    autoinit_mass = .TRUE.
400    autoinit_NH = .NOT. hydrostatic
401    w(:,:) = 0
402
403    SELECT CASE (TRIM(etat0_type))
404    CASE ('isothermal')
405       CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
406       CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
407    CASE ('temperature_profile')
408       CALL compute_etat0_temperature(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
409       CALL compute_etat0_temperature(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
410    CASE('jablonowsky06')
411       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
412       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
413    CASE('dcmip1')
414       CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
415       CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
416    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
417       CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
418       CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)     
419    CASE('dcmip3')
420       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
421       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
422       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot
423    CASE('dcmip4')
424       CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
425       CALL compute_dcmip4(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
426    CASE('dcmip5')
427       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
428       CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
429    CASE('bubble')
430       CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
431       CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
432!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot
433    CASE('williamson91.6')
434       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1))
435       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))
436       autoinit_mass = .FALSE. ! do not overwrite mass
437    CASE('dcmip2016_baroclinic_wave')
438       CALL compute_dcmip2016_baroclinic_wave(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
439       CALL compute_dcmip2016_baroclinic_wave(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
440    CASE('dcmip2016_cyclone')
441       CALL compute_dcmip2016_cyclone(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
442       CALL compute_dcmip2016_cyclone(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
443    CASE('dcmip2016_supercell')
444       CALL compute_dcmip2016_supercell(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
445       CALL compute_dcmip2016_supercell(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
446    END SELECT
447
448    IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
449    IF(autoinit_NH) THEN
450       geopot(:,1) = phis(:) ! surface geopotential
451       DO l = 1, llm
452          DO ij=1,iim*jjm
453             ! hybrid pressure coordinate
454             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
455             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
456             ! v=R.T/p, R=kappa*cpp
457             v_ik = kappa*cpp*temp_i(ij,l)/p_ik
458             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g
459          END DO
460       END DO
461    END IF
462
463    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
464
465  END SUBROUTINE compute_etat0_collocated
466
467!----------------------------- Resting isothermal state --------------------------------
468
469  SUBROUTINE getin_etat0_isothermal
470    etat0_temp=300
471    CALL getin("etat0_isothermal_temp",etat0_temp)
472  END SUBROUTINE getin_etat0_isothermal
473
474  SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
475    INTEGER, INTENT(IN)    :: ngrid
476    REAL(rstd),INTENT(OUT) :: phis(ngrid)
477    REAL(rstd),INTENT(OUT) :: ps(ngrid)
478    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
479    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
480    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
481    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
482    phis(:)=0
483    ps(:)=preff
484    temp(:,:)=etat0_temp
485    ulon(:,:)=0
486    ulat(:,:)=0
487    q(:,:,:)=0
488  END SUBROUTINE compute_etat0_isothermal
489
490END MODULE etat0_mod
Note: See TracBrowser for help on using the repository browser.