source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/etat0.f90 @ 260

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

Implement restartability for dynamico

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.