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

Last change on this file since 592 was 592, checked in by dubos, 7 years ago

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