1 | !---------------------------------------------------------------------- |
---|
2 | ! forcing_les = .T. : Impose a constant cooling |
---|
3 | ! forcing_radconv = .T. : Pure radiative-convective equilibrium: |
---|
4 | !---------------------------------------------------------------------- |
---|
5 | |
---|
6 | if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then |
---|
7 | |
---|
8 | !---------------------------------------------------------------------- |
---|
9 | ! Read profiles from files: prof.inp.001 and lscale.inp.001 |
---|
10 | ! (repris de readlesfiles) |
---|
11 | !---------------------------------------------------------------------- |
---|
12 | |
---|
13 | call readprofiles(nlev_max,kmax,height, |
---|
14 | . tttprof,qtprof,uprof,vprof, |
---|
15 | . e12prof,ugprof,vgprof, |
---|
16 | . wfls,dqtdxls,dqtdyls,dqtdtls, |
---|
17 | . thlpcar) |
---|
18 | |
---|
19 | ! compute altitudes of play levels. |
---|
20 | zlay(1) =zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) |
---|
21 | do l = 2,llm |
---|
22 | zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf) |
---|
23 | enddo |
---|
24 | |
---|
25 | !---------------------------------------------------------------------- |
---|
26 | ! Interpolation of the profiles given on the input file to |
---|
27 | ! model levels |
---|
28 | !---------------------------------------------------------------------- |
---|
29 | zlay(1) = zsurf + rd*tsurf*(psurf-play(1))/(rg*psurf) |
---|
30 | do l=1,llm |
---|
31 | ! Above the max altutide of the input file |
---|
32 | |
---|
33 | if (zlay(l)<height(kmax)) mxcalc=l |
---|
34 | |
---|
35 | frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1)) |
---|
36 | ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1)) |
---|
37 | if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile |
---|
38 | temp(l) = ttt*(play(l)/pzero)**rkappa |
---|
39 | teta(l) = ttt |
---|
40 | else |
---|
41 | temp(l) = ttt |
---|
42 | teta(l) = ttt*(pzero/play(l))**rkappa |
---|
43 | endif |
---|
44 | print *,' temp,teta ',l,temp(l),teta(l) |
---|
45 | q(l,1) = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1)) |
---|
46 | u(l) = uprof(kmax)-frac*( uprof(kmax)- uprof(kmax-1)) |
---|
47 | v(l) = vprof(kmax)-frac*( vprof(kmax)- vprof(kmax-1)) |
---|
48 | ug(l) = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1)) |
---|
49 | vg(l) = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1)) |
---|
50 | omega(l)= wfls(kmax)-frac*( wfls(kmax)- wfls(kmax-1)) |
---|
51 | |
---|
52 | dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1)) |
---|
53 | dt_cooling(l) |
---|
54 | . =thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1)) |
---|
55 | do k=2,kmax |
---|
56 | frac = (height(k)-zlay(l))/(height(k)-height(k-1)) |
---|
57 | if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k) |
---|
58 | if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then |
---|
59 | ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1)) |
---|
60 | if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile |
---|
61 | temp(l) = ttt*(play(l)/pzero)**rkappa |
---|
62 | teta(l) = ttt |
---|
63 | else |
---|
64 | temp(l) = ttt |
---|
65 | teta(l) = ttt*(pzero/play(l))**rkappa |
---|
66 | endif |
---|
67 | print *,' temp,teta ',l,temp(l),teta(l) |
---|
68 | q(l,1) = qtprof(k)-frac*( qtprof(k)- qtprof(k-1)) |
---|
69 | u(l) = uprof(k)-frac*( uprof(k)- uprof(k-1)) |
---|
70 | v(l) = vprof(k)-frac*( vprof(k)- vprof(k-1)) |
---|
71 | ug(l) = ugprof(k)-frac*( ugprof(k)- ugprof(k-1)) |
---|
72 | vg(l) = vgprof(k)-frac*( vgprof(k)- vgprof(k-1)) |
---|
73 | omega(l)= wfls(k)-frac*( wfls(k)- wfls(k-1)) |
---|
74 | dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1)) |
---|
75 | dt_cooling(l) |
---|
76 | . =thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1)) |
---|
77 | elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1) |
---|
78 | ttt =tttprof(1) |
---|
79 | if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile |
---|
80 | temp(l) = ttt*(play(l)/pzero)**rkappa |
---|
81 | teta(l) = ttt |
---|
82 | else |
---|
83 | temp(l) = ttt |
---|
84 | teta(l) = ttt*(pzero/play(l))**rkappa |
---|
85 | endif |
---|
86 | q(l,1) = qtprof(1) |
---|
87 | u(l) = uprof(1) |
---|
88 | v(l) = vprof(1) |
---|
89 | ug(l) = ugprof(1) |
---|
90 | vg(l) = vgprof(1) |
---|
91 | omega(l)= wfls(1) |
---|
92 | dq_dyn(l,1) =dqtdtls(1) |
---|
93 | dt_cooling(l)=thlpcar(1) |
---|
94 | endif |
---|
95 | enddo |
---|
96 | |
---|
97 | temp(l)=max(min(temp(l),350.),150.) |
---|
98 | rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
99 | if (l .lt. llm) then |
---|
100 | zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l)) |
---|
101 | endif |
---|
102 | omega2(l)=-rho(l)*omega(l) |
---|
103 | omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s |
---|
104 | if (l>1) then |
---|
105 | if(zlay(l-1)>height(kmax)) then |
---|
106 | omega(l)=0.0 |
---|
107 | omega2(l)=0.0 |
---|
108 | endif |
---|
109 | endif |
---|
110 | if(q(l,1)<0.) q(l,1)=0.0 |
---|
111 | q(l,2) = 0.0 |
---|
112 | enddo |
---|
113 | |
---|
114 | endif ! forcing_les .or. forcing_GCSSold |
---|
115 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
116 | !--------------------------------------------------------------------- |
---|
117 | ! Forcing for GCSSold: |
---|
118 | !--------------------------------------------------------------------- |
---|
119 | if (forcing_GCSSold) then |
---|
120 | fich_gcssold_ctl = './forcing.ctl' |
---|
121 | fich_gcssold_dat = './forcing8.dat' |
---|
122 | call copie(llm,play,psurf,fich_gcssold_ctl) |
---|
123 | call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat, |
---|
124 | : ht_gcssold,hq_gcssold,hw_gcssold, |
---|
125 | : hu_gcssold,hv_gcssold, |
---|
126 | : hthturb_gcssold,hqturb_gcssold,Ts_gcssold, |
---|
127 | : imp_fcg_gcssold,ts_fcg_gcssold, |
---|
128 | : Tp_fcg_gcssold,Turb_fcg_gcssold) |
---|
129 | print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold |
---|
130 | endif ! forcing_GCSSold |
---|
131 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
132 | !--------------------------------------------------------------------- |
---|
133 | ! Forcing for RICO: |
---|
134 | !--------------------------------------------------------------------- |
---|
135 | if (forcing_rico) then |
---|
136 | |
---|
137 | ! call writefield_phy('omega', omega,llm+1) |
---|
138 | fich_rico = 'rico.txt' |
---|
139 | call read_rico(fich_rico,nlev_rico,ps_rico,play |
---|
140 | : ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico |
---|
141 | : ,dth_rico,dqh_rico) |
---|
142 | print*, ' on a lu et prepare RICO' |
---|
143 | |
---|
144 | mxcalc=llm |
---|
145 | print *, airefi, ' airefi ' |
---|
146 | do l = 1, llm |
---|
147 | rho(l) = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l))) |
---|
148 | temp(l) = t_rico(l) |
---|
149 | q(l,1) = q_rico(l) |
---|
150 | q(l,2) = 0.0 |
---|
151 | u(l) = u_rico(l) |
---|
152 | v(l) = v_rico(l) |
---|
153 | ug(l)=u_rico(l) |
---|
154 | vg(l)=v_rico(l) |
---|
155 | omega(l) = -w_rico(l)*rg |
---|
156 | omega2(l) = omega(l)/rg*airefi |
---|
157 | enddo |
---|
158 | endif |
---|
159 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
160 | !--------------------------------------------------------------------- |
---|
161 | ! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) : |
---|
162 | !--------------------------------------------------------------------- |
---|
163 | |
---|
164 | if (forcing_toga) then |
---|
165 | |
---|
166 | ! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps): |
---|
167 | fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt' |
---|
168 | CALL read_togacoare(fich_toga,nlev_toga,nt_toga |
---|
169 | : ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga |
---|
170 | : ,ht_toga,vt_toga,hq_toga,vq_toga) |
---|
171 | |
---|
172 | write(*,*) 'Forcing TOGA lu' |
---|
173 | |
---|
174 | ! time interpolation for initial conditions: |
---|
175 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 |
---|
176 | CALL interp_toga_time(daytime,day1,annee_ref |
---|
177 | i ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga |
---|
178 | i ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga |
---|
179 | i ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga |
---|
180 | o ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof |
---|
181 | o ,ht_prof,vt_prof,hq_prof,vq_prof) |
---|
182 | |
---|
183 | ! vertical interpolation: |
---|
184 | CALL interp_toga_vertical(play,nlev_toga,plev_prof |
---|
185 | : ,t_prof,q_prof,u_prof,v_prof,w_prof |
---|
186 | : ,ht_prof,vt_prof,hq_prof,vq_prof |
---|
187 | : ,t_mod,q_mod,u_mod,v_mod,w_mod |
---|
188 | : ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) |
---|
189 | write(*,*) 'Profil initial forcing TOGA interpole' |
---|
190 | |
---|
191 | ! initial and boundary conditions : |
---|
192 | tsurf = ts_prof |
---|
193 | write(*,*) 'SST initiale: ',tsurf |
---|
194 | do l = 1, llm |
---|
195 | temp(l) = t_mod(l) |
---|
196 | q(l,1) = q_mod(l) |
---|
197 | q(l,2) = 0.0 |
---|
198 | u(l) = u_mod(l) |
---|
199 | v(l) = v_mod(l) |
---|
200 | omega(l) = w_mod(l) |
---|
201 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
202 | !? rho(l) = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))) |
---|
203 | !? omega2(l)=-rho(l)*omega(l) |
---|
204 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
205 | d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l)) |
---|
206 | d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l)) |
---|
207 | d_q_adv(l,2) = 0.0 |
---|
208 | enddo |
---|
209 | |
---|
210 | endif ! forcing_toga |
---|
211 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
212 | !--------------------------------------------------------------------- |
---|
213 | ! Forcing from TWPICE experiment (Shaocheng et al. 2010) : |
---|
214 | !--------------------------------------------------------------------- |
---|
215 | |
---|
216 | if (forcing_twpice) then |
---|
217 | !read TWP-ICE forcings |
---|
218 | fich_twpice= |
---|
219 | : 'd_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf' |
---|
220 | call read_twpice(fich_twpice,nlev_twpi,nt_twpi |
---|
221 | : ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi |
---|
222 | : ,ht_twpi,vt_twpi,hq_twpi,vq_twpi) |
---|
223 | |
---|
224 | write(*,*) 'Forcing TWP-ICE lu' |
---|
225 | !Time interpolation for initial conditions using TOGA interpolation routine |
---|
226 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1 |
---|
227 | CALL interp_toga_time(daytime,day1,annee_ref |
---|
228 | i ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi |
---|
229 | i ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi |
---|
230 | i ,ht_twpi,vt_twpi,hq_twpi,vq_twpi |
---|
231 | o ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp |
---|
232 | o ,u_proftwp,v_proftwp,w_proftwp |
---|
233 | o ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp) |
---|
234 | |
---|
235 | ! vertical interpolation using TOGA interpolation routine: |
---|
236 | ! write(*,*)'avant interp vert', t_proftwp |
---|
237 | CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp |
---|
238 | : ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp |
---|
239 | : ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp |
---|
240 | : ,t_mod,q_mod,u_mod,v_mod,w_mod |
---|
241 | : ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc) |
---|
242 | ! write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod |
---|
243 | |
---|
244 | ! initial and boundary conditions : |
---|
245 | ! tsurf = ts_proftwp |
---|
246 | write(*,*) 'SST initiale: ',tsurf |
---|
247 | do l = 1, llm |
---|
248 | temp(l) = t_mod(l) |
---|
249 | q(l,1) = q_mod(l) |
---|
250 | q(l,2) = 0.0 |
---|
251 | u(l) = u_mod(l) |
---|
252 | v(l) = v_mod(l) |
---|
253 | omega(l) = w_mod(l) |
---|
254 | omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq |
---|
255 | |
---|
256 | alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) |
---|
257 | !on applique le forcage total au premier pas de temps |
---|
258 | !attention: signe different de toga |
---|
259 | d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l)) |
---|
260 | d_q_adv(l,1) = (hq_mod(l)+vq_mod(l)) |
---|
261 | d_q_adv(l,2) = 0.0 |
---|
262 | enddo |
---|
263 | |
---|
264 | endif !forcing_twpice |
---|
265 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
266 | !--------------------------------------------------------------------- |
---|
267 | ! Forcing from Arm_Cu case |
---|
268 | ! For this case, ifa_armcu.txt contains sensible, latent heat fluxes |
---|
269 | ! large scale advective forcing,radiative forcing |
---|
270 | ! and advective tendency of theta and qt to be applied |
---|
271 | !--------------------------------------------------------------------- |
---|
272 | |
---|
273 | if (forcing_armcu) then |
---|
274 | ! read armcu forcing : |
---|
275 | write(*,*) 'Avant lecture Forcing Arm_Cu' |
---|
276 | fich_armcu = './ifa_armcu.txt' |
---|
277 | CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu, |
---|
278 | : sens_armcu,flat_armcu,adv_theta_armcu, |
---|
279 | : rad_theta_armcu,adv_qt_armcu) |
---|
280 | write(*,*) 'Forcing Arm_Cu lu' |
---|
281 | |
---|
282 | !---------------------------------------------------------------------- |
---|
283 | ! Read profiles from file: prof.inp.19 or prof.inp.40 |
---|
284 | ! For this case, profiles are given for two vertical resolution |
---|
285 | ! 19 or 40 levels |
---|
286 | ! |
---|
287 | ! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html |
---|
288 | ! Note that the initial profiles contain no liquid water! |
---|
289 | ! (so potential temperature can be interpreted as liquid water |
---|
290 | ! potential temperature and water vapor as total water) |
---|
291 | ! profiles are given at full levels |
---|
292 | !---------------------------------------------------------------------- |
---|
293 | |
---|
294 | call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod, |
---|
295 | . v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp) |
---|
296 | |
---|
297 | ! time interpolation for initial conditions: |
---|
298 | write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1 |
---|
299 | |
---|
300 | print *,'Avant interp_armcu_time' |
---|
301 | print *,'daytime=',daytime |
---|
302 | print *,'day1=',day1 |
---|
303 | print *,'annee_ref=',annee_ref |
---|
304 | print *,'year_ini_armcu=',year_ini_armcu |
---|
305 | print *,'day_ju_ini_armcu=',day_ju_ini_armcu |
---|
306 | print *,'nt_armcu=',nt_armcu |
---|
307 | print *,'dt_armcu=',dt_armcu |
---|
308 | print *,'nlev_armcu=',nlev_armcu |
---|
309 | CALL interp_armcu_time(daytime,day1,annee_ref |
---|
310 | i ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu |
---|
311 | i ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu |
---|
312 | i ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof |
---|
313 | i ,adv_theta_prof,rad_theta_prof,adv_qt_prof) |
---|
314 | write(*,*) 'Forcages interpoles dans temps' |
---|
315 | |
---|
316 | ! No vertical interpolation if nlev imposed to 19 or 40 |
---|
317 | ! The vertical grid stops at 4000m # 600hPa |
---|
318 | mxcalc=llm |
---|
319 | |
---|
320 | ! initial and boundary conditions : |
---|
321 | ! tsurf = ts_prof |
---|
322 | ! tsurf read in lmdz1d.def |
---|
323 | write(*,*) 'Tsurf initiale: ',tsurf |
---|
324 | do l = 1, llm |
---|
325 | play(l)=play_mod(l)*100. |
---|
326 | presnivs(l)=play(l) |
---|
327 | zlay(l)=height(l) |
---|
328 | temp(l) = t_mod(l) |
---|
329 | teta(l)=theta_mod(l) |
---|
330 | q(l,1) = qv_mod(l)/1000. |
---|
331 | ! No liquid water in the initial profil |
---|
332 | q(l,2) = 0. |
---|
333 | u(l) = u_mod(l) |
---|
334 | ug(l)= u_mod(l) |
---|
335 | v(l) = v_mod(l) |
---|
336 | vg(l)= v_mod(l) |
---|
337 | ! Advective forcings are given in K or g/kg ... per HOUR |
---|
338 | ! IF(height(l).LT.1000) THEN |
---|
339 | ! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600. |
---|
340 | ! d_q_adv(l,1) = adv_qt_prof/1000./3600. |
---|
341 | ! d_q_adv(l,2) = 0.0 |
---|
342 | ! ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN |
---|
343 | ! d_th_adv(l) = (adv_theta_prof + rad_theta_prof)* |
---|
344 | ! : (1-(height(l)-1000.)/2000.) |
---|
345 | ! d_th_adv(l) = d_th_adv(l)/3600. |
---|
346 | ! d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.) |
---|
347 | ! d_q_adv(l,1) = d_q_adv(l,1)/1000./3600. |
---|
348 | ! d_q_adv(l,2) = 0.0 |
---|
349 | ! ELSE |
---|
350 | ! d_th_adv(l) = 0.0 |
---|
351 | ! d_q_adv(l,1) = 0.0 |
---|
352 | ! d_q_adv(l,2) = 0.0 |
---|
353 | ! ENDIF |
---|
354 | enddo |
---|
355 | ! plev at half levels is given in proh.inp.19 or proh.inp.40 files |
---|
356 | plev(1)= ap(llm+1)+bp(llm+1)*psurf |
---|
357 | do l = 1, llm |
---|
358 | plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf |
---|
359 | print *,'Read_forc: l height play plev zlay temp', |
---|
360 | : l,height(l),play(l),plev(l),zlay(l),temp(l) |
---|
361 | enddo |
---|
362 | ! For this case, fluxes are imposed |
---|
363 | fsens=-1*sens_prof |
---|
364 | flat=-1*flat_prof |
---|
365 | |
---|
366 | endif ! forcing_armcu |
---|
367 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
368 | |
---|