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

Last change on this file since 283 was 266, checked in by ymipsl, 10 years ago

Synchronize trunk and Saturn branch.
Merge modification from Saturn branch to trunk

YM
File size: 8.5 KB
Line 
1MODULE etat0_mod
2  USE icosa
3  PRIVATE
4
5    CHARACTER(len=255),SAVE :: etat0_type
6!$OMP THREADPRIVATE(etat0_type)
7
8    REAL(rstd) :: etat0_temp
9
10    PUBLIC :: etat0, etat0_type
11
12CONTAINS
13 
14  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
15    USE mpipara, ONLY : is_mpi_root
16    USE disvert_mod
17    ! New interface
18    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
19    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
20    ! Old interface
21    USE etat0_academic_mod, ONLY : etat0_academic=>etat0 
22    USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0
23    USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0
24    USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0 
25    USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 
26    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 
27    USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 
28    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>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
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    SELECT CASE (TRIM(etat0_type))
55       !------------------- New interface ---------------------
56    CASE ('isothermal')
57       CALL getin_etat0_isothermal
58       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
59    CASE ('jablonowsky06')
60       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
61     CASE ('dcmip5')
62        CALL getin_etat0_dcmip5
63        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
64    CASE ('williamson91.6')
65       init_mass=.FALSE.
66       CALL getin_etat0_williamson
67       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
68       !------------------- Old interface --------------------
69    CASE ('start_file')
70       CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
71    CASE ('academic')
72       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
73    CASE ('held_suarez')
74       PRINT *,"Held & Suarez (1994) test case"
75       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
76    CASE ('dcmip1')
77       CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
78    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
79       CALL etat0_dcmip2(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
80    CASE ('dcmip3')
81       CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
82     CASE ('dcmip4')
83        IF(nqtot<2) THEN
84           IF (is_mpi_root)  THEN
85              PRINT *, "nqtot must be at least 2 for test case DCMIP4"
86           END IF
87           STOP
88        END IF
89       CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
90     CASE ('readnf_start') 
91          print*,"readnf_start used"   
92       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
93        CASE ('readnf_hz') 
94          print*,"readnf_hz used"
95       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
96   CASE DEFAULT
97       PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
98            '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
99       STOP
100    END SELECT
101
102    IF(init_mass) THEN ! initialize mass distribution using ps
103!       !$OMP BARRIER
104       DO ind=1,ndomain
105          IF (.NOT. assigned_domain(ind)) CYCLE
106          CALL swap_dimensions(ind)
107          CALL swap_geometry(ind)
108          mass=f_mass(ind); ps=f_ps(ind)
109          CALL compute_rhodz(.TRUE., ps, mass)
110       END DO
111    END IF
112
113  END SUBROUTINE etat0
114
115  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
116    IMPLICIT NONE
117    TYPE(t_field),POINTER :: f_ps(:)
118    TYPE(t_field),POINTER :: f_mass(:)
119    TYPE(t_field),POINTER :: f_phis(:)
120    TYPE(t_field),POINTER :: f_theta_rhodz(:)
121    TYPE(t_field),POINTER :: f_u(:)
122    TYPE(t_field),POINTER :: f_q(:)
123 
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 :: u(:,:)
129    REAL(rstd),POINTER :: q(:,:,:)
130    INTEGER :: ind
131
132    DO ind=1,ndomain
133      IF (.NOT. assigned_domain(ind)) CYCLE
134      CALL swap_dimensions(ind)
135      CALL swap_geometry(ind)
136      ps=f_ps(ind)
137      mass=f_mass(ind)
138      phis=f_phis(ind)
139      theta_rhodz=f_theta_rhodz(ind)
140      u=f_u(ind)
141      q=f_q(ind)
142      CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
143    ENDDO
144  END SUBROUTINE etat0_collocated
145
146  SUBROUTINE compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
147    USE theta2theta_rhodz_mod
148    USE wind_mod
149    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
150    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
151    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
152    IMPLICIT NONE
153    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
154    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
155    REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
156    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
157    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
158    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
159
160    REAL(rstd) :: lon_i(iim*jjm)
161    REAL(rstd) :: lat_i(iim*jjm)
162    REAL(rstd) :: temp_i(iim*jjm,llm)
163    REAL(rstd) :: ulon_i(iim*jjm,llm)
164    REAL(rstd) :: ulat_i(iim*jjm,llm)
165
166    REAL(rstd) :: lon_e(3*iim*jjm)
167    REAL(rstd) :: lat_e(3*iim*jjm)
168    REAL(rstd) :: ps_e(3*iim*jjm)
169    REAL(rstd) :: mass_e(3*iim*jjm,llm)
170    REAL(rstd) :: phis_e(3*iim*jjm)
171    REAL(rstd) :: temp_e(3*iim*jjm,llm)
172    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
173    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
174    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
175
176    INTEGER :: l,i,j,ij
177
178    DO l=1,llm
179       DO j=jj_begin-1,jj_end+1
180          DO i=ii_begin-1,ii_end+1
181             ij=(j-1)*iim+i
182             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon_i(ij),lat_i(ij))
183             CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon_e(ij+u_right),lat_e(ij+u_right))
184             CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon_e(ij+u_lup),lat_e(ij+u_lup))
185             CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon_e(ij+u_ldown),lat_e(ij+u_ldown))
186          END DO
187       END DO
188    END DO
189
190    SELECT CASE (TRIM(etat0_type))
191    CASE ('isothermal')
192       CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
193       CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
194    CASE('jablonowsky06')
195       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
196       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
197    CASE('dcmip5')
198       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
199       CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
200    CASE('williamson91.6')
201       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), theta_rhodz(:,1), ulon_i(:,1), ulat_i(:,1))
202       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))
203    END SELECT
204
205    SELECT CASE (TRIM(etat0_type))
206    CASE('williamson91.6')
207       ! Do nothing
208    CASE DEFAULT
209       CALL compute_temperature2theta_rhodz(ps,temp_i,theta_rhodz,0)   
210    END SELECT
211 
212    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
213
214  END SUBROUTINE compute_etat0_collocated
215
216!----------------------------- Resting isothermal state --------------------------------
217
218  SUBROUTINE getin_etat0_isothermal
219    etat0_temp=300
220    CALL getin("etat0_isothermal_temp",etat0_temp)
221  END SUBROUTINE getin_etat0_isothermal
222
223  SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
224    IMPLICIT NONE 
225    INTEGER, INTENT(IN)    :: ngrid
226    REAL(rstd),INTENT(OUT) :: phis(ngrid)
227    REAL(rstd),INTENT(OUT) :: ps(ngrid)
228    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
229    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
230    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
231    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
232    phis(:)=0
233    ps(:)=preff
234    temp(:,:)=etat0_temp
235    ulon(:,:)=0
236    ulat(:,:)=0
237    q(:,:,:)=0
238  END SUBROUTINE compute_etat0_isothermal
239
240END MODULE etat0_mod
Note: See TracBrowser for help on using the repository browser.