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

Last change on this file since 714 was 714, checked in by dubos, 6 years ago

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