1 | ! |
---|
2 | ! $Id: 1D_read_forc_cases.h 2720 2016-11-30 12:28:41Z asima $ |
---|
3 | ! |
---|
4 | !---------------------------------------------------------------------- |
---|
5 | ! forcing_les = .T. : Impose a constant cooling |
---|
6 | ! forcing_radconv = .T. : Pure radiative-convective equilibrium: |
---|
7 | !---------------------------------------------------------------------- |
---|
8 | |
---|
9 | |
---|
10 | nq1=0 |
---|
11 | nq2=0 |
---|
12 | |
---|
13 | if (forcing_les .or. forcing_radconv & |
---|
14 | & .or. forcing_GCSSold .or. forcing_fire) then |
---|
15 | |
---|
16 | if (forcing_fire) then |
---|
17 | !---------------------------------------------------------------------- |
---|
18 | !read fire forcings from fire.nc |
---|
19 | !---------------------------------------------------------------------- |
---|
20 | fich_fire='fire.nc' |
---|
21 | call read_fire(fich_fire,nlev_fire,nt_fire & |
---|
22 | & ,height,tttprof,qtprof,uprof,vprof,e12prof & |
---|
23 | & ,ugprof,vgprof,wfls,dqtdxls & |
---|
24 | & ,dqtdyls,dqtdtls,thlpcar) |
---|
25 | write(*,*) 'Forcing FIRE lu' |
---|
26 | kmax=120 ! nombre de niveaux dans les profils et forcages |
---|
27 | else |
---|
28 | !---------------------------------------------------------------------- |
---|
29 | ! Read profiles from files: prof.inp.001 and lscale.inp.001 |
---|
30 | ! (repris de readlesfiles) |
---|
31 | !---------------------------------------------------------------------- |
---|
32 | |
---|
33 | call readprofiles(nlev_max,kmax,nqtot,height, & |
---|
34 | & tttprof,qtprof,uprof,vprof, & |
---|
35 | & e12prof,ugprof,vgprof, & |
---|
36 | & wfls,dqtdxls,dqtdyls,dqtdtls, & |
---|
37 | & thlpcar,qprof,nq1,nq2) |
---|
38 | endif |
---|
39 | |
---|
40 | ! compute altitudes of play levels. |
---|
41 | zlay(1) =zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) |
---|
42 | do l = 2,llm |
---|
43 | zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf) |
---|
44 | enddo |
---|
45 | |
---|
46 | !---------------------------------------------------------------------- |
---|
47 | ! Interpolation of the profiles given on the input file to |
---|
48 | ! model levels |
---|
49 | !---------------------------------------------------------------------- |
---|
50 | zlay(1) = zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) |
---|
51 | do l=1,llm |
---|
52 | ! Above the max altutide of the input file |
---|
53 | |
---|
54 | if (zlay(l)<height(kmax)) mxcalc=l |
---|
55 | |
---|
56 | frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1)) |
---|
57 | ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1)) |
---|
58 | if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile |
---|
59 | temp(l) = ttt*(play(l)/pzero)**rkappa |
---|
60 | teta(l) = ttt |
---|
61 | else |
---|
62 | temp(l) = ttt |
---|
63 | teta(l) = ttt*(pzero/play(l))**rkappa |
---|
64 | endif |
---|
65 | print *,' temp,teta ',l,temp(l),teta(l) |
---|
66 | q(l,1) = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1)) |
---|
67 | u(l) = uprof(kmax)-frac*( uprof(kmax)- uprof(kmax-1)) |
---|
68 | v(l) = vprof(kmax)-frac*( vprof(kmax)- vprof(kmax-1)) |
---|
69 | ug(l) = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1)) |
---|
70 | vg(l) = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1)) |
---|
71 | IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2) & |
---|
72 | & -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2)) |
---|
73 | omega(l)= wfls(kmax)-frac*( wfls(kmax)- wfls(kmax-1)) |
---|
74 | |
---|
75 | dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1)) |
---|
76 | dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1)) |
---|
77 | do k=2,kmax |
---|
78 | frac = (height(k)-zlay(l))/(height(k)-height(k-1)) |
---|
79 | if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k) |
---|
80 | if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then |
---|
81 | ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1)) |
---|
82 | if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile |
---|
83 | temp(l) = ttt*(play(l)/pzero)**rkappa |
---|
84 | teta(l) = ttt |
---|
85 | else |
---|
86 | temp(l) = ttt |
---|
87 | teta(l) = ttt*(pzero/play(l))**rkappa |
---|
88 | endif |
---|
89 | print *,' temp,teta ',l,temp(l),teta(l) |
---|
90 | q(l,1) = qtprof(k)-frac*( qtprof(k)- qtprof(k-1)) |
---|
91 | u(l) = uprof(k)-frac*( uprof(k)- uprof(k-1)) |
---|
92 | v(l) = vprof(k)-frac*( vprof(k)- vprof(k-1)) |
---|
93 | ug(l) = ugprof(k)-frac*( ugprof(k)- ugprof(k-1)) |
---|
94 | vg(l) = vgprof(k)-frac*( vgprof(k)- vgprof(k-1)) |
---|
95 | IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2) & |
---|
96 | & -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2)) |
---|
97 | omega(l)= wfls(k)-frac*( wfls(k)- wfls(k-1)) |
---|
98 | dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1)) |
---|
99 | dt_cooling(l)=thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1)) |
---|
100 | elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1) |
---|
101 | ttt =tttprof(1) |
---|
102 | if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile |
---|
103 | temp(l) = ttt*(play(l)/pzero)**rkappa |
---|
104 | teta(l) = ttt |
---|
105 | else |
---|
106 | temp(l) = ttt |
---|
107 | teta(l) = ttt*(pzero/play(l))**rkappa |
---|
108 | endif |
---|
109 | q(l,1) = qtprof(1) |
---|
110 | u(l) = uprof(1) |
---|
111 | v(l) = vprof(1) |
---|
112 | ug(l) = ugprof(1) |
---|
113 | vg(l) = vgprof(1) |
---|
114 | omega(l)= wfls(1) |
---|
115 | IF (nq2>0) q(l,nq1:nq2)=qprof(1,nq1:nq2) |
---|
116 | dq_dyn(l,1) =dqtdtls(1) |
---|
117 | dt_cooling(l)=thlpcar(1) |
---|
118 | endif |
---|
119 | enddo |
---|
120 | |
---|
121 | temp(l)=max(min(temp(l),350.),150.) |
---|
122 | rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
123 | if (l .lt. llm) then |
---|
124 | zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l)) |
---|
125 | endif |
---|
126 | omega2(l)=-rho(l)*omega(l) |
---|
127 | omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s |
---|
128 | if (l>1) then |
---|
129 | if(zlay(l-1)>height(kmax)) then |
---|
130 | omega(l)=0.0 |
---|
131 | omega2(l)=0.0 |
---|
132 | endif |
---|
133 | endif |
---|
134 | if(q(l,1)<0.) q(l,1)=0.0 |
---|
135 | q(l,2) = 0.0 |
---|
136 | enddo |
---|
137 | |
---|
138 | endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire |
---|
139 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
140 | !--------------------------------------------------------------------- |
---|
141 | ! Forcing for GCSSold: |
---|
142 | !--------------------------------------------------------------------- |
---|
143 | if (forcing_GCSSold) then |
---|
144 | fich_gcssold_ctl = './forcing.ctl' |
---|
145 | fich_gcssold_dat = './forcing8.dat' |
---|
146 | call copie(llm,play,psurf,fich_gcssold_ctl) |
---|
147 | call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat, & |
---|
148 | & ht_gcssold,hq_gcssold,hw_gcssold, & |
---|
149 | & hu_gcssold,hv_gcssold, & |
---|
150 | & hthturb_gcssold,hqturb_gcssold,Ts_gcssold, & |
---|
151 | & imp_fcg_gcssold,ts_fcg_gcssold, & |
---|
152 | & Tp_fcg_gcssold,Turb_fcg_gcssold) |
---|
153 | print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold |
---|
154 | endif ! forcing_GCSSold |
---|
155 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
156 | !--------------------------------------------------------------------- |
---|
157 | ! Forcing for RICO: |
---|
158 | !--------------------------------------------------------------------- |
---|
159 | if (forcing_rico) then |
---|
160 | |
---|
161 | ! call writefield_phy('omega', omega,llm+1) |
---|
162 | fich_rico = 'rico.txt' |
---|
163 | call read_rico(fich_rico,nlev_rico,ps_rico,play & |
---|
164 | & ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico & |
---|
165 | & ,dth_rico,dqh_rico) |
---|
166 | print*, ' on a lu et prepare RICO' |
---|
167 | |
---|
168 | mxcalc=llm |
---|
169 | print *, airefi, ' airefi ' |
---|
170 | do l = 1, llm |
---|
171 | rho(l) = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l))) |
---|
172 | temp(l) = t_rico(l) |
---|
173 | q(l,1) = q_rico(l) |
---|
174 | q(l,2) = 0.0 |
---|
175 | u(l) = u_rico(l) |
---|
176 | v(l) = v_rico(l) |
---|
177 | ug(l)=u_rico(l) |
---|
178 | vg(l)=v_rico(l) |
---|
179 | omega(l) = -w_rico(l)*rg |
---|
180 | omega2(l) = omega(l)/rg*airefi |
---|
181 | enddo |
---|
182 | endif |
---|
183 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
184 | !--------------------------------------------------------------------- |
---|
185 | ! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) : |
---|
186 | !--------------------------------------------------------------------- |
---|
187 | |
---|
188 | if (forcing_toga) then |
---|
189 | |
---|
190 | ! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps): |
---|
191 | fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt' |
---|
192 | CALL read_togacoare(fich_toga,nlev_toga,nt_toga & |
---|
193 | & ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga & |
---|
194 | & ,ht_toga,vt_toga,hq_toga,vq_toga) |
---|
195 | |
---|
196 | write(*,*) 'Forcing TOGA lu' |
---|
197 | |
---|
198 | ! time interpolation for initial conditions: |
---|
199 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 |
---|
200 | CALL interp_toga_time(daytime,day1,annee_ref & |
---|
201 | & ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga & |
---|
202 | & ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga & |
---|
203 | & ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga & |
---|
204 | & ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof & |
---|
205 | & ,ht_prof,vt_prof,hq_prof,vq_prof) |
---|
206 | |
---|
207 | ! vertical interpolation: |
---|
208 | CALL interp_toga_vertical(play,nlev_toga,plev_prof & |
---|
209 | & ,t_prof,q_prof,u_prof,v_prof,w_prof & |
---|
210 | & ,ht_prof,vt_prof,hq_prof,vq_prof & |
---|
211 | & ,t_mod,q_mod,u_mod,v_mod,w_mod & |
---|
212 | & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) |
---|
213 | write(*,*) 'Profil initial forcing TOGA interpole' |
---|
214 | |
---|
215 | ! initial and boundary conditions : |
---|
216 | tsurf = ts_prof |
---|
217 | write(*,*) 'SST initiale: ',tsurf |
---|
218 | do l = 1, llm |
---|
219 | temp(l) = t_mod(l) |
---|
220 | q(l,1) = q_mod(l) |
---|
221 | q(l,2) = 0.0 |
---|
222 | u(l) = u_mod(l) |
---|
223 | v(l) = v_mod(l) |
---|
224 | omega(l) = w_mod(l) |
---|
225 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
226 | !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
227 | !? omega2(l)=-rho(l)*omega(l) |
---|
228 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
229 | d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l)) |
---|
230 | d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l)) |
---|
231 | d_q_adv(l,2) = 0.0 |
---|
232 | enddo |
---|
233 | |
---|
234 | endif ! forcing_toga |
---|
235 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
236 | !--------------------------------------------------------------------- |
---|
237 | ! Forcing from TWPICE experiment (Shaocheng et al. 2010) : |
---|
238 | !--------------------------------------------------------------------- |
---|
239 | |
---|
240 | if (forcing_twpice) then |
---|
241 | !read TWP-ICE forcings |
---|
242 | fich_twpice='d_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf' |
---|
243 | call read_twpice(fich_twpice,nlev_twpi,nt_twpi & |
---|
244 | & ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi & |
---|
245 | & ,ht_twpi,vt_twpi,hq_twpi,vq_twpi) |
---|
246 | |
---|
247 | write(*,*) 'Forcing TWP-ICE lu' |
---|
248 | !Time interpolation for initial conditions using TOGA interpolation routine |
---|
249 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 |
---|
250 | CALL interp_toga_time(daytime,day1,annee_ref & |
---|
251 | & ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi & |
---|
252 | & ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi & |
---|
253 | & ,ht_twpi,vt_twpi,hq_twpi,vq_twpi & |
---|
254 | & ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp & |
---|
255 | & ,u_proftwp,v_proftwp,w_proftwp & |
---|
256 | & ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp) |
---|
257 | |
---|
258 | ! vertical interpolation using TOGA interpolation routine: |
---|
259 | ! write(*,*)'avant interp vert', t_proftwp |
---|
260 | CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp & |
---|
261 | & ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp & |
---|
262 | & ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp & |
---|
263 | & ,t_mod,q_mod,u_mod,v_mod,w_mod & |
---|
264 | & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) |
---|
265 | ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod |
---|
266 | |
---|
267 | ! initial and boundary conditions : |
---|
268 | ! tsurf = ts_proftwp |
---|
269 | write(*,*) 'SST initiale: ',tsurf |
---|
270 | do l = 1, llm |
---|
271 | temp(l) = t_mod(l) |
---|
272 | q(l,1) = q_mod(l) |
---|
273 | q(l,2) = 0.0 |
---|
274 | u(l) = u_mod(l) |
---|
275 | v(l) = v_mod(l) |
---|
276 | omega(l) = w_mod(l) |
---|
277 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
278 | |
---|
279 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
280 | !on applique le forcage total au premier pas de temps |
---|
281 | !attention: signe different de toga |
---|
282 | d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l)) |
---|
283 | d_q_adv(l,1) = (hq_mod(l)+vq_mod(l)) |
---|
284 | d_q_adv(l,2) = 0.0 |
---|
285 | enddo |
---|
286 | |
---|
287 | endif !forcing_twpice |
---|
288 | |
---|
289 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
290 | !--------------------------------------------------------------------- |
---|
291 | ! Forcing from AMMA experiment (Couvreux et al. 2010) : |
---|
292 | !--------------------------------------------------------------------- |
---|
293 | |
---|
294 | if (forcing_amma) then |
---|
295 | |
---|
296 | call read_1D_cases |
---|
297 | |
---|
298 | write(*,*) 'Forcing AMMA lu' |
---|
299 | |
---|
300 | !champs initiaux: |
---|
301 | do k=1,nlev_amma |
---|
302 | th_ammai(k)=th_amma(k) |
---|
303 | q_ammai(k)=q_amma(k) |
---|
304 | u_ammai(k)=u_amma(k) |
---|
305 | v_ammai(k)=v_amma(k) |
---|
306 | vitw_ammai(k)=vitw_amma(k,12) |
---|
307 | ht_ammai(k)=ht_amma(k,12) |
---|
308 | hq_ammai(k)=hq_amma(k,12) |
---|
309 | vt_ammai(k)=0. |
---|
310 | vq_ammai(k)=0. |
---|
311 | enddo |
---|
312 | omega(:)=0. |
---|
313 | omega2(:)=0. |
---|
314 | rho(:)=0. |
---|
315 | ! vertical interpolation using TOGA interpolation routine: |
---|
316 | ! write(*,*)'avant interp vert', t_proftwp |
---|
317 | CALL interp_toga_vertical(play,nlev_amma,plev_amma & |
---|
318 | & ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai & |
---|
319 | & ,ht_ammai,vt_ammai,hq_ammai,vq_ammai & |
---|
320 | & ,t_mod,q_mod,u_mod,v_mod,w_mod & |
---|
321 | & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) |
---|
322 | ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod |
---|
323 | |
---|
324 | ! initial and boundary conditions : |
---|
325 | ! tsurf = ts_proftwp |
---|
326 | write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc |
---|
327 | do l = 1, llm |
---|
328 | ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp |
---|
329 | ! temp(l) = t_mod(l)*(play(l)/pzero)**rkappa |
---|
330 | temp(l) = t_mod(l) |
---|
331 | q(l,1) = q_mod(l) |
---|
332 | q(l,2) = 0.0 |
---|
333 | ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) |
---|
334 | u(l) = u_mod(l) |
---|
335 | v(l) = v_mod(l) |
---|
336 | rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
337 | omega(l) = w_mod(l)*(-rg*rho(l)) |
---|
338 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
339 | |
---|
340 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
341 | !on applique le forcage total au premier pas de temps |
---|
342 | !attention: signe different de toga |
---|
343 | d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) |
---|
344 | !forcage en th |
---|
345 | ! d_th_adv(l) = ht_mod(l) |
---|
346 | d_q_adv(l,1) = hq_mod(l) |
---|
347 | d_q_adv(l,2) = 0.0 |
---|
348 | dt_cooling(l)=0. |
---|
349 | enddo |
---|
350 | write(*,*) 'Prof initeforcing AMMA interpole temp39',temp(39) |
---|
351 | |
---|
352 | |
---|
353 | ! ok_flux_surf=.false. |
---|
354 | fsens=-1.*sens_amma(12) |
---|
355 | flat=-1.*lat_amma(12) |
---|
356 | |
---|
357 | endif !forcing_amma |
---|
358 | |
---|
359 | |
---|
360 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
361 | !--------------------------------------------------------------------- |
---|
362 | ! Forcing from DICE experiment (see file DICE_protocol_vn2-3.pdf) |
---|
363 | !--------------------------------------------------------------------- |
---|
364 | |
---|
365 | if (forcing_dice) then |
---|
366 | !read DICE forcings |
---|
367 | fich_dice='dice_driver.nc' |
---|
368 | call read_dice(fich_dice,nlev_dice,nt_dice & |
---|
369 | & ,zz_dice,plev_dice,t_dice,qv_dice,u_dice,v_dice,o3_dice & |
---|
370 | & ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice& |
---|
371 | & ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice & |
---|
372 | & ,hu_dice,hv_dice,w_dice,omega_dice) |
---|
373 | |
---|
374 | write(*,*) 'Forcing DICE lu' |
---|
375 | |
---|
376 | !champs initiaux: |
---|
377 | do k=1,nlev_dice |
---|
378 | t_dicei(k)=t_dice(k) |
---|
379 | qv_dicei(k)=qv_dice(k) |
---|
380 | u_dicei(k)=u_dice(k) |
---|
381 | v_dicei(k)=v_dice(k) |
---|
382 | o3_dicei(k)=o3_dice(k) |
---|
383 | ht_dicei(k)=ht_dice(k,1) |
---|
384 | hq_dicei(k)=hq_dice(k,1) |
---|
385 | hu_dicei(k)=hu_dice(k,1) |
---|
386 | hv_dicei(k)=hv_dice(k,1) |
---|
387 | w_dicei(k)=w_dice(k,1) |
---|
388 | omega_dicei(k)=omega_dice(k,1) |
---|
389 | enddo |
---|
390 | omega(:)=0. |
---|
391 | omega2(:)=0. |
---|
392 | rho(:)=0. |
---|
393 | ! vertical interpolation using TOGA interpolation routine: |
---|
394 | ! write(*,*)'avant interp vert', t_proftwp |
---|
395 | ! |
---|
396 | ! CALL interp_dice_time(daytime,day1,annee_ref |
---|
397 | ! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice |
---|
398 | ! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice |
---|
399 | ! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice |
---|
400 | ! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice |
---|
401 | ! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof |
---|
402 | ! o ,ustar_prof,psurf_prof,ug_profd,vg_profd |
---|
403 | ! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd |
---|
404 | ! o ,omega_profd) |
---|
405 | |
---|
406 | CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice & |
---|
407 | & ,t_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei & |
---|
408 | & ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei& |
---|
409 | & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & |
---|
410 | & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc) |
---|
411 | |
---|
412 | ! Pour tester les advections horizontales de T et Q, on met w_mod et omega_mod ?? zero (MPL 20131108) |
---|
413 | ! w_mod(:,:)=0. |
---|
414 | ! omega_mod(:,:)=0. |
---|
415 | |
---|
416 | ! write(*,*) 'Profil initial forcing DICE interpole',t_mod |
---|
417 | ! Les forcages DICE sont donnes /jour et non /seconde ! |
---|
418 | ht_mod(:)=ht_mod(:)/86400. |
---|
419 | hq_mod(:)=hq_mod(:)/86400. |
---|
420 | hu_mod(:)=hu_mod(:)/86400. |
---|
421 | hv_mod(:)=hv_mod(:)/86400. |
---|
422 | |
---|
423 | ! initial and boundary conditions : |
---|
424 | write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc |
---|
425 | do l = 1, llm |
---|
426 | ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp |
---|
427 | ! temp(l) = th_mod(l)*(play(l)/pzero)**rkappa |
---|
428 | temp(l) = t_mod(l) |
---|
429 | q(l,1) = qv_mod(l) |
---|
430 | q(l,2) = 0.0 |
---|
431 | ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) |
---|
432 | u(l) = u_mod(l) |
---|
433 | v(l) = v_mod(l) |
---|
434 | ug(l)=ug_dice(1) |
---|
435 | vg(l)=vg_dice(1) |
---|
436 | rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
437 | ! omega(l) = w_mod(l)*(-rg*rho(l)) |
---|
438 | omega(l) = omega_mod(l) |
---|
439 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
440 | |
---|
441 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
442 | !on applique le forcage total au premier pas de temps |
---|
443 | !attention: signe different de toga |
---|
444 | d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) |
---|
445 | !forcage en th |
---|
446 | ! d_th_adv(l) = ht_mod(l) |
---|
447 | d_q_adv(l,1) = hq_mod(l) |
---|
448 | d_q_adv(l,2) = 0.0 |
---|
449 | dt_cooling(l)=0. |
---|
450 | enddo |
---|
451 | write(*,*) 'Profil initial forcing DICE interpole temp39',temp(39) |
---|
452 | |
---|
453 | |
---|
454 | ! ok_flux_surf=.false. |
---|
455 | fsens=-1.*shf_dice(1) |
---|
456 | flat=-1.*lhf_dice(1) |
---|
457 | ! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par |
---|
458 | ! le coefficient de trainee en surface cd**2=ustar*vent(k=1) |
---|
459 | ! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface |
---|
460 | ! MPL 05082013 |
---|
461 | ust=ustar_dice(1) |
---|
462 | tg=tg_dice(1) |
---|
463 | print *,'ust= ',ust |
---|
464 | IF (tsurf .LE. 0.) THEN |
---|
465 | tsurf= tg_dice(1) |
---|
466 | ENDIF |
---|
467 | psurf= psurf_dice(1) |
---|
468 | solsw_in = (1.-albedo)/albedo*swup_dice(1) |
---|
469 | sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1) |
---|
470 | PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in |
---|
471 | endif !forcing_dice |
---|
472 | |
---|
473 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
474 | !--------------------------------------------------------------------- |
---|
475 | ! Forcing from GABLS4 experiment |
---|
476 | !--------------------------------------------------------------------- |
---|
477 | |
---|
478 | !!!! Si la temperature de surface n'est pas impos??e: |
---|
479 | |
---|
480 | if (forcing_gabls4) then |
---|
481 | !read GABLS4 forcings |
---|
482 | |
---|
483 | fich_gabls4='gabls4_driver.nc' |
---|
484 | |
---|
485 | |
---|
486 | call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 & |
---|
487 | & ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4) |
---|
488 | |
---|
489 | write(*,*) 'Forcing GABLS4 lu' |
---|
490 | |
---|
491 | !champs initiaux: |
---|
492 | do k=1,nlev_gabls4 |
---|
493 | t_gabi(k)=t_gabls4(k) |
---|
494 | qv_gabi(k)=qv_gabls4(k) |
---|
495 | u_gabi(k)=u_gabls4(k) |
---|
496 | v_gabi(k)=v_gabls4(k) |
---|
497 | poub(k)=0. |
---|
498 | ht_gabi(k)=ht_gabls4(k,1) |
---|
499 | hq_gabi(k)=hq_gabls4(k,1) |
---|
500 | ug_gabi(k)=ug_gabls4(k,1) |
---|
501 | vg_gabi(k)=vg_gabls4(k,1) |
---|
502 | enddo |
---|
503 | |
---|
504 | omega(:)=0. |
---|
505 | omega2(:)=0. |
---|
506 | rho(:)=0. |
---|
507 | ! vertical interpolation using TOGA interpolation routine: |
---|
508 | ! write(*,*)'avant interp vert', t_proftwp |
---|
509 | ! |
---|
510 | ! CALL interp_dice_time(daytime,day1,annee_ref |
---|
511 | ! i ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice |
---|
512 | ! i ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice |
---|
513 | ! i ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice |
---|
514 | ! i ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice |
---|
515 | ! o ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof |
---|
516 | ! o ,ustar_prof,psurf_prof,ug_profd,vg_profd |
---|
517 | ! o ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd |
---|
518 | ! o ,omega_profd) |
---|
519 | |
---|
520 | CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4 & |
---|
521 | & ,t_gabi,qv_gabi,u_gabi,v_gabi,poub & |
---|
522 | & ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,poub,poub & |
---|
523 | & ,t_mod,qv_mod,u_mod,v_mod,o3_mod & |
---|
524 | & ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc) |
---|
525 | |
---|
526 | ! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024 |
---|
527 | ! ht_mod(:)=ht_mod(:)/86400. |
---|
528 | ! hq_mod(:)=hq_mod(:)/86400. |
---|
529 | |
---|
530 | ! initial and boundary conditions : |
---|
531 | write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc |
---|
532 | do l = 1, llm |
---|
533 | ! Ligne du dessous ?? decommenter si on lit theta au lieu de temp |
---|
534 | ! temp(l) = th_mod(l)*(play(l)/pzero)**rkappa |
---|
535 | temp(l) = t_mod(l) |
---|
536 | q(l,1) = qv_mod(l) |
---|
537 | q(l,2) = 0.0 |
---|
538 | ! print *,'read_forc: l,temp,q=',l,temp(l),q(l,1) |
---|
539 | u(l) = u_mod(l) |
---|
540 | v(l) = v_mod(l) |
---|
541 | ug(l)=ug_mod(l) |
---|
542 | vg(l)=vg_mod(l) |
---|
543 | |
---|
544 | ! |
---|
545 | ! tg=tsurf |
---|
546 | ! |
---|
547 | |
---|
548 | print *,'***** tsurf=',tsurf |
---|
549 | rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
550 | ! omega(l) = w_mod(l)*(-rg*rho(l)) |
---|
551 | omega(l) = omega_mod(l) |
---|
552 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
553 | |
---|
554 | |
---|
555 | |
---|
556 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
557 | !on applique le forcage total au premier pas de temps |
---|
558 | !attention: signe different de toga |
---|
559 | ! d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) |
---|
560 | !forcage en th |
---|
561 | d_th_adv(l) = ht_mod(l) |
---|
562 | d_q_adv(l,1) = hq_mod(l) |
---|
563 | d_q_adv(l,2) = 0.0 |
---|
564 | dt_cooling(l)=0. |
---|
565 | enddo |
---|
566 | |
---|
567 | !--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024--------------- |
---|
568 | ! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par |
---|
569 | ! le coefficient de trainee en surface cd**2=ustar*vent(k=1) |
---|
570 | ! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface |
---|
571 | ! MPL 05082013 |
---|
572 | ! ust=ustar_dice(1) |
---|
573 | ! tg=tg_dice(1) |
---|
574 | ! print *,'ust= ',ust |
---|
575 | ! IF (tsurf .LE. 0.) THEN |
---|
576 | ! tsurf= tg_dice(1) |
---|
577 | ! ENDIF |
---|
578 | ! psurf= psurf_dice(1) |
---|
579 | ! solsw_in = (1.-albedo)/albedo*swup_dice(1) |
---|
580 | ! sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1) |
---|
581 | ! PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in |
---|
582 | !-------------------------------------------------------------------------------------- |
---|
583 | endif !forcing_gabls4 |
---|
584 | |
---|
585 | |
---|
586 | |
---|
587 | ! Forcing from Arm_Cu case |
---|
588 | ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes |
---|
589 | ! large scale advective forcing,radiative forcing |
---|
590 | ! and advective tendency of theta and qt to be applied |
---|
591 | !--------------------------------------------------------------------- |
---|
592 | |
---|
593 | if (forcing_armcu) then |
---|
594 | ! read armcu forcing : |
---|
595 | write(*,*) 'Avant lecture Forcing Arm_Cu' |
---|
596 | fich_armcu = './ifa_armcu.txt' |
---|
597 | CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu, & |
---|
598 | & sens_armcu,flat_armcu,adv_theta_armcu, & |
---|
599 | & rad_theta_armcu,adv_qt_armcu) |
---|
600 | write(*,*) 'Forcing Arm_Cu lu' |
---|
601 | |
---|
602 | !---------------------------------------------------------------------- |
---|
603 | ! Read profiles from file: prof.inp.19 or prof.inp.40 |
---|
604 | ! For this case, profiles are given for two vertical resolution |
---|
605 | ! 19 or 40 levels |
---|
606 | ! |
---|
607 | ! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html |
---|
608 | ! Note that the initial profiles contain no liquid water! |
---|
609 | ! (so potential temperature can be interpreted as liquid water |
---|
610 | ! potential temperature and water vapor as total water) |
---|
611 | ! profiles are given at full levels |
---|
612 | !---------------------------------------------------------------------- |
---|
613 | |
---|
614 | call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod, & |
---|
615 | & v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp) |
---|
616 | |
---|
617 | ! time interpolation for initial conditions: |
---|
618 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 |
---|
619 | |
---|
620 | print *,'Avant interp_armcu_time' |
---|
621 | print *,'daytime=',daytime |
---|
622 | print *,'day1=',day1 |
---|
623 | print *,'annee_ref=',annee_ref |
---|
624 | print *,'year_ini_armcu=',year_ini_armcu |
---|
625 | print *,'day_ju_ini_armcu=',day_ju_ini_armcu |
---|
626 | print *,'nt_armcu=',nt_armcu |
---|
627 | print *,'dt_armcu=',dt_armcu |
---|
628 | print *,'nlev_armcu=',nlev_armcu |
---|
629 | CALL interp_armcu_time(daytime,day1,annee_ref & |
---|
630 | & ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu & |
---|
631 | & ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu & |
---|
632 | & ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof & |
---|
633 | & ,adv_theta_prof,rad_theta_prof,adv_qt_prof) |
---|
634 | write(*,*) 'Forcages interpoles dans temps' |
---|
635 | |
---|
636 | ! No vertical interpolation if nlev imposed to 19 or 40 |
---|
637 | ! The vertical grid stops at 4000m # 600hPa |
---|
638 | mxcalc=llm |
---|
639 | |
---|
640 | ! initial and boundary conditions : |
---|
641 | ! tsurf = ts_prof |
---|
642 | ! tsurf read in lmdz1d.def |
---|
643 | write(*,*) 'Tsurf initiale: ',tsurf |
---|
644 | do l = 1, llm |
---|
645 | play(l)=play_mod(l)*100. |
---|
646 | presnivs(l)=play(l) |
---|
647 | zlay(l)=height(l) |
---|
648 | temp(l) = t_mod(l) |
---|
649 | teta(l)=theta_mod(l) |
---|
650 | q(l,1) = qv_mod(l)/1000. |
---|
651 | ! No liquid water in the initial profil |
---|
652 | q(l,2) = 0. |
---|
653 | u(l) = u_mod(l) |
---|
654 | ug(l)= u_mod(l) |
---|
655 | v(l) = v_mod(l) |
---|
656 | vg(l)= v_mod(l) |
---|
657 | ! Advective forcings are given in K or g/kg ... per HOUR |
---|
658 | ! IF(height(l).LT.1000) THEN |
---|
659 | ! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600. |
---|
660 | ! d_q_adv(l,1) = adv_qt_prof/1000./3600. |
---|
661 | ! d_q_adv(l,2) = 0.0 |
---|
662 | ! ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN |
---|
663 | ! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)* |
---|
664 | ! : (1-(height(l)-1000.)/2000.) |
---|
665 | ! d_th_adv(l) = d_th_adv(l)/3600. |
---|
666 | ! d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.) |
---|
667 | ! d_q_adv(l,1) = d_q_adv(l,1)/1000./3600. |
---|
668 | ! d_q_adv(l,2) = 0.0 |
---|
669 | ! ELSE |
---|
670 | ! d_th_adv(l) = 0.0 |
---|
671 | ! d_q_adv(l,1) = 0.0 |
---|
672 | ! d_q_adv(l,2) = 0.0 |
---|
673 | ! ENDIF |
---|
674 | enddo |
---|
675 | ! plev at half levels is given in proh.inp.19 or proh.inp.40 files |
---|
676 | plev(1)= ap(llm+1)+bp(llm+1)*psurf |
---|
677 | do l = 1, llm |
---|
678 | plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf |
---|
679 | print *,'Read_forc: l height play plev zlay temp', & |
---|
680 | & l,height(l),play(l),plev(l),zlay(l),temp(l) |
---|
681 | enddo |
---|
682 | ! For this case, fluxes are imposed |
---|
683 | fsens=-1*sens_prof |
---|
684 | flat=-1*flat_prof |
---|
685 | |
---|
686 | endif ! forcing_armcu |
---|
687 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
688 | !--------------------------------------------------------------------- |
---|
689 | ! Forcing from transition case of Irina Sandu |
---|
690 | !--------------------------------------------------------------------- |
---|
691 | |
---|
692 | if (forcing_sandu) then |
---|
693 | write(*,*) 'Avant lecture Forcing SANDU' |
---|
694 | |
---|
695 | ! read sanduref forcing : |
---|
696 | fich_sandu = './ifa_sanduref.txt' |
---|
697 | CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu) |
---|
698 | |
---|
699 | write(*,*) 'Forcing SANDU lu' |
---|
700 | |
---|
701 | !---------------------------------------------------------------------- |
---|
702 | ! Read profiles from file: prof.inp.001 |
---|
703 | !---------------------------------------------------------------------- |
---|
704 | |
---|
705 | call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs, & |
---|
706 | & thl_profs,q_profs,u_profs,v_profs, & |
---|
707 | & w_profs,omega_profs,o3mmr_profs) |
---|
708 | |
---|
709 | ! time interpolation for initial conditions: |
---|
710 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 |
---|
711 | ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! |
---|
712 | ! revoir 1DUTILS.h et les arguments |
---|
713 | |
---|
714 | print *,'Avant interp_sandu_time' |
---|
715 | print *,'daytime=',daytime |
---|
716 | print *,'day1=',day1 |
---|
717 | print *,'annee_ref=',annee_ref |
---|
718 | print *,'year_ini_sandu=',year_ini_sandu |
---|
719 | print *,'day_ju_ini_sandu=',day_ju_ini_sandu |
---|
720 | print *,'nt_sandu=',nt_sandu |
---|
721 | print *,'dt_sandu=',dt_sandu |
---|
722 | print *,'nlev_sandu=',nlev_sandu |
---|
723 | CALL interp_sandu_time(daytime,day1,annee_ref & |
---|
724 | & ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu & |
---|
725 | & ,nlev_sandu & |
---|
726 | & ,ts_sandu,ts_prof) |
---|
727 | |
---|
728 | ! vertical interpolation: |
---|
729 | print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu |
---|
730 | CALL interp_sandu_vertical(play,nlev_sandu,plev_profs & |
---|
731 | & ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs & |
---|
732 | & ,omega_profs,o3mmr_profs & |
---|
733 | & ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod & |
---|
734 | & ,omega_mod,o3mmr_mod,mxcalc) |
---|
735 | write(*,*) 'Profil initial forcing SANDU interpole' |
---|
736 | |
---|
737 | ! initial and boundary conditions : |
---|
738 | tsurf = ts_prof |
---|
739 | write(*,*) 'SST initiale: ',tsurf |
---|
740 | do l = 1, llm |
---|
741 | temp(l) = t_mod(l) |
---|
742 | tetal(l)=thl_mod(l) |
---|
743 | q(l,1) = q_mod(l) |
---|
744 | q(l,2) = 0.0 |
---|
745 | u(l) = u_mod(l) |
---|
746 | v(l) = v_mod(l) |
---|
747 | w(l) = w_mod(l) |
---|
748 | omega(l) = omega_mod(l) |
---|
749 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
750 | !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
751 | !? omega2(l)=-rho(l)*omega(l) |
---|
752 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
753 | ! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) |
---|
754 | ! d_q_adv(l,1) = vq_mod(l) |
---|
755 | d_th_adv(l) = alpha*omega(l)/rcpd |
---|
756 | d_q_adv(l,1) = 0.0 |
---|
757 | d_q_adv(l,2) = 0.0 |
---|
758 | enddo |
---|
759 | |
---|
760 | endif ! forcing_sandu |
---|
761 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
762 | !--------------------------------------------------------------------- |
---|
763 | ! Forcing from Astex case |
---|
764 | !--------------------------------------------------------------------- |
---|
765 | |
---|
766 | if (forcing_astex) then |
---|
767 | write(*,*) 'Avant lecture Forcing Astex' |
---|
768 | |
---|
769 | ! read astex forcing : |
---|
770 | fich_astex = './ifa_astex.txt' |
---|
771 | CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex, & |
---|
772 | & ug_astex,vg_astex,ufa_astex,vfa_astex) |
---|
773 | |
---|
774 | write(*,*) 'Forcing Astex lu' |
---|
775 | |
---|
776 | !---------------------------------------------------------------------- |
---|
777 | ! Read profiles from file: prof.inp.001 |
---|
778 | !---------------------------------------------------------------------- |
---|
779 | |
---|
780 | call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa, & |
---|
781 | & thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa, & |
---|
782 | & w_profa,tke_profa,o3mmr_profa) |
---|
783 | |
---|
784 | ! time interpolation for initial conditions: |
---|
785 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 |
---|
786 | ! ATTENTION, cet appel ne convient pas pour le cas SANDU !! |
---|
787 | ! revoir 1DUTILS.h et les arguments |
---|
788 | |
---|
789 | print *,'Avant interp_astex_time' |
---|
790 | print *,'daytime=',daytime |
---|
791 | print *,'day1=',day1 |
---|
792 | print *,'annee_ref=',annee_ref |
---|
793 | print *,'year_ini_astex=',year_ini_astex |
---|
794 | print *,'day_ju_ini_astex=',day_ju_ini_astex |
---|
795 | print *,'nt_astex=',nt_astex |
---|
796 | print *,'dt_astex=',dt_astex |
---|
797 | print *,'nlev_astex=',nlev_astex |
---|
798 | CALL interp_astex_time(daytime,day1,annee_ref & |
---|
799 | & ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex & |
---|
800 | & ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex & |
---|
801 | & ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof & |
---|
802 | & ,ufa_prof,vfa_prof) |
---|
803 | |
---|
804 | ! vertical interpolation: |
---|
805 | print *,'Avant interp_vertical: nlev_astex=',nlev_astex |
---|
806 | CALL interp_astex_vertical(play,nlev_astex,plev_profa & |
---|
807 | & ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa & |
---|
808 | & ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa & |
---|
809 | & ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod & |
---|
810 | & ,tke_mod,o3mmr_mod,mxcalc) |
---|
811 | write(*,*) 'Profil initial forcing Astex interpole' |
---|
812 | |
---|
813 | ! initial and boundary conditions : |
---|
814 | tsurf = ts_prof |
---|
815 | write(*,*) 'SST initiale: ',tsurf |
---|
816 | do l = 1, llm |
---|
817 | temp(l) = t_mod(l) |
---|
818 | tetal(l)=thl_mod(l) |
---|
819 | q(l,1) = qv_mod(l) |
---|
820 | q(l,2) = ql_mod(l) |
---|
821 | u(l) = u_mod(l) |
---|
822 | v(l) = v_mod(l) |
---|
823 | w(l) = w_mod(l) |
---|
824 | omega(l) = w_mod(l) |
---|
825 | ! omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
826 | ! rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
827 | ! omega2(l)=-rho(l)*omega(l) |
---|
828 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
829 | ! d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) |
---|
830 | ! d_q_adv(l,1) = vq_mod(l) |
---|
831 | d_th_adv(l) = alpha*omega(l)/rcpd |
---|
832 | d_q_adv(l,1) = 0.0 |
---|
833 | d_q_adv(l,2) = 0.0 |
---|
834 | enddo |
---|
835 | |
---|
836 | endif ! forcing_astex |
---|
837 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
838 | !--------------------------------------------------------------------- |
---|
839 | ! Forcing from standard case : |
---|
840 | !--------------------------------------------------------------------- |
---|
841 | |
---|
842 | if (forcing_case) then |
---|
843 | |
---|
844 | write(*,*),'avant call read_1D_cas' |
---|
845 | call read_1D_cas |
---|
846 | write(*,*) 'Forcing read' |
---|
847 | |
---|
848 | !Time interpolation for initial conditions using TOGA interpolation routine |
---|
849 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 |
---|
850 | CALL interp_case_time(day,day1,annee_ref & |
---|
851 | ! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas & |
---|
852 | & ,nt_cas,nlev_cas & |
---|
853 | & ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas & |
---|
854 | & ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas & |
---|
855 | & ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas & |
---|
856 | & ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas & |
---|
857 | & ,uw_cas,vw_cas,q1_cas,q2_cas & |
---|
858 | & ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas & |
---|
859 | & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas & |
---|
860 | & ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas & |
---|
861 | & ,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas,ustar_prof_cas & |
---|
862 | & ,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas) |
---|
863 | |
---|
864 | ! vertical interpolation using TOGA interpolation routine: |
---|
865 | ! write(*,*)'avant interp vert', t_prof |
---|
866 | CALL interp_case_vertical(play,nlev_cas,plev_prof_cas & |
---|
867 | & ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas & |
---|
868 | & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & |
---|
869 | & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & |
---|
870 | & ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas & |
---|
871 | & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & |
---|
872 | & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc) |
---|
873 | ! write(*,*) 'Profil initial forcing case interpole',t_mod |
---|
874 | |
---|
875 | ! initial and boundary conditions : |
---|
876 | ! tsurf = ts_prof_cas |
---|
877 | ts_cur = ts_prof_cas |
---|
878 | psurf=plev_prof_cas(1) |
---|
879 | write(*,*) 'SST initiale: ',tsurf |
---|
880 | do l = 1, llm |
---|
881 | temp(l) = t_mod_cas(l) |
---|
882 | q(l,1) = q_mod_cas(l) |
---|
883 | q(l,2) = 0.0 |
---|
884 | u(l) = u_mod_cas(l) |
---|
885 | v(l) = v_mod_cas(l) |
---|
886 | omega(l) = w_mod_cas(l) |
---|
887 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
888 | |
---|
889 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
890 | !on applique le forcage total au premier pas de temps |
---|
891 | !attention: signe different de toga |
---|
892 | d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) |
---|
893 | d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l)) |
---|
894 | d_q_adv(l,2) = 0.0 |
---|
895 | d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l)) |
---|
896 | d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l)) |
---|
897 | enddo |
---|
898 | |
---|
899 | ! In case fluxes are imposed |
---|
900 | IF (ok_flux_surf) THEN |
---|
901 | fsens=sens_prof_cas |
---|
902 | flat=lat_prof_cas |
---|
903 | ENDIF |
---|
904 | IF (ok_prescr_ust) THEN |
---|
905 | ust=ustar_prof_cas |
---|
906 | print *,'ust=',ust |
---|
907 | ENDIF |
---|
908 | |
---|
909 | endif !forcing_case |
---|
910 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
911 | !--------------------------------------------------------------------- |
---|
912 | ! Forcing from standard case : |
---|
913 | !--------------------------------------------------------------------- |
---|
914 | |
---|
915 | if (forcing_case2) then |
---|
916 | |
---|
917 | write(*,*),'avant call read2_1D_cas' |
---|
918 | call read2_1D_cas |
---|
919 | write(*,*) 'Forcing read' |
---|
920 | |
---|
921 | !Time interpolation for initial conditions using interpolation routine |
---|
922 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 |
---|
923 | CALL interp2_case_time(daytime,day1,annee_ref & |
---|
924 | ! & ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas & |
---|
925 | & ,nt_cas,nlev_cas & |
---|
926 | & ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas & |
---|
927 | & ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas & |
---|
928 | & ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas & |
---|
929 | & ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas & |
---|
930 | & ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas & |
---|
931 | ! |
---|
932 | & ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas & |
---|
933 | & ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas & |
---|
934 | & ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & |
---|
935 | & ,du_prof_cas,hu_prof_cas,vu_prof_cas & |
---|
936 | & ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas & |
---|
937 | & ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & |
---|
938 | & ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas & |
---|
939 | & ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas) |
---|
940 | |
---|
941 | do l = 1, nlev_cas |
---|
942 | print *,'apres 1ere interp: plev_cas, plev_prof_cas=',l,plev_cas(l,1),plev_prof_cas(l) |
---|
943 | enddo |
---|
944 | |
---|
945 | ! vertical interpolation using interpolation routine: |
---|
946 | ! write(*,*)'avant interp vert', t_prof |
---|
947 | CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas & |
---|
948 | & ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas & |
---|
949 | & ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas & |
---|
950 | & ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas & |
---|
951 | & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas & |
---|
952 | & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas & |
---|
953 | & ,dth_prof_cas,hth_prof_cas,vth_prof_cas & |
---|
954 | ! |
---|
955 | & ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas & |
---|
956 | & ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas & |
---|
957 | & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas & |
---|
958 | & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas & |
---|
959 | & ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc) |
---|
960 | |
---|
961 | ! write(*,*) 'Profil initial forcing case interpole',t_mod |
---|
962 | |
---|
963 | ! initial and boundary conditions : |
---|
964 | ! tsurf = ts_prof_cas |
---|
965 | ts_cur = ts_prof_cas |
---|
966 | psurf=plev_prof_cas(1) |
---|
967 | write(*,*) 'SST initiale: ',tsurf |
---|
968 | do l = 1, llm |
---|
969 | temp(l) = t_mod_cas(l) |
---|
970 | q(l,1) = qv_mod_cas(l) |
---|
971 | q(l,2) = ql_mod_cas(l) |
---|
972 | u(l) = u_mod_cas(l) |
---|
973 | ug(l)= u_mod_cas(l) |
---|
974 | v(l) = v_mod_cas(l) |
---|
975 | vg(l)= v_mod_cas(l) |
---|
976 | omega(l) = w_mod_cas(l) |
---|
977 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
978 | |
---|
979 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
980 | !on applique le forcage total au premier pas de temps |
---|
981 | !attention: signe different de toga |
---|
982 | d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) |
---|
983 | d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) |
---|
984 | ! d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l)) |
---|
985 | d_q_adv(l,1) = dq_mod_cas(l) |
---|
986 | d_q_adv(l,2) = 0.0 |
---|
987 | ! d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l)) |
---|
988 | d_u_adv(l) = du_mod_cas(l) |
---|
989 | ! d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l)) |
---|
990 | d_u_adv(l) = dv_mod_cas(l) |
---|
991 | enddo |
---|
992 | |
---|
993 | ! Faut-il multiplier par -1 ? (MPL 20160713) |
---|
994 | IF (ok_flux_surf) THEN |
---|
995 | fsens=-1.*sens_prof_cas |
---|
996 | flat=-1.*lat_prof_cas |
---|
997 | ENDIF |
---|
998 | ! |
---|
999 | IF (ok_prescr_ust) THEN |
---|
1000 | ust=ustar_prof_cas |
---|
1001 | print *,'ust=',ust |
---|
1002 | ENDIF |
---|
1003 | |
---|
1004 | endif !forcing_case2 |
---|
1005 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1006 | |
---|