source: codes/icosagcm/trunk/src/observable.f90 @ 434

Last change on this file since 434 was 434, checked in by dubos, 8 years ago

Dry/moist output cleanup
File size: 12.0 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_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_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
39    USE disvert_mod
40    USE wind_mod
41    USE output_field_mod
42    USE omp_para
43    USE time_mod
44    USE xios
45    USE earth_const
46    USE pression_mod
47    USE vertical_interp_mod
48    USE theta2theta_rhodz_mod
49    USE omega_mod
50    LOGICAL, INTENT(IN) :: init
51    INTEGER :: l
52    REAL :: scalar(1)
53    REAL :: mid_ap(llm)
54    REAL :: mid_bp(llm)
55
56    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
57!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
58
59    CALL transfert_request(f_ps,req_i1)
60   
61    IF(init) THEN
62       scalar(1)=dt
63       CALL xios_send_field("timestep", scalar)
64       scalar(1)=preff
65       CALL xios_send_field("preff", scalar)
66       CALL xios_send_field("ap",ap)
67       CALL xios_send_field("bp",bp)
68       DO l=1,llm
69          mid_ap(l)=(ap(l)+ap(l+1))/2
70          mid_bp(l)=(bp(l)+bp(l+1))/2
71       ENDDO
72       CALL xios_send_field("mid_ap",mid_ap)
73       CALL xios_send_field("mid_bp",mid_bp)
74
75       CALL output_field("phis",f_phis)
76       CALL output_field("Ai",geom%Ai)       
77    END IF
78
79    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)
80    IF(init) THEN
81       CALL output_field("theta_init",f_buf_i)
82    ELSE
83       CALL output_field("theta",f_buf_i)
84    END IF
85
86    IF(nqdyn>1) THEN
87       CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i)
88       IF(init) THEN
89          CALL output_field("dyn_q_init",f_buf_i)
90       ELSE
91          CALL output_field("dyn_q",f_buf_i)
92       END IF
93    END IF
94
95!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i)
96!    CALL Tv2T(f_buf_i,f_q,f_buf1_i)
97    CALL diagnose_temperature(f_ps, f_theta_rhodz, f_q, f_buf_i)
98
99    IF(init) THEN
100       CALL output_field("temp_init",f_buf_i)
101    ELSE
102       CALL output_field("temp",f_buf_i)
103       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
104       CALL output_field("t850",f_buf_s)
105       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
106       CALL output_field("t500",f_buf_s)
107       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,preff)
108       CALL output_field("SST",f_buf_s)       
109    END IF
110   
111    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
112    CALL transfert_request(f_buf_uh,req_e1_vect) 
113    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
114    CALL pression_mid(f_ps, f_pmid)
115    IF(init) THEN
116       CALL output_field("uz_init",f_buf_i)
117       CALL output_field("ulon_init",f_buf_ulon)
118       CALL output_field("ulat_init",f_buf_ulat)
119       CALL output_field("p_init",f_pmid)
120       CALL output_field("ps_init",f_ps)
121       CALL output_field("mass_init",f_mass)
122       CALL output_field("geopot_init",f_geopot)
123       CALL output_field("q_init",f_q)
124    ELSE
125       CALL output_field("uz",f_buf_i)
126       CALL output_field("ulon",f_buf_ulon)
127       CALL output_field("ulat",f_buf_ulat)
128       CALL output_field("p",f_pmid)
129       CALL output_field("ps",f_ps)
130       CALL output_field("mass",f_mass)
131       CALL output_field("geopot",f_geopot)
132       CALL output_field("q",f_q)
133
134       !       CALL output_field("exner",f_pk)
135       !       CALL output_field("pv",f_qv)
136       
137       CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,85000.)
138       CALL output_field("u850",f_buf_s)
139       CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,50000.)
140       CALL output_field("u500",f_buf_s)
141       
142       CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,85000.)
143       CALL output_field("v850",f_buf_s)
144       CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,50000.)
145       CALL output_field("v500",f_buf_s)
146
147       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
148       CALL output_field("w850",f_buf_s)
149       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
150       CALL output_field("w500",f_buf_s)   
151
152       CALL w_omega(f_ps, f_u, f_buf_i)
153       CALL output_field("omega",f_buf_i)
154       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
155       CALL output_field("omega850",f_buf_s)
156       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
157       CALL output_field("omega500",f_buf_s)
158    END IF
159
160  END SUBROUTINE write_output_fields_basic
161 
162 !------------------- Conversion from prognostic to observable variables ------------------
163
164  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
165    USE disvert_mod
166    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
167         f_u(:), f_W(:), f_uz(:), &  ! IN
168         f_uh(:)                         ! OUT
169    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:)
170    INTEGER :: ind
171   
172    DO ind=1,ndomain
173       IF (.NOT. assigned_domain(ind)) CYCLE
174       CALL swap_dimensions(ind)
175       CALL swap_geometry(ind)
176       geopot = f_geopot(ind)
177       rhodz = f_rhodz(ind)
178       u = f_u(ind)
179       W = f_W(ind)
180       uh  = f_uh(ind)
181       IF(caldyn_eta==eta_mass) THEN
182          ps=f_ps(ind)
183          CALL compute_rhodz(.TRUE., ps, rhodz)
184       END IF
185       uz = f_uz(ind)
186       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W,uh,uz)
187    END DO
188  END SUBROUTINE progonostic_vel_to_horiz
189 
190  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, uh, uz)
191    USE omp_para
192    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
193    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
194    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
195    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
196    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
197    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
198    INTEGER :: ij,l
199    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
200    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
201    ! NB : u and uh are not in DEC form, they are normal components   
202    ! => we must divide by de
203    IF(hydrostatic) THEN
204       uh(:,:)=u(:,:)
205       uz(:,:)=0.
206    ELSE
207       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
208          DO ij=ij_begin_ext, ij_end_ext
209             ! Compute on edge 'right'
210             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
211             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
212             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
213             ! Compute on edge 'lup'
214             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
215             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
216             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
217             ! Compute on edge 'ldown'
218             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
219             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
220             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
221          END DO
222       END DO
223
224       DO l=ll_begin, ll_end ! compute on k levels (full levels)
225          DO ij=ij_begin_ext, ij_end_ext
226             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
227             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
228             ! uh = u-w.grad(Phi) = u - uz.grad(z)
229             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))
230             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))
231             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))
232          END DO
233       END DO
234
235    END IF
236  END SUBROUTINE compute_prognostic_vel_to_horiz
237
238  SUBROUTINE diagnose_temperature(f_ps,f_theta_rhodz,f_q,f_temp)
239    USE icosa
240    USE pression_mod
241    IMPLICIT NONE
242    TYPE(t_field), POINTER :: f_ps(:)           ! IN
243    TYPE(t_field), POINTER :: f_theta_rhodz(:)  ! IN
244    TYPE(t_field), POINTER :: f_q(:)            ! IN
245    TYPE(t_field), POINTER :: f_temp(:)         ! OUT
246   
247    REAL(rstd), POINTER :: ps(:)
248    REAL(rstd), POINTER :: theta_rhodz(:,:,:)
249    REAL(rstd), POINTER :: q(:,:,:)
250    REAL(rstd), POINTER :: temp(:,:)
251    INTEGER :: ind
252   
253    DO ind=1,ndomain
254       IF (.NOT. assigned_domain(ind)) CYCLE
255       CALL swap_dimensions(ind)
256       CALL swap_geometry(ind)
257       ps=f_ps(ind)
258       theta_rhodz=f_theta_rhodz(ind)
259       q=f_q(ind)
260       temp=f_temp(ind)
261       CALL compute_diagnose_temp(ps,theta_rhodz,q,temp)
262    END DO
263  END SUBROUTINE diagnose_temperature
264 
265  SUBROUTINE compute_diagnose_temp(ps,theta_rhodz,q,temp)
266    USE omp_para
267    USE pression_mod
268    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
269    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm,nqdyn)
270    REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot)
271    REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm)
272
273    REAL(rstd)  :: p(iim*jjm,llm+1)
274    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix
275    INTEGER :: ij,l
276
277    Rd = kappa*cpp
278    CALL compute_pression(ps,p,0)
279    DO l=ll_begin,ll_end
280       DO ij=ij_begin,ij_end
281          p_ik = .5*(p(ij,l)+p(ij,l+1))
282          theta_ik = g*theta_rhodz(ij,l,1)/(p(ij,l)-p(ij,l+1))
283          qv = q(ij,l,1) ! water vaper mixing ratio = mv/md
284          SELECT CASE(caldyn_thermo)
285          CASE(thermo_theta)
286             temp_ik = theta_ik*((p_ik/preff)**kappa)
287          CASE(thermo_entropy)
288             temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
289          CASE(thermo_moist)
290             Rmix = Rd+qv*Rv
291             chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
292             temp_ik = Treff*exp(chi)
293          END SELECT
294          IF(physics_thermo==thermo_fake_moist) temp_ik=temp_ik/(1+0.608*qv)
295          temp(ij,l)=temp_ik
296       END DO
297    END DO
298  END SUBROUTINE compute_diagnose_temp
299
300  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
301    TYPE(t_field), POINTER :: f_TV(:)
302    TYPE(t_field), POINTER :: f_q(:)
303    TYPE(t_field), POINTER :: f_T(:)
304   
305    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
306    INTEGER :: ind
307   
308    DO ind=1,ndomain
309       IF (.NOT. assigned_domain(ind)) CYCLE
310       CALL swap_dimensions(ind)
311       CALL swap_geometry(ind)
312       Tv=f_Tv(ind)
313       T=f_T(ind)
314       SELECT CASE(physics_thermo)
315       CASE(thermo_dry)
316          T=Tv
317       CASE(thermo_fake_moist)
318          q=f_q(ind)
319          T=Tv/(1+0.608*q(:,:,1))
320       END SELECT
321    END DO
322  END SUBROUTINE Tv2T
323
324  SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
325    INTEGER, INTENT(IN) :: iq
326    TYPE(t_field), POINTER :: f_mass(:), f_theta_rhodz(:), f_theta(:)
327    REAL(rstd), POINTER :: mass(:,:), theta_rhodz(:,:,:), theta(:,:)
328    INTEGER :: ind
329    DO ind=1,ndomain
330       IF (.NOT. assigned_domain(ind)) CYCLE
331       CALL swap_dimensions(ind)
332       CALL swap_geometry(ind)
333       mass=f_mass(ind)
334       theta_rhodz=f_theta_rhodz(ind)
335       theta=f_theta(ind)
336       theta(:,:) = theta_rhodz(:,:,iq) / mass(:,:)
337    END DO
338  END SUBROUTINE divide_by_mass
339   
340END MODULE observable_mod
Note: See TracBrowser for help on using the repository browser.