source: codes/icosagcm/devel/src/diagnostics/observable.f90 @ 913

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

devel : compute_pression for unstructured mesh
File size: 14.2 KB
Line 
1MODULE observable_mod
2  USE icosa
3  USE caldyn_vars_mod
4  USE diagflux_mod
5  USE output_field_mod
6  IMPLICIT NONE
7  PRIVATE
8
9  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), &
10       f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
11       f_buf_ulon(:), f_buf_ulat(:)
12  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
13  TYPE(t_field),POINTER, SAVE :: f_pmid(:)
14
15! temporary shared variable for caldyn
16  TYPE(t_field),POINTER, SAVE :: f_theta(:)
17
18  PUBLIC init_observable, write_output_fields_basic, & 
19       f_theta, f_buf_i, f_buf_ulon, f_buf_ulat
20
21CONTAINS
22 
23  SUBROUTINE init_observable
24    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
25    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
26    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
27    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
28    CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
29    CALL allocate_field(f_buf_Fel, field_u,type_real,llm+1, name="buf_F_el")
30    CALL allocate_field(f_buf_v,   field_z,type_real,llm, name="buf_v")
31    CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
32
33    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta') ! potential temperature
34    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure
35  END SUBROUTINE init_observable
36
37  SUBROUTINE write_output_fields_basic(init, f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
38    USE xios_mod
39    USE disvert_mod
40    USE wind_mod
41    USE omp_para
42    USE time_mod
43    USE xios_mod
44    USE earth_const
45    USE compute_pression_mod, ONLY : pression_mid
46    USE vertical_interp_mod
47    USE theta2theta_rhodz_mod
48    USE omega_mod
49    USE kinetic_mod
50    USE grid_param
51    LOGICAL, INTENT(IN) :: init
52    INTEGER :: l
53    REAL :: scalar(1)
54    REAL :: mid_ap(llm)
55    REAL :: mid_bp(llm)
56
57    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
58!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
59
60    CALL transfert_request(f_ps,req_i1)
61   
62    IF(init) THEN
63       IF(is_master) PRINT *, 'Creating output files ...'
64       scalar(1)=dt
65       IF (is_omp_master) CALL xios_send_field("timestep", scalar)
66       scalar(1)=preff
67       IF (is_omp_master) CALL xios_send_field("preff", scalar)
68       IF (is_omp_master) CALL xios_send_field("ap",ap)
69       IF (is_omp_master) CALL xios_send_field("bp",bp)
70       DO l=1,llm
71          mid_ap(l)=(ap(l)+ap(l+1))/2
72          mid_bp(l)=(bp(l)+bp(l+1))/2
73       ENDDO
74       IF (is_omp_master) CALL xios_send_field("mid_ap",mid_ap)
75       IF (is_omp_master) CALL xios_send_field("mid_bp",mid_bp)
76
77       CALL output_field("phis",f_phis)
78       CALL output_field("Ai",geom%Ai) 
79       IF(is_master) PRINT *, '... done creating output files. Writing initial condition ...'
80    END IF
81
82    IF(nqdyn>1) THEN
83       CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i)
84       IF(init) THEN
85          CALL output_field("dyn_q_init",f_buf_i)
86       ELSE
87          CALL output_field("dyn_q",f_buf_i)
88       END IF
89    END IF
90
91    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)
92    IF(init) THEN
93       CALL output_field("theta_init",f_buf_i)
94    ELSE
95       CALL output_field("theta",f_buf_i)
96    END IF
97
98    CALL pression_mid(f_ps, f_pmid)
99    CALL diagnose_temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T
100
101    IF(init) THEN
102       CALL output_field("p_init",f_pmid)
103       CALL output_field("ps_init",f_ps)
104       CALL output_field("mass_init",f_mass)
105       CALL output_field("geopot_init",f_geopot)
106       CALL output_field("q_init",f_q)
107
108       CALL output_field("temp_init",f_buf_i)
109    ELSE
110       CALL output_field("p",f_pmid)
111       CALL output_field("ps",f_ps)
112       CALL output_field("mass",f_mass)
113       CALL output_field("geopot",f_geopot)
114       CALL output_field("q",f_q)
115
116       CALL output_field("temp",f_buf_i)
117       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
118       CALL output_field("t850",f_buf_s)
119       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
120       CALL output_field("t500",f_buf_s)
121       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,preff)
122       CALL output_field("SST",f_buf_s)       
123    END IF
124   
125    IF(grid_type == grid_unst) RETURN
126
127    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
128    CALL transfert_request(f_buf_uh,req_e1_vect) 
129    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
130    IF(init) THEN
131       CALL output_field("uz_init",f_buf_i)
132       CALL output_field("ulon_init",f_buf_ulon)
133       CALL output_field("ulat_init",f_buf_ulat)
134       IF(is_master) PRINT *, 'Done writing initial condition ...'
135    ELSE
136       CALL output_field("uz",f_buf_i)
137       CALL output_field("ulon",f_buf_ulon)
138       CALL output_field("ulat",f_buf_ulat)
139
140       !       CALL output_field("exner",f_pk)
141       !       CALL output_field("pv",f_qv)
142       
143       CALL vertical_interp(f_pmid,f_buf_ulon,f_buf_s,85000.)
144       CALL output_field("u850",f_buf_s)
145       CALL vertical_interp(f_pmid,f_buf_ulon,f_buf_s,50000.)
146       CALL output_field("u500",f_buf_s)
147       
148       CALL vertical_interp(f_pmid,f_buf_ulat,f_buf_s,85000.)
149       CALL output_field("v850",f_buf_s)
150       CALL vertical_interp(f_pmid,f_buf_ulat,f_buf_s,50000.)
151       CALL output_field("v500",f_buf_s)
152
153       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
154       CALL output_field("w850",f_buf_s)
155       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
156       CALL output_field("w500",f_buf_s)   
157
158       CALL w_omega(f_ps, f_u, f_buf_i)
159       CALL output_field("omega",f_buf_i)
160       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
161       CALL output_field("omega850",f_buf_s)
162       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
163       CALL output_field("omega500",f_buf_s)
164    END IF
165
166    CALL kinetic(f_u, f_buf_i)
167    IF(init) THEN
168       CALL output_field("kinetic_trisk_init",f_buf_i)
169    ELSE
170       CALL output_field("kinetic_trisk",f_buf_i)
171    END IF
172
173    CALL kinetic_new(f_u, f_buf_v, f_buf_i)
174    IF(init) THEN
175       CALL output_field("kinetic_init",f_buf_i)
176    ELSE
177       CALL output_field("kinetic",f_buf_i)
178    END IF
179
180    IF(.NOT. init) THEN
181       IF(diagflux_on) THEN
182          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat)
183          CALL output_field("mass_t", f_masst)
184          CALL output_field("massflux_lon",f_buf_ulon)
185          CALL output_field("massflux_lat",f_buf_ulat)
186
187          CALL output_energyflux(f_ulont, f_ulonfluxt, "ulon_t", "ulonflux_lon", "ulonflux_lat")
188          CALL output_energyflux(f_thetat, f_thetafluxt, "theta_t", "thetaflux_lon", "thetaflux_lat")
189          CALL output_energyflux(f_epot, f_epotfluxt, "epot_t", "epotflux_lon", "epotflux_lat")
190          CALL output_energyflux(f_ekin, f_ekinfluxt, "ekin_t", "ekinflux_lon", "ekinflux_lat")
191          CALL output_energyflux(f_enthalpy, f_enthalpyfluxt, "enthalpy_t", "enthalpyflux_lon", "enthalpyflux_lat")
192
193          CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
194          CALL output_field("qmass_t", f_qmasst)
195          CALL output_field("qflux_lon",f_qfluxt_lon)
196          CALL output_field("qflux_lat",f_qfluxt_lat)
197          CALL zero_qfluxt  ! restart accumulating fluxes
198       END IF
199    END IF
200  END SUBROUTINE write_output_fields_basic
201
202  SUBROUTINE output_energyflux(f_energy, f_flux, name_energy, name_fluxlon, name_fluxlat)
203    TYPE(t_field), POINTER :: f_energy(:), f_flux(:)
204    CHARACTER(*), INTENT(IN) :: name_energy, name_fluxlon, name_fluxlat
205    CALL transfert_request(f_flux,req_e1_vect)
206    CALL flux_centered_lonlat(1./(itau_out*dt) , f_flux, f_buf_ulon, f_buf_ulat)
207    CALL output_field(name_energy,  f_energy)
208    CALL output_field(name_fluxlon, f_buf_ulon)
209    CALL output_field(name_fluxlat, f_buf_ulat)
210  END SUBROUTINE output_energyflux
211 
212 !------------------- Conversion from prognostic to observable variables ------------------
213
214  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
215    USE disvert_mod, ONLY : caldyn_eta, eta_mass
216    USE compute_diagnostics_mod, ONLY : compute_rhodz
217    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
218         f_u(:), f_W(:), f_uz(:), &  ! IN
219         f_uh(:)                         ! OUT
220    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:), F_el(:,:)
221    INTEGER :: ind
222   
223    DO ind=1,ndomain
224       IF (.NOT. assigned_domain(ind)) CYCLE
225       CALL swap_dimensions(ind)
226       CALL swap_geometry(ind)
227       geopot = f_geopot(ind)
228       rhodz = f_rhodz(ind)
229       u = f_u(ind)
230       W = f_W(ind)
231       uh  = f_uh(ind)
232       F_el  = f_buf_Fel(ind)
233       IF(caldyn_eta==eta_mass) THEN
234          ps=f_ps(ind)
235          CALL compute_rhodz(.TRUE., ps, rhodz)
236       END IF
237       uz = f_uz(ind)
238       !$OMP BARRIER
239       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W, F_el, uh,uz)
240       !$OMP BARRIER
241    END DO
242  END SUBROUTINE progonostic_vel_to_horiz
243 
244  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, F_el, uh, uz)
245    USE omp_para
246    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
247    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
248    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
249    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
250    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
251    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
252    INTEGER :: ij,l
253    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
254    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
255    ! NB : u and uh are not in DEC form, they are normal components   
256    ! => we must divide by de
257    IF(hydrostatic) THEN
258       uh(:,:)=u(:,:)
259       uz(:,:)=0.
260    ELSE
261       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
262          DO ij=ij_begin_ext, ij_end_ext
263             ! Compute on edge 'right'
264             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
265             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
266             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
267             ! Compute on edge 'lup'
268             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
269             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
270             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
271             ! Compute on edge 'ldown'
272             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
273             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
274             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
275          END DO
276       END DO
277
278       ! We need a barrier here because we compute F_el above and do a vertical average below
279       !$OMP BARRIER
280
281       DO l=ll_begin, ll_end ! compute on k levels (full levels)
282          DO ij=ij_begin_ext, ij_end_ext
283             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
284             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
285             ! uh = u-w.grad(Phi) = u - uz.grad(z)
286             uh(ij+u_right,l) = u(ij+u_right,l) - (F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
287             uh(ij+u_lup,l)   = u(ij+u_lup,l)   - (F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))     / (rhodz(ij,l)+rhodz(ij+t_lup,l))
288             uh(ij+u_ldown,l) = u(ij+u_ldown,l) - (F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
289          END DO
290       END DO
291
292    END IF
293  END SUBROUTINE compute_prognostic_vel_to_horiz
294
295  SUBROUTINE diagnose_temperature(f_pmid,f_q,f_temp)
296    USE icosa
297    IMPLICIT NONE
298    TYPE(t_field), POINTER :: f_pmid(:)         ! IN
299    TYPE(t_field), POINTER :: f_q(:)            ! IN
300    TYPE(t_field), POINTER :: f_temp(:)         ! INOUT
301   
302    REAL(rstd), POINTER :: pmid(:,:)
303    REAL(rstd), POINTER :: q(:,:,:)
304    REAL(rstd), POINTER :: temp(:,:)
305    INTEGER :: ind
306   
307    DO ind=1,ndomain
308       IF (.NOT. assigned_domain(ind)) CYCLE
309       CALL swap_dimensions(ind)
310       CALL swap_geometry(ind)
311       pmid=f_pmid(ind)
312       q=f_q(ind)
313       temp=f_temp(ind)
314       CALL compute_diagnose_temp(pmid,q,temp)
315    END DO
316  END SUBROUTINE diagnose_temperature
317 
318  SUBROUTINE compute_diagnose_temp(pmid,q,temp)
319    USE omp_para
320    REAL(rstd),INTENT(IN)    :: pmid(iim*jjm,llm)
321    REAL(rstd),INTENT(IN)    :: q(iim*jjm,llm,nqtot)
322    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)
323
324    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix
325    INTEGER :: ij,l
326
327    Rd = kappa*cpp
328    DO l=ll_begin,ll_end
329       DO ij=ij_begin,ij_end
330          p_ik = pmid(ij,l)
331          theta_ik = temp(ij,l)
332          qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
333          SELECT CASE(caldyn_thermo)
334          CASE(thermo_theta)
335             temp_ik = theta_ik*((p_ik/preff)**kappa)
336          CASE(thermo_entropy)
337             temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
338          CASE(thermo_moist)
339             Rmix = Rd+qv*Rv
340             chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
341             temp_ik = Treff*exp(chi)
342          END SELECT
343          IF(physics_thermo==thermo_fake_moist) temp_ik=temp_ik/(1+0.608*qv)
344          temp(ij,l)=temp_ik
345       END DO
346    END DO
347  END SUBROUTINE compute_diagnose_temp
348
349  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
350    TYPE(t_field), POINTER :: f_TV(:)
351    TYPE(t_field), POINTER :: f_q(:)
352    TYPE(t_field), POINTER :: f_T(:)
353   
354    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
355    INTEGER :: ind
356   
357    DO ind=1,ndomain
358       IF (.NOT. assigned_domain(ind)) CYCLE
359       CALL swap_dimensions(ind)
360       CALL swap_geometry(ind)
361       Tv=f_Tv(ind)
362       T=f_T(ind)
363       SELECT CASE(physics_thermo)
364       CASE(thermo_dry)
365          T=Tv
366       CASE(thermo_fake_moist)
367          q=f_q(ind)
368          T=Tv/(1+0.608*q(:,:,1))
369       END SELECT
370    END DO
371  END SUBROUTINE Tv2T
372
373  SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
374    INTEGER, INTENT(IN) :: iq
375    TYPE(t_field), POINTER :: f_mass(:), f_theta_rhodz(:), f_theta(:)
376    REAL(rstd), POINTER :: mass(:,:), theta_rhodz(:,:,:), theta(:,:)
377    INTEGER :: ind
378    DO ind=1,ndomain
379       IF (.NOT. assigned_domain(ind)) CYCLE
380       CALL swap_dimensions(ind)
381       CALL swap_geometry(ind)
382       mass=f_mass(ind)
383       theta_rhodz=f_theta_rhodz(ind)
384       theta=f_theta(ind)
385       theta(:,:) = theta_rhodz(:,:,iq) / mass(:,:)
386    END DO
387  END SUBROUTINE divide_by_mass
388   
389END MODULE observable_mod
Note: See TracBrowser for help on using the repository browser.