source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phy1d/1D_interp_cases.h @ 4400

Last change on this file since 4400 was 1780, checked in by Laurent Fairhead, 11 years ago

Inclusion des cas AMMA et Sandu dans la configuration 1D
MPL


AMMA and Sandu cases included in 1D configuration
MPL

File size: 17.8 KB
Line 
1!---------------------------------------------------------------------
2! Interpolation forcing in time and onto model levels
3!---------------------------------------------------------------------
4      if (forcing_GCSSold) then
5
6       call get_uvd(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,
7     :               ht_gcssold,hq_gcssold,hw_gcssold,
8     :               hu_gcssold,hv_gcssold,
9     :               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,
10     :               imp_fcg_gcssold,ts_fcg_gcssold,
11     :               Tp_fcg_gcssold,Turb_fcg_gcssold)
12       if (prt_level.ge.1) then
13         print *,' get_uvd -> hqturb_gcssold ',it,hqturb_gcssold
14       endif
15! large-scale forcing :
16!!!      tsurf = ts_gcssold
17      do l = 1, llm
18!       u(l) = hu_gcssold(l) !  on prescrit le vent
19!       v(l) = hv_gcssold(l)    !  on prescrit le vent
20!       omega(l) = hw_gcssold(l)
21!       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
22!       omega2(l)=-rho(l)*omega(l)
23       omega(l) = hw_gcssold(l)
24       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
25
26       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
27       d_th_adv(l) = ht_gcssold(l)
28       d_q_adv(l,1) = hq_gcssold(l)
29       dt_cooling(l) = 0.0
30      enddo
31
32      endif ! forcing_GCSSold
33!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34!---------------------------------------------------------------------
35! Interpolation Toga forcing
36!---------------------------------------------------------------------
37      if (forcing_toga) then
38
39       if (prt_level.ge.1) then
40        print*,
41     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_toga=',
42     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_toga
43       endif
44
45! time interpolation:
46        CALL interp_toga_time(daytime,day1,annee_ref
47     i             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga
48     i             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga
49     i             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga
50     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
51     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
52
53        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
54
55! vertical interpolation:
56      CALL interp_toga_vertical(play,nlev_toga,plev_prof
57     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
58     :         ,ht_prof,vt_prof,hq_prof,vq_prof
59     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
60     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
61
62! large-scale forcing :
63      tsurf = ts_prof
64      do l = 1, llm
65       u(l) = u_mod(l) ! sb: on prescrit le vent
66       v(l) = v_mod(l) ! sb: on prescrit le vent
67!       omega(l) = w_prof(l)
68!       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
69!       omega2(l)=-rho(l)*omega(l)
70       omega(l) = w_mod(l)
71       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
72
73       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
74       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
75       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
76       dt_cooling(l) = 0.0
77      enddo
78
79      endif ! forcing_toga
80!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
81!---------------------------------------------------------------------
82! Interpolation forcing TWPice
83!---------------------------------------------------------------------
84      if (forcing_twpice) then
85
86        print*,
87     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_twpi=',
88     :    daytime,day1,(daytime-day1)*86400.,
89     :    (daytime-day1)*86400/dt_twpi
90
91! time interpolation:
92        CALL interp_toga_time(daytime,day1,annee_ref
93     i       ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi
94     i       ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
95     i       ,ht_twpi,vt_twpi,hq_twpi,vq_twpi
96     o       ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp,u_proftwp
97     o       ,v_proftwp,w_proftwp
98     o       ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
99
100! vertical interpolation:
101      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp
102     :         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp
103     :         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp
104     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
105     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
106
107
108!calcul de l'advection verticale a partir du omega
109cCalcul des gradients verticaux
110cinitialisation
111      d_t_z(:)=0.
112      d_q_z(:)=0.
113      d_t_dyn_z(:)=0.
114      d_q_dyn_z(:)=0.
115      DO l=2,llm-1
116       d_t_z(l)=(temp(l+1)-temp(l-1))
117     &          /(play(l+1)-play(l-1))
118       d_q_z(l)=(q(l+1,1)-q(l-1,1))
119     &          /(play(l+1)-play(l-1))
120      ENDDO
121      d_t_z(1)=d_t_z(2)
122      d_q_z(1)=d_q_z(2)
123      d_t_z(llm)=d_t_z(llm-1)
124      d_q_z(llm)=d_q_z(llm-1)
125
126cCalcul de l advection verticale
127      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
128      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
129
130!wind nudging above 500m with a 2h time scale
131        do l=1,llm
132        if (nudge_wind) then
133!           if (phi(l).gt.5000.) then
134        if (phi(l).gt.0.) then
135        u(l)=u(l)
136     .      +timestep*(u_mod(l)-u(l))/(2.*3600.)
137        v(l)=v(l) 
138     .      +timestep*(v_mod(l)-v(l))/(2.*3600.)
139           endif   
140        else
141        u(l) = u_mod(l) 
142        v(l) = v_mod(l)
143        endif
144        enddo
145
146!CR:nudging of q and theta with a 6h time scale above 15km
147        if (nudge_thermo) then
148        do l=1,llm
149           zz(l)=phi(l)/9.8
150           if ((zz(l).le.16000.).and.(zz(l).gt.15000.)) then
151             zfact=(zz(l)-15000.)/1000.
152        q(l,1)=q(l,1)
153     .      +timestep*(q_mod(l)-q(l,1))/(6.*3600.)*zfact
154        temp(l)=temp(l) 
155     .      +timestep*(t_mod(l)-temp(l))/(6.*3600.)*zfact
156           else if (zz(l).gt.16000.) then
157        q(l,1)=q(l,1)
158     .      +timestep*(q_mod(l)-q(l,1))/(6.*3600.)
159        temp(l)=temp(l) 
160     .      +timestep*(t_mod(l)-temp(l))/(6.*3600.)
161           endif
162        enddo   
163        endif
164
165      do l = 1, llm
166       omega(l) = w_mod(l)
167       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
168       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
169!calcul de l'advection totale
170        if (cptadvw) then
171        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
172!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
173        d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
174!        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
175        else
176        d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
177        d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
178        endif
179       dt_cooling(l) = 0.0
180      enddo
181
182      endif ! forcing_twpice
183
184!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185!---------------------------------------------------------------------
186! Interpolation forcing AMMA
187!---------------------------------------------------------------------
188
189       if (forcing_amma) then
190
191        print*,
192     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_amma=',
193     :    daytime,day1,(daytime-day1)*86400.,
194     :    (daytime-day1)*86400/dt_amma
195
196! time interpolation using TOGA interpolation routine
197        CALL interp_amma_time(daytime,day1,annee_ref
198     i       ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma
199     i       ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
200     o       ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma
201     :       ,sens_profamma)
202
203      print*,'apres interpolation temporelle AMMA'
204
205      do k=1,nlev_amma
206         th_profamma(k)=0.
207         q_profamma(k)=0.
208         u_profamma(k)=0.
209         v_profamma(k)=0.
210         vt_profamma(k)=0.
211         vq_profamma(k)=0.
212       enddo
213! vertical interpolation using TOGA interpolation routine:
214!      write(*,*)'avant interp vert', t_proftwp
215      CALL interp_toga_vertical(play,nlev_amma,plev_amma
216     :         ,th_profamma,q_profamma,u_profamma,v_profamma
217     :         ,vitw_profamma
218     :         ,ht_profamma,vt_profamma,hq_profamma,vq_profamma
219     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
220     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
221       write(*,*) 'Profil initial forcing AMMA interpole'
222
223
224!calcul de l'advection verticale a partir du omega
225cCalcul des gradients verticaux
226cinitialisation
227      do l=1,llm
228      d_t_z(l)=0.
229      d_q_z(l)=0.
230      enddo
231
232      DO l=2,llm-1
233       d_t_z(l)=(temp(l+1)-temp(l-1))
234     &          /(play(l+1)-play(l-1))
235       d_q_z(l)=(q(l+1,1)-q(l-1,1))
236     &          /(play(l+1)-play(l-1))
237      ENDDO
238      d_t_z(1)=d_t_z(2)
239      d_q_z(1)=d_q_z(2)
240      d_t_z(llm)=d_t_z(llm-1)
241      d_q_z(llm)=d_q_z(llm-1)
242
243
244      do l = 1, llm
245       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
246       omega(l) = w_mod(l)*(-rg*rho(l))
247       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
248       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
249!calcul de l'advection totale
250!        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
251!attention: on impose dth
252        d_th_adv(l) = alpha*omega(l)/rcpd+
253     &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
254!        d_th_adv(l) = 0.
255!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
256        d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l)
257!        d_q_adv(l,1) = 0.
258!        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
259   
260       dt_cooling(l) = 0.0
261      enddo
262
263
264!     ok_flux_surf=.false.
265      fsens=-1.*sens_profamma
266      flat=-1.*lat_profamma
267
268      endif ! forcing_amma
269
270!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
271!---------------------------------------------------------------------
272! Interpolation forcing Rico
273!---------------------------------------------------------------------
274      if (forcing_rico) then
275!       call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,
276!     :  q,temp,u,v,play)
277       call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,
278     :  q,temp,u,v,play)
279
280        do l=1,llm
281       d_th_adv(l) =  (dth_rico(l) +  dt_dyn(l))
282       d_q_adv(l,1) = (dqh_rico(l) +  dq_dyn(l,1))
283       d_q_adv(l,2) = 0.
284        enddo
285      endif  ! forcing_rico
286!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
287!---------------------------------------------------------------------
288! Interpolation forcing Arm_cu
289!---------------------------------------------------------------------
290      if (forcing_armcu) then
291
292        print*,
293     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_armcu=',
294     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_armcu
295
296! time interpolation:
297! ATTENTION, cet appel ne convient pas pour TOGA !!
298! revoir 1DUTILS.h et les arguments
299      CALL interp_armcu_time(daytime,day1,annee_ref
300     i            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu
301     i            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu
302     i            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof
303     i            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
304
305! vertical interpolation:
306! No vertical interpolation if nlev imposed to 19 or 40
307
308! For this case, fluxes are imposed
309       fsens=-1*sens_prof
310       flat=-1*flat_prof
311
312! Advective forcings are given in K or g/kg ... BY HOUR
313      do l = 1, llm
314       ug(l)= u_mod(l)
315       vg(l)= v_mod(l)
316       IF((phi(l)/RG).LT.1000) THEN
317         d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
318         d_q_adv(l,1) = adv_qt_prof/1000./3600.
319         d_q_adv(l,2) = 0.0
320!        print *,'INF1000: phi dth dq1 dq2',
321!    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
322       ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN
323         fact=((phi(l)/RG)-1000.)/2000.
324         fact=1-fact
325         d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
326         d_q_adv(l,1) = adv_qt_prof*fact/1000./3600.
327         d_q_adv(l,2) = 0.0
328!        print *,'SUP1000: phi fact dth dq1 dq2',
329!    :  phi(l)/RG,fact,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
330       ELSE
331         d_th_adv(l) = 0.0
332         d_q_adv(l,1) = 0.0
333         d_q_adv(l,2) = 0.0
334!        print *,'SUP3000: phi dth dq1 dq2',
335!    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
336       ENDIF
337      dt_cooling(l) = 0.0 
338!     print *,'Interp armcu: phi dth dq1 dq2',
339!    :  l,phi(l),d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
340      enddo
341      endif ! forcing_armcu
342!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
343!---------------------------------------------------------------------
344! Interpolation forcing in time and onto model levels
345!---------------------------------------------------------------------
346      if (forcing_sandu) then
347
348        print*,
349     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=',
350     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu
351
352! time interpolation:
353! ATTENTION, cet appel ne convient pas pour TOGA !!
354! revoir 1DUTILS.h et les arguments
355      CALL interp_sandu_time(daytime,day1,annee_ref
356     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
357     i             ,nlev_sandu
358     i             ,ts_sandu,ts_prof)
359
360        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
361
362! vertical interpolation:
363      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
364     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
365     :         ,omega_profs,o3mmr_profs
366     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
367     :         ,omega_mod,o3mmr_mod,mxcalc)
368!calcul de l'advection verticale
369cCalcul des gradients verticaux
370cinitialisation
371      d_t_z(:)=0.
372      d_q_z(:)=0.
373      d_t_dyn_z(:)=0.
374      d_q_dyn_z(:)=0.
375! schema centre
376!     DO l=2,llm-1
377!      d_t_z(l)=(temp(l+1)-temp(l-1))
378!    &          /(play(l+1)-play(l-1))
379!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
380!    &          /(play(l+1)-play(l-1))
381! schema amont
382      DO l=2,llm-1
383       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
384       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
385!     print *,'l temp2 temp0 play2 play0 omega_mod',
386!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
387      ENDDO
388      d_t_z(1)=d_t_z(2)
389      d_q_z(1)=d_q_z(2)
390      d_t_z(llm)=d_t_z(llm-1)
391      d_q_z(llm)=d_q_z(llm-1)
392
393!  calcul de l advection verticale
394! Confusion w (m/s) et omega (Pa/s) !!
395      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
396      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
397!     do l=1,llm
398!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
399!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
400!     enddo
401
402
403! large-scale forcing : pour le cas Sandu ces forcages sont la SST
404! et une divergence constante -> profil de omega
405      tsurf = ts_prof
406      write(*,*) 'SST suivante: ',tsurf
407      do l = 1, llm
408       omega(l) = omega_mod(l)
409       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
410
411       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
412!
413!      d_th_adv(l) = 0.0
414!      d_q_adv(l,1) = 0.0
415!CR:test advection=0
416!calcul de l'advection verticale
417        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
418!        print*,'temp adv',l,-d_t_dyn_z(l)
419        d_q_adv(l,1) = -d_q_dyn_z(l)
420!        print*,'q adv',l,-d_q_dyn_z(l)
421       dt_cooling(l) = 0.0
422      enddo
423      endif ! forcing_sandu
424!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
425!---------------------------------------------------------------------
426! Interpolation forcing in time and onto model levels
427!---------------------------------------------------------------------
428      if (forcing_astex) then
429
430        print*,
431     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=',
432     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex
433
434! time interpolation:
435! ATTENTION, cet appel ne convient pas pour TOGA !!
436! revoir 1DUTILS.h et les arguments
437      CALL interp_astex_time(daytime,day1,annee_ref
438     i             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex
439     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
440     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
441     i             ,ufa_prof,vfa_prof)
442
443        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
444
445! vertical interpolation:
446      CALL interp_astex_vertical(play,nlev_astex,plev_profa
447     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
448     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
449     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
450     :         ,tke_mod,o3mmr_mod,mxcalc)
451!calcul de l'advection verticale
452!Calcul des gradients verticaux
453!initialisation
454      d_t_z(:)=0.
455      d_q_z(:)=0.
456      d_t_dyn_z(:)=0.
457      d_q_dyn_z(:)=0.
458! schema centre
459!     DO l=2,llm-1
460!      d_t_z(l)=(temp(l+1)-temp(l-1))
461!    &          /(play(l+1)-play(l-1))
462!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
463!    &          /(play(l+1)-play(l-1))
464! schema amont
465      DO l=2,llm-1
466       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
467       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
468!     print *,'l temp2 temp0 play2 play0 omega_mod',
469!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
470      ENDDO
471      d_t_z(1)=d_t_z(2)
472      d_q_z(1)=d_q_z(2)
473      d_t_z(llm)=d_t_z(llm-1)
474      d_q_z(llm)=d_q_z(llm-1)
475
476!  calcul de l advection verticale
477! Confusion w (m/s) et omega (Pa/s) !!
478      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
479      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
480!     do l=1,llm
481!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
482!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
483!     enddo
484
485
486! large-scale forcing : pour le cas Astex ces forcages sont la SST
487! la divergence,ug,vg,ufa,vfa
488      tsurf = ts_prof
489      write(*,*) 'SST suivante: ',tsurf
490      do l = 1, llm
491       omega(l) = w_mod(l)
492       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
493
494       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
495!
496!      d_th_adv(l) = 0.0
497!      d_q_adv(l,1) = 0.0
498!CR:test advection=0
499!calcul de l'advection verticale
500        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
501!        print*,'temp adv',l,-d_t_dyn_z(l)
502        d_q_adv(l,1) = -d_q_dyn_z(l)
503!        print*,'q adv',l,-d_q_dyn_z(l)
504       dt_cooling(l) = 0.0
505      enddo
506      endif ! forcing_astex
507!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
508
Note: See TracBrowser for help on using the repository browser.