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

Last change on this file since 499 was 499, checked in by dubos, 8 years ago

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