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

Last change on this file was 1052, checked in by dubos, 4 years ago

devel : diagnose divergence and vorticity
File size: 8.9 KB
RevLine 
[354]1MODULE observable_mod
2  USE icosa
[732]3  USE caldyn_vars_mod
[603]4  USE diagflux_mod
5  USE output_field_mod
[354]6  IMPLICIT NONE
7  PRIVATE
8
[374]9  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), &
[579]10       f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
[598]11       f_buf_ulon(:), f_buf_ulat(:)
[354]12  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
[397]13  TYPE(t_field),POINTER, SAVE :: f_pmid(:)
[354]14
15! temporary shared variable for caldyn
16  TYPE(t_field),POINTER, SAVE :: f_theta(:)
17
[714]18  PUBLIC init_observable, write_output_fields_basic, & 
19       f_theta, f_buf_i, f_buf_ulon, f_buf_ulat
[413]20
[354]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")
[374]28    CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
[579]29    CALL allocate_field(f_buf_Fel, field_u,type_real,llm+1, name="buf_F_el")
[374]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")
[354]32
[404]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
[354]35  END SUBROUTINE init_observable
[413]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)
[482]38    USE xios_mod
[413]39    USE disvert_mod
[354]40    USE wind_mod
41    USE omp_para
[397]42    USE time_mod
[482]43    USE xios_mod
[397]44    USE earth_const
[955]45    USE compute_pression_mod, ONLY : pression_mid, hydrostatic_pressure
[914]46    USE compute_temperature_mod
47    USE compute_velocity_mod
[1052]48    USE compute_vorticity_mod
49    USE compute_divergence_mod
[397]50    USE vertical_interp_mod
51    USE theta2theta_rhodz_mod
[919]52    USE compute_omega_mod, ONLY : w_omega
[728]53    USE kinetic_mod
[868]54    USE grid_param
[413]55    LOGICAL, INTENT(IN) :: init
56    INTEGER :: l
[397]57    REAL :: scalar(1)
58    REAL :: mid_ap(llm)
59    REAL :: mid_bp(llm)
60
[413]61    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
62!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
[403]63
[417]64    CALL transfert_request(f_ps,req_i1)
65   
[413]66    IF(init) THEN
[579]67       IF(is_master) PRINT *, 'Creating output files ...'
[413]68       scalar(1)=dt
[470]69       IF (is_omp_master) CALL xios_send_field("timestep", scalar)
[413]70       scalar(1)=preff
[470]71       IF (is_omp_master) CALL xios_send_field("preff", scalar)
72       IF (is_omp_master) CALL xios_send_field("ap",ap)
73       IF (is_omp_master) CALL xios_send_field("bp",bp)
[413]74       DO l=1,llm
75          mid_ap(l)=(ap(l)+ap(l+1))/2
76          mid_bp(l)=(bp(l)+bp(l+1))/2
77       ENDDO
[470]78       IF (is_omp_master) CALL xios_send_field("mid_ap",mid_ap)
79       IF (is_omp_master) CALL xios_send_field("mid_bp",mid_bp)
[413]80
81       CALL output_field("phis",f_phis)
[579]82       CALL output_field("Ai",geom%Ai) 
83       IF(is_master) PRINT *, '... done creating output files. Writing initial condition ...'
[413]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
[529]95    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)
96    IF(init) THEN
97       CALL output_field("theta_init",f_buf_i)
98    ELSE
99       CALL output_field("theta",f_buf_i)
100    END IF
[434]101
[955]102    CALL hydrostatic_pressure(f_mass, f_theta_rhodz, f_ps, f_pmid) 
[914]103    CALL temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T
[529]104
[413]105    IF(init) THEN
[906]106       CALL output_field("p_init",f_pmid)
107       CALL output_field("ps_init",f_ps)
108       CALL output_field("mass_init",f_mass)
109       CALL output_field("geopot_init",f_geopot)
110       CALL output_field("q_init",f_q)
111
[413]112       CALL output_field("temp_init",f_buf_i)
[1052]113       CALL output_field("vort_init",f_buf_v)
[413]114    ELSE
[906]115       CALL output_field("p",f_pmid)
116       CALL output_field("ps",f_ps)
117       CALL output_field("mass",f_mass)
118       CALL output_field("geopot",f_geopot)
119       CALL output_field("q",f_q)
120
[413]121       CALL output_field("temp",f_buf_i)
[436]122       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
[413]123       CALL output_field("t850",f_buf_s)
[436]124       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
[413]125       CALL output_field("t500",f_buf_s)
[436]126       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,preff)
[413]127       CALL output_field("SST",f_buf_s)       
128    END IF
[1052]129
130    CALL vorticity(f_u, f_buf_v)
131    CALL divergence(f_u, f_buf_i)
132    IF(init) THEN
133       CALL output_field("vort_init",f_buf_v)
134       CALL output_field("div_init",f_buf_i)
135    ELSE
136       CALL output_field("vort",f_buf_v)
137       CALL output_field("div",f_buf_i)
138    END IF
[397]139   
[906]140    IF(grid_type == grid_unst) RETURN
141
[914]142    CALL velocity(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_Fel, f_buf_uh, f_buf_i)
[374]143    CALL transfert_request(f_buf_uh,req_e1_vect) 
144    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
[413]145    IF(init) THEN
146       CALL output_field("uz_init",f_buf_i)
147       CALL output_field("ulon_init",f_buf_ulon)
148       CALL output_field("ulat_init",f_buf_ulat)
[579]149       IF(is_master) PRINT *, 'Done writing initial condition ...'
[413]150    ELSE
151       CALL output_field("uz",f_buf_i)
152       CALL output_field("ulon",f_buf_ulon)
153       CALL output_field("ulat",f_buf_ulat)
[397]154
[413]155       !       CALL output_field("exner",f_pk)
156       !       CALL output_field("pv",f_qv)
157       
[436]158       CALL vertical_interp(f_pmid,f_buf_ulon,f_buf_s,85000.)
[413]159       CALL output_field("u850",f_buf_s)
[436]160       CALL vertical_interp(f_pmid,f_buf_ulon,f_buf_s,50000.)
[413]161       CALL output_field("u500",f_buf_s)
162       
[436]163       CALL vertical_interp(f_pmid,f_buf_ulat,f_buf_s,85000.)
[413]164       CALL output_field("v850",f_buf_s)
[436]165       CALL vertical_interp(f_pmid,f_buf_ulat,f_buf_s,50000.)
[413]166       CALL output_field("v500",f_buf_s)
[397]167
[436]168       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
[413]169       CALL output_field("w850",f_buf_s)
[436]170       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
[413]171       CALL output_field("w500",f_buf_s)   
[397]172
[413]173       CALL w_omega(f_ps, f_u, f_buf_i)
174       CALL output_field("omega",f_buf_i)
[436]175       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
[413]176       CALL output_field("omega850",f_buf_s)
[436]177       CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,50000.)
[413]178       CALL output_field("omega500",f_buf_s)
179    END IF
[397]180
[728]181    CALL kinetic(f_u, f_buf_i)
182    IF(init) THEN
183       CALL output_field("kinetic_trisk_init",f_buf_i)
184    ELSE
185       CALL output_field("kinetic_trisk",f_buf_i)
186    END IF
187
188    CALL kinetic_new(f_u, f_buf_v, f_buf_i)
189    IF(init) THEN
190       CALL output_field("kinetic_init",f_buf_i)
191    ELSE
192       CALL output_field("kinetic",f_buf_i)
193    END IF
194
[588]195    IF(.NOT. init) THEN
196       IF(diagflux_on) THEN
[594]197          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat)
[592]198          CALL output_field("mass_t", f_masst)
199          CALL output_field("massflux_lon",f_buf_ulon)
200          CALL output_field("massflux_lat",f_buf_ulat)
[595]201
[603]202          CALL output_energyflux(f_ulont, f_ulonfluxt, "ulon_t", "ulonflux_lon", "ulonflux_lat")
203          CALL output_energyflux(f_thetat, f_thetafluxt, "theta_t", "thetaflux_lon", "thetaflux_lat")
204          CALL output_energyflux(f_epot, f_epotfluxt, "epot_t", "epotflux_lon", "epotflux_lat")
205          CALL output_energyflux(f_ekin, f_ekinfluxt, "ekin_t", "ekinflux_lon", "ekinflux_lat")
206          CALL output_energyflux(f_enthalpy, f_enthalpyfluxt, "enthalpy_t", "enthalpyflux_lon", "enthalpyflux_lat")
[595]207
[594]208          CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
209          CALL output_field("qmass_t", f_qmasst)
[588]210          CALL output_field("qflux_lon",f_qfluxt_lon)
211          CALL output_field("qflux_lat",f_qfluxt_lat)
212          CALL zero_qfluxt  ! restart accumulating fluxes
213       END IF
214    END IF
[354]215  END SUBROUTINE write_output_fields_basic
[603]216
217  SUBROUTINE output_energyflux(f_energy, f_flux, name_energy, name_fluxlon, name_fluxlat)
218    TYPE(t_field), POINTER :: f_energy(:), f_flux(:)
219    CHARACTER(*), INTENT(IN) :: name_energy, name_fluxlon, name_fluxlat
220    CALL transfert_request(f_flux,req_e1_vect)
221    CALL flux_centered_lonlat(1./(itau_out*dt) , f_flux, f_buf_ulon, f_buf_ulat)
222    CALL output_field(name_energy,  f_energy)
223    CALL output_field(name_fluxlon, f_buf_ulon)
224    CALL output_field(name_fluxlat, f_buf_ulat)
225  END SUBROUTINE output_energyflux
[354]226 
[413]227  SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
228    INTEGER, INTENT(IN) :: iq
229    TYPE(t_field), POINTER :: f_mass(:), f_theta_rhodz(:), f_theta(:)
230    REAL(rstd), POINTER :: mass(:,:), theta_rhodz(:,:,:), theta(:,:)
231    INTEGER :: ind
232    DO ind=1,ndomain
233       IF (.NOT. assigned_domain(ind)) CYCLE
234       CALL swap_dimensions(ind)
235       CALL swap_geometry(ind)
236       mass=f_mass(ind)
237       theta_rhodz=f_theta_rhodz(ind)
238       theta=f_theta(ind)
239       theta(:,:) = theta_rhodz(:,:,iq) / mass(:,:)
240    END DO
241  END SUBROUTINE divide_by_mass
242   
[354]243END MODULE observable_mod
Note: See TracBrowser for help on using the repository browser.