source: codes/icosagcm/devel/src/initial/etat0_collocated.f90 @ 884

Last change on this file since 884 was 884, checked in by dubos, 5 years ago

devel/unstructured : fixed etat0 + renamed compute => compute_caldyn
File size: 10.6 KB
Line 
1MODULE etat0_collocated_mod
2  USE icosa
3  USE grid_param
4  USE omp_para, ONLY : is_omp_level_master
5  USE caldyn_vars_mod, ONLY : hydrostatic
6  IMPLICIT NONE         
7  PRIVATE
8
9  CHARACTER(len=255),SAVE :: etat0_type
10!$OMP THREADPRIVATE(etat0_type)
11
12    PUBLIC :: etat0_type, etat0_collocated
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 etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
28    USE theta2theta_rhodz_mod
29    TYPE(t_field),POINTER :: f_ps(:)
30    TYPE(t_field),POINTER :: f_mass(:)
31    TYPE(t_field),POINTER :: f_phis(:)
32    TYPE(t_field),POINTER :: f_theta_rhodz(:)
33    TYPE(t_field),POINTER :: f_u(:)
34    TYPE(t_field),POINTER :: f_geopot(:)
35    TYPE(t_field),POINTER :: f_W(:)
36    TYPE(t_field),POINTER :: f_q(:)
37 
38    TYPE(t_field),POINTER,SAVE :: f_temp(:) ! SAVE => all threads see the same f_temp
39    REAL(rstd),POINTER :: ps(:)
40    REAL(rstd),POINTER :: mass(:,:)
41    REAL(rstd),POINTER :: phis(:)
42    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
43    REAL(rstd),POINTER :: temp(:,:)
44    REAL(rstd),POINTER :: u(:,:)
45    REAL(rstd),POINTER :: geopot(:,:)
46    REAL(rstd),POINTER :: W(:,:)
47    REAL(rstd),POINTER :: q(:,:,:)
48    INTEGER :: ind
49
50    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp')
51
52    DO ind=1,ndomain
53      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
54      CALL swap_dimensions(ind)
55      CALL swap_geometry(ind)
56      ps=f_ps(ind)
57      mass=f_mass(ind)
58      phis=f_phis(ind)
59      theta_rhodz=f_theta_rhodz(ind)
60      temp=f_temp(ind)
61      u=f_u(ind)
62      geopot=f_geopot(ind)
63      w=f_w(ind)
64      q=f_q(ind)
65
66      SELECT CASE(grid_type)
67      CASE(grid_ico)
68         CALL compute_etat0_collocated_hex(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q)
69      CASE(grid_unst)
70         CALL compute_etat0_collocated_unst(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q)
71      CASE DEFAULT
72         STOP 'Unexpected value of grid_type encountered in etat0_collocated.'
73      END SELECT
74
75    ENDDO
76   
77    CALL deallocate_field(f_temp)
78   
79  END SUBROUTINE etat0_collocated
80
81  SUBROUTINE compute_etat0_collocated_hex(phis,ps,mass,theta_rhodz,u, geopot,W, q)
82    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
83    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
84    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
85    REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
86    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
87    REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1)
88    REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1)
89    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
90
91    REAL(rstd) :: ps_e(3*iim*jjm)
92    REAL(rstd) :: phis_e(3*iim*jjm)
93    REAL(rstd) :: mass_e(3*iim*jjm,llm)
94    REAL(rstd) :: geopot_e(3*iim*jjm,llm+1)
95    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
96
97    w(:,:) = 0
98    CALL compute_etat0_collocated(iim*jjm  , lon_i, lat_i, phis,   ps,   mass,   theta_rhodz, geopot,   q)
99    CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u)
100
101  END SUBROUTINE compute_etat0_collocated_hex
102
103  SUBROUTINE compute_etat0_collocated_unst(phis,ps,mass,theta_rhodz,u, geopot,W, q)
104    REAL(rstd),INTENT(INOUT) :: ps(primal_num)
105    REAL(rstd),INTENT(INOUT) :: mass(llm,primal_num)
106    REAL(rstd),INTENT(OUT) :: theta_rhodz(llm,primal_num)
107    REAL(rstd),INTENT(OUT) :: phis(primal_num)
108    REAL(rstd),INTENT(OUT) :: u(llm, edge_num)
109    REAL(rstd),INTENT(OUT) :: W(llm+1, primal_num)
110    REAL(rstd),INTENT(OUT) :: geopot(llm+1, primal_num)
111    REAL(rstd),INTENT(OUT) :: q(llm, primal_num, nqtot)
112    ! The 2D/3D arrays below have the shape expected by compute_etat0_collocated
113    REAL(rstd) :: ps_e(edge_num)
114    REAL(rstd) :: phis_e(edge_num)
115    REAL(rstd) :: u_e(edge_num, llm)
116    REAL(rstd) :: ep(edge_num,3)
117    REAL(rstd) :: mass_i(primal_num, llm), theta_rhodz_i(primal_num, llm), mass_e(edge_num, llm)
118    REAL(rstd) :: geopot_i(edge_num, llm+1), geopot_e(edge_num, llm+1)
119    REAL(rstd) :: q_i(primal_num, llm, nqtot), q_e(edge_num, llm, nqtot)
120    INTEGER :: iq
121
122    w(:,:) = 0
123    ep = TRANSPOSE(ep_e)
124    CALL compute_etat0_collocated(primal_num  , lon_i, lat_i, phis,   ps,   mass_i,   theta_rhodz_i, geopot_i,   q_i)
125    CALL compute_etat0_collocated(edge_num, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep, u_e)
126    mass = TRANSPOSE(mass_i)
127    theta_rhodz = TRANSPOSE(theta_rhodz_i)
128    geopot = TRANSPOSE(geopot_i)
129    u = TRANSPOSE(u_e)
130    DO iq=1,nqtot
131       q(:,:,iq) = TRANSPOSE(q_i(:,:,iq))
132    END DO
133  END SUBROUTINE compute_etat0_collocated_unst
134 
135  SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, theta_rhodz, geopot, q, ep, uperp)
136    USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0
137    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
138    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
139    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0
140    USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
141    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
142    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
143    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0 
144    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
145    USE etat0_temperature_mod, ONLY: compute_temperature => compute_etat0
146    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0
147    USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0
148    USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0
149    USE disvert_mod, ONLY : ptop, ap, bp, mass_ak, mass_bk, mass_dak, mass_dbk
150    INTEGER :: ngrid
151    REAL(rstd),INTENT(IN)    :: lon(ngrid), lat(ngrid)
152    REAL(rstd),INTENT(INOUT) :: ps(ngrid)
153    REAL(rstd),INTENT(INOUT) :: mass(ngrid,llm)
154    REAL(rstd),INTENT(OUT)   :: theta_rhodz(ngrid,llm)
155    REAL(rstd),INTENT(OUT)   :: phis(ngrid)
156    REAL(rstd),INTENT(OUT)   :: geopot(ngrid,llm+1)
157    REAL(rstd),INTENT(OUT)   :: q(ngrid,llm,nqtot)
158    REAL(rstd),INTENT(OUT), OPTIONAL   :: ep(ngrid,3), uperp(ngrid,llm)
159
160    REAL(rstd)   :: ulon(ngrid,llm)
161    REAL(rstd)   :: ulat(ngrid,llm)
162    REAL(rstd)   :: temp(ngrid,llm)
163
164    LOGICAL :: autoinit_mass, autoinit_NH
165    INTEGER :: ij,l
166    REAL(rstd) :: clat, slat, clon, slon, u3d(3), p_ik, mass_ik, v_ik
167    REAL(rstd) :: q_ik, r_ik, chi, entropy, theta, log_p_preff
168
169    ! For NH, geopotential and vertical momentum must be initialized.
170    ! Unless autoinit_NH is set to .FALSE. , they will be initialized
171    ! to hydrostatic geopotential and zero
172    autoinit_mass = .TRUE.
173    autoinit_NH = .NOT. hydrostatic
174
175    SELECT CASE (TRIM(etat0_type))
176    CASE ('isothermal')
177       CALL compute_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
178    CASE ('temperature_profile')
179       CALL compute_temperature(ngrid, phis, ps, temp, ulon, ulat, q)
180    CASE('jablonowsky06')
181       CALL compute_jablonowsky06(ngrid, lon, lat, phis, ps, temp, ulon, ulat)
182    CASE('dcmip1')
183       CALL compute_dcmip1(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
184    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
185       CALL compute_dcmip2(ngrid, lon, lat, phis, ps, temp, ulon, ulat)
186    CASE('dcmip3')
187       CALL compute_dcmip3(ngrid, lon, lat, phis, ps, temp, ulon, ulat, geopot, q)
188       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot
189    CASE('dcmip4')
190       CALL compute_dcmip4(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
191    CASE('dcmip5')
192       CALL compute_dcmip5(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
193    CASE('bubble')
194       CALL compute_bubble(ngrid, lon, lat, phis, ps, temp, ulon, ulat, geopot, q)
195!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot
196    CASE('williamson91.6')
197       CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), theta_rhodz(:,1), ulon(:,1), ulat(:,1))
198       autoinit_mass = .FALSE. ! do not overwrite mass
199    CASE('dcmip2016_baroclinic_wave')
200       CALL compute_dcmip2016_baroclinic_wave(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
201    CASE('dcmip2016_cyclone')
202       CALL compute_dcmip2016_cyclone(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
203    CASE('dcmip2016_supercell')
204       CALL compute_dcmip2016_supercell(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
205    END SELECT
206
207    IF(PRESENT(uperp)) THEN
208       DO l = 1, llm
209          DO ij=1,ngrid
210             clat = COS(lat(ij))
211             slat = SIN(lat(ij))
212             clon = COS(lon(ij))
213             slon = SIN(lon(ij))
214             u3d(1) = -clon*slat*ulat(ij,l) - slon*ulon(ij,l)
215             u3d(2) = -slon*slat*ulat(ij,l) + clon*ulon(ij,l)
216             u3d(3) =       clat*ulat(ij,l)   
217             uperp(ij,l) = SUM(u3d*ep(ij,:))
218          END DO
219       END DO
220    END IF
221
222    IF(autoinit_mass) THEN
223       DO l = 1, llm
224          DO ij=1,ngrid
225             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
226             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
227             mass(ij,l) = mass_ik
228             SELECT CASE(caldyn_thermo)
229             CASE(thermo_theta)
230                theta = temp(ij,l)*(p_ik/preff)**(-kappa)
231                theta_rhodz(ij,l) = mass_ik * theta
232             CASE(thermo_entropy)
233                log_p_preff = log(p_ik/preff)
234                chi = log(temp(ij,l)/Treff)
235                entropy = cpp*chi-Rd*log_p_preff
236                theta_rhodz(ij,l) = mass_ik * entropy
237             CASE(thermo_moist)
238                q_ik=q(ij,l,1)
239                r_ik=1.-q_ik
240                mass_ik = mass_ik*(1.-q_ik) ! dry mass
241                log_p_preff = log(p_ik/preff)
242                chi = log(temp(ij,l)/Treff)
243                entropy = r_ik*(cpp*chi-Rd*nu) + q_ik*(cppv*chi-Rv*log_p_preff)
244                theta_rhodz(ij,l) = mass_ik * entropy
245             CASE DEFAULT
246                STOP
247             END SELECT
248          END DO
249       END DO
250    END IF
251
252    IF(autoinit_NH) THEN
253       geopot(:,1) = phis(:) ! surface geopotential
254       DO l = 1, llm
255          DO ij=1,ngrid
256             ! hybrid pressure coordinate
257             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
258             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
259             ! v=R.T/p, R=kappa*cpp
260             v_ik = kappa*cpp*temp(ij,l)/p_ik
261             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g
262          END DO
263       END DO
264    END IF
265
266  END SUBROUTINE compute_etat0_collocated
267
268END MODULE etat0_collocated_mod
Note: See TracBrowser for help on using the repository browser.