source: dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0_database.f90 @ 3887

Last change on this file since 3887 was 3887, checked in by ymipsl, 9 years ago
  • Initial state from file in dynamico can be properly generated by switching flag in run_icosa.def
  • bug fix for openMP when generate initial state from file

YM
File size: 5.2 KB
Line 
1MODULE etat0_database_mod
2
3
4CONTAINS
5
6  SUBROUTINE init_etat0_database
7  USE xios
8  IMPLICIT NONE
9 
10    CALL xios_set_fieldgroup_attr("read_fields",enabled=.TRUE.)
11    CALL xios_set_filegroup_attr("read_files",enabled=.TRUE.)
12
13  END SUBROUTINE init_etat0_database
14
15  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
16  USE icosa
17  USE restart_mod
18  USE wind_mod
19  USE write_field_mod
20  USE time_mod
21  USE transfert_mod
22  USE xios_mod
23  USE write_field_mod
24  USE vertical_remap_mod
25  USE theta2theta_rhodz_mod
26  USE qsat_mod
27  USE pression_mod
28  USE write_etat0_mod
29  USE omp_para
30  IMPLICIT NONE
31    TYPE(t_field),POINTER :: f_ps(:)
32    TYPE(t_field),POINTER :: f_phis(:)
33    TYPE(t_field),POINTER :: f_theta_rhodz(:)
34    TYPE(t_field),POINTER :: f_u(:)
35    TYPE(t_field),POINTER :: f_q(:)
36
37    TYPE(t_field),POINTER,SAVE :: f_ulon_reg(:)
38    TYPE(t_field),POINTER,SAVE :: f_ulat_reg(:)
39    TYPE(t_field),POINTER,SAVE :: f_temp_reg(:)
40    TYPE(t_field),POINTER,SAVE :: f_q_reg(:)
41
42    TYPE(t_field),POINTER,SAVE :: f_ts(:)
43    TYPE(t_field),POINTER,SAVE :: f_z(:)
44    TYPE(t_field),POINTER,SAVE :: f_ulon(:)
45    TYPE(t_field),POINTER,SAVE :: f_ulat(:)
46    TYPE(t_field),POINTER,SAVE :: f_temp(:)
47    TYPE(t_field),POINTER,SAVE :: f_q1(:)
48    TYPE(t_field),POINTER,SAVE :: f_qsat(:)
49    TYPE(t_field),POINTER,SAVE :: f_p(:)
50    INTEGER :: nb_level
51    REAL,ALLOCATABLE:: levels(:)
52    INTEGER :: ind
53
54    CALL xios_read_field("relief",f_phis)
55   
56    CALL writeField("relief_out",f_phis,once=.TRUE.)
57
58    DO ind=1,ndomain
59      IF (.NOT. assigned_domain(ind)) CYCLE
60      f_phis(ind)%rval2d(:)=f_phis(ind)%rval2d(:)*g     
61    ENDDO
62
63
64    IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn",n_glo=nb_level)
65    CALL bcast_omp(nb_level)
66    ALLOCATE(levels(nb_level))
67
68    IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn",value=levels)
69    CALL bcast_omp(levels)
70   
71    levels=levels*100  ! hectoPascal -> Pascal
72 
73    CALL allocate_field(f_ts, field_t, type_real, name="ts")
74    CALL allocate_field(f_z, field_t, type_real, name="z")
75    CALL allocate_field(f_ulon_reg, field_t, type_real,nb_level)
76    CALL allocate_field(f_ulat_reg, field_t, type_real,nb_level)
77    CALL allocate_field(f_temp_reg, field_t, type_real,nb_level)
78    CALL allocate_field(f_q_reg,    field_t, type_real,nb_level)
79   
80    CALL allocate_field(f_q1,    field_t, type_real,llm)
81    CALL allocate_field(f_qsat,  field_t, type_real,llm)
82    CALL allocate_field(f_p,  field_t, type_real,llm+1)
83    CALL allocate_field(f_temp,  field_t, type_real,llm)
84    CALL allocate_field(f_ulon,  field_t, type_real,llm)
85    CALL allocate_field(f_ulat,  field_t, type_real,llm)
86
87    CALL xios_read_field("z",f_z)
88    CALL xios_read_field("ps",f_ps)
89    CALL xios_read_field("ts",f_ts)
90    CALL writeField("ps_out",f_ps)
91
92!$OMP BARRIER
93   
94!    CALL writeField("phis_out",f_phis,once=.TRUE.)
95!    CALL writeField("ts_out",f_ts,once=.TRUE.)
96
97! make correction to ps due to relief at higher resolution
98! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY
99    DO ind=1,ndomain
100      IF (.NOT. assigned_domain(ind)) CYCLE
101      f_ps(ind)%rval2d(:)=f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:)-f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:))
102    ENDDO
103    CALL transfert_request(f_ps,req_i0)
104    CALL writeField("ps_out",f_ps)
105   
106   
107
108    CALL xios_read_field("temp",f_temp_reg)
109    CALL vertical_remap(levels,f_temp_reg,f_ps,f_temp)
110    CALL transfert_request(f_temp,req_i0)
111    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
112
113    CALL xios_read_field("u",f_ulon_reg)
114    CALL vertical_remap(levels,f_ulon_reg,f_ps,f_ulon)
115    CALL xios_read_field("v",f_ulat_reg)
116    CALL vertical_remap(levels,f_ulat_reg,f_ps,f_ulat)
117    CALL transfert_request(f_ulat,req_i0)
118    CALL transfert_request(f_ulon,req_i0)
119    CALL ulonlat2un(f_ulon, f_ulat,f_u)
120
121    CALL xios_read_field("q",f_q_reg)
122    CALL vertical_remap(levels,f_q_reg,f_ps,f_q1)
123
124    CALL pression(f_ps,f_p)
125! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2   
126    CALL qsat(f_temp,f_p,f_qsat)
127    CALL transfert_request(f_qsat,req_i0)
128
129    DO ind=1,ndomain
130      IF (.NOT. assigned_domain(ind)) CYCLE
131      f_q(ind)%rval4d(:,:,:)=0
132      f_q(ind)%rval4d(:,:,1)=f_q1(ind)%rval3d(:,:)*f_qsat(ind)%rval3d(:,:)*0.01
133      WHERE(f_q(ind)%rval4d(:,:,1)<0) f_q(ind)%rval4d(:,:,1)=0
134    ENDDO
135
136   
137    CALL writeField("ulon_out",f_ulon,once=.TRUE.)
138    CALL writeField("ulat_out",f_ulat,once=.TRUE.)
139    CALL writeField("temp_out",f_temp,once=.TRUE.)
140    CALL writeField("temp2_out",f_theta_rhodz,once=.TRUE.)
141    CALL writeField("f_qsat",f_qsat,once=.TRUE.)
142    CALL writeField("f_q1",f_q1,once=.TRUE.)
143    CALL writeField("f_q",f_q,once=.TRUE.)
144    CALL write_etat0(0,f_ps,f_phis,f_theta_rhodz,f_u, f_q)
145   
146
147
148    CALL deallocate_field(f_ts)
149    CALL deallocate_field(f_z)
150    CALL deallocate_field(f_ulon_reg)
151    CALL deallocate_field(f_ulat_reg)
152    CALL deallocate_field(f_temp_reg)
153    CALL deallocate_field(f_q_reg)
154   
155    CALL deallocate_field(f_q1)
156    CALL deallocate_field(f_qsat)
157    CALL deallocate_field(f_p)
158    CALL deallocate_field(f_temp)
159    CALL deallocate_field(f_ulon)
160    CALL deallocate_field(f_ulat)   
161 
162  END SUBROUTINE etat0
163
164END MODULE etat0_database_mod
Note: See TracBrowser for help on using the repository browser.