source: LMDZ5/branches/IPSLCM6.0.10/libf/phylmd/dyn1d/1D_interp_cases.h @ 5445

Last change on this file since 5445 was 2720, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2664:2719 into testing branch

  • Property svn:keywords set to Id
File size: 38.3 KB
Line 
1!
2! $Id: 1D_interp_cases.h 2720 2016-11-30 12:28:41Z fhourdin $
3!
4!---------------------------------------------------------------------
5! Forcing_LES case: constant dq_dyn
6!---------------------------------------------------------------------
7      if (forcing_LES) then
8        DO l = 1,llm
9          d_q_adv(l,1) = dq_dyn(l,1)
10        ENDDO
11      endif ! forcing_LES
12!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13!---------------------------------------------------------------------
14! Interpolation forcing in time and onto model levels
15!---------------------------------------------------------------------
16      if (forcing_GCSSold) then
17
18       call get_uvd(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,              &
19     &               ht_gcssold,hq_gcssold,hw_gcssold,                          &
20     &               hu_gcssold,hv_gcssold,                                     &
21     &               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,                 &
22     &               imp_fcg_gcssold,ts_fcg_gcssold,                            &
23     &               Tp_fcg_gcssold,Turb_fcg_gcssold)
24       if (prt_level.ge.1) then
25         print *,' get_uvd -> hqturb_gcssold ',it,hqturb_gcssold
26       endif
27! large-scale forcing :
28!!!      tsurf = ts_gcssold
29      do l = 1, llm
30!       u(l) = hu_gcssold(l) !  on prescrit le vent
31!       v(l) = hv_gcssold(l)    !  on prescrit le vent
32!       omega(l) = hw_gcssold(l)
33!       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
34!       omega2(l)=-rho(l)*omega(l)
35       omega(l) = hw_gcssold(l)
36       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
37
38       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
39       d_th_adv(l) = ht_gcssold(l)
40       d_q_adv(l,1) = hq_gcssold(l)
41       dt_cooling(l) = 0.0
42      enddo
43
44      endif ! forcing_GCSSold
45!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
46!---------------------------------------------------------------------
47! Interpolation Toga forcing
48!---------------------------------------------------------------------
49      if (forcing_toga) then
50
51       if (prt_level.ge.1) then
52        print*,                                                             &
53     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_toga=',     &
54     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_toga
55       endif
56
57! time interpolation:
58        CALL interp_toga_time(daytime,day1,annee_ref                        &
59     &             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga           &
60     &             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga        &
61     &             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga           &
62     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof    &
63     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
64
65        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
66
67! vertical interpolation:
68      CALL interp_toga_vertical(play,nlev_toga,plev_prof                    &
69     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                          &
70     &         ,ht_prof,vt_prof,hq_prof,vq_prof                             &
71     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
72     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
73
74! large-scale forcing :
75      tsurf = ts_prof
76      do l = 1, llm
77       u(l) = u_mod(l) ! sb: on prescrit le vent
78       v(l) = v_mod(l) ! sb: on prescrit le vent
79!       omega(l) = w_prof(l)
80!       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
81!       omega2(l)=-rho(l)*omega(l)
82       omega(l) = w_mod(l)
83       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
84
85       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
86       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
87       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
88       dt_cooling(l) = 0.0
89      enddo
90
91      endif ! forcing_toga
92!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
93! Interpolation DICE forcing
94!---------------------------------------------------------------------
95      if (forcing_dice) then
96
97       if (prt_level.ge.1) then
98        print*,'#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_dice=',&
99     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_dice
100       endif
101
102! time interpolation:
103      CALL interp_dice_time(daytime,day1,annee_ref                    &
104     &             ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice     & 
105     &             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice   &
106     &             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice     &
107     &             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &
108     &             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof     &
109     &             ,ustar_prof,psurf_prof,ug_profd,vg_profd           &
110     &             ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd       &
111     &             ,omega_profd)
112!     do l = 1, llm
113!     print *,'llm l omega_profd',llm,l,omega_profd(l)
114!     enddo
115
116        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
117
118! vertical interpolation:
119      CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice        &
120     &         ,t_dice,qv_dice,u_dice,v_dice,o3_dice                   &
121     &         ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd,omega_profd &
122     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                        &
123     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
124!     do l = 1, llm
125!      print *,'llm l omega_mod',llm,l,omega_mod(l)
126!     enddo
127
128! Les forcages DICE sont donnes /jour et non /seconde !
129      ht_mod(:)=ht_mod(:)/86400.
130      hq_mod(:)=hq_mod(:)/86400.
131      hu_mod(:)=hu_mod(:)/86400.
132      hv_mod(:)=hv_mod(:)/86400.
133
134!calcul de l'advection verticale a partir du omega (repris cas TWPICE, MPL 05082013)
135!Calcul des gradients verticaux
136!initialisation
137      d_t_z(:)=0.
138      d_q_z(:)=0.
139      d_u_z(:)=0.
140      d_v_z(:)=0.
141      DO l=2,llm-1
142       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
143       d_q_z(l)=(q(l+1,1)-q(l-1,1)) /(play(l+1)-play(l-1))
144       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
145       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
146      ENDDO
147      d_t_z(1)=d_t_z(2)
148      d_q_z(1)=d_q_z(2)
149!     d_u_z(1)=u(2)/(play(2)-psurf)/5.
150!     d_v_z(1)=v(2)/(play(2)-psurf)/5.
151      d_u_z(1)=0.
152      d_v_z(1)=0.
153      d_t_z(llm)=d_t_z(llm-1)
154      d_q_z(llm)=d_q_z(llm-1)
155      d_u_z(llm)=d_u_z(llm-1)
156      d_v_z(llm)=d_v_z(llm-1)
157
158!Calcul de l advection verticale: 
159! utiliser omega (Pa/s) et non w (m/s) !! MP 20131108
160      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
161      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
162      d_u_dyn_z(:)=omega_mod(:)*d_u_z(:)
163      d_v_dyn_z(:)=omega_mod(:)*d_v_z(:)
164
165! large-scale forcing :
166!     tsurf = tg_prof    MPL 20130925 commente
167      psurf = psurf_prof
168! For this case, fluxes are imposed
169      fsens=-1*shf_prof
170      flat=-1*lhf_prof
171      ust=ustar_prof
172      tg=tg_prof
173      print *,'ust= ',ust
174      do l = 1, llm
175       ug(l)= ug_profd
176       vg(l)= vg_profd
177!       omega(l) = w_prof(l)
178!      rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
179!       omega2(l)=-rho(l)*omega(l)
180!      omega(l) = w_mod(l)*(-rg*rho(l))
181       omega(l) = omega_mod(l)
182       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
183
184       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
185       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
186       d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
187       d_u_adv(l) = hu_mod(l)-d_u_dyn_z(l)
188       d_v_adv(l) = hv_mod(l)-d_v_dyn_z(l)
189       dt_cooling(l) = 0.0
190      enddo
191
192      endif ! forcing_dice
193!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194! Interpolation gabls4 forcing
195!---------------------------------------------------------------------
196      if (forcing_gabls4 ) then
197
198       if (prt_level.ge.1) then
199        print*,'#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_gabls4=',&
200     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_gabls4
201       endif
202
203! time interpolation:
204      CALL interp_gabls4_time(daytime,day1,annee_ref                                     &
205     &             ,year_ini_gabls4,day_ju_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4  & 
206     &             ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4                            &
207     &             ,ug_profg,vg_profg,ht_profg,hq_profg,tg_profg)
208
209        if (type_ts_forcing.eq.1) ts_cur = tg_prof ! SST used in read_tsurf1d
210
211! vertical interpolation:
212! on re-utilise le programme interp_dice_vertical: les transformations sur
213! plev_gabls4,th_gabls4,qv_gabls4,u_gabls4,v_gabls4 ne sont pas prises en compte.
214! seules celles sur ht_profg,hq_profg,ug_profg,vg_profg sont prises en compte.
215
216      CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4         &
217!    &         ,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,poub            &
218     &         ,poub,poub,poub,poub,poub                             &
219     &         ,ht_profg,hq_profg,ug_profg,vg_profg,poub,poub        &
220     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                      &
221     &         ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
222
223      do l = 1, llm
224       ug(l)= ug_mod(l)
225       vg(l)= vg_mod(l)
226       d_th_adv(l)=ht_mod(l)
227       d_q_adv(l,1)=hq_mod(l)
228      enddo
229
230      endif ! forcing_gabls4
231!---------------------------------------------------------------------
232
233!---------------------------------------------------------------------
234!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
235!---------------------------------------------------------------------
236! Interpolation forcing TWPice
237!---------------------------------------------------------------------
238      if (forcing_twpice) then
239
240        print*,                                                             &
241     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_twpi=',     &
242     &    daytime,day1,(daytime-day1)*86400.,                               &
243     &    (daytime-day1)*86400/dt_twpi
244
245! time interpolation:
246        CALL interp_toga_time(daytime,day1,annee_ref                        &
247     &       ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi       &
248     &       ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi          &
249     &       ,ht_twpi,vt_twpi,hq_twpi,vq_twpi                               &
250     &       ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp,u_proftwp         &
251     &       ,v_proftwp,w_proftwp                                           &
252     &       ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
253
254! vertical interpolation:
255      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp                 &
256     &         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp           &
257     &         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp                 &
258     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
259     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
260
261
262!calcul de l'advection verticale a partir du omega
263!Calcul des gradients verticaux
264!initialisation
265      d_t_z(:)=0.
266      d_q_z(:)=0.
267      d_t_dyn_z(:)=0.
268      d_q_dyn_z(:)=0.
269      DO l=2,llm-1
270       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
271       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
272      ENDDO
273      d_t_z(1)=d_t_z(2)
274      d_q_z(1)=d_q_z(2)
275      d_t_z(llm)=d_t_z(llm-1)
276      d_q_z(llm)=d_q_z(llm-1)
277
278!Calcul de l advection verticale
279      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
280      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
281
282!wind nudging above 500m with a 2h time scale
283        do l=1,llm
284        if (nudge_wind) then
285!           if (phi(l).gt.5000.) then
286        if (phi(l).gt.0.) then
287        u(l)=u(l)+timestep*(u_mod(l)-u(l))/(2.*3600.)
288        v(l)=v(l)+timestep*(v_mod(l)-v(l))/(2.*3600.)
289           endif   
290        else
291        u(l) = u_mod(l) 
292        v(l) = v_mod(l)
293        endif
294        enddo
295
296!CR:nudging of q and theta with a 6h time scale above 15km
297        if (nudge_thermo) then
298        do l=1,llm
299           zz(l)=phi(l)/9.8
300           if ((zz(l).le.16000.).and.(zz(l).gt.15000.)) then
301             zfact=(zz(l)-15000.)/1000.
302        q(l,1)=q(l,1)+timestep*(q_mod(l)-q(l,1))/(6.*3600.)*zfact
303        temp(l)=temp(l)+timestep*(t_mod(l)-temp(l))/(6.*3600.)*zfact
304           else if (zz(l).gt.16000.) then
305        q(l,1)=q(l,1)+timestep*(q_mod(l)-q(l,1))/(6.*3600.)
306        temp(l)=temp(l)+timestep*(t_mod(l)-temp(l))/(6.*3600.)
307           endif
308        enddo   
309        endif
310
311      do l = 1, llm
312       omega(l) = w_mod(l)
313       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
314       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
315!calcul de l'advection totale
316        if (cptadvw) then
317        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
318!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
319        d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
320!        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
321        else
322        d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
323        d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
324        endif
325       dt_cooling(l) = 0.0
326      enddo
327
328      endif ! forcing_twpice
329
330!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
331!---------------------------------------------------------------------
332! Interpolation forcing AMMA
333!---------------------------------------------------------------------
334
335       if (forcing_amma) then
336
337        print*,                                                             &
338     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_amma=',     &
339     &    daytime,day1,(daytime-day1)*86400.,                               &
340     &    (daytime-day1)*86400/dt_amma
341
342! time interpolation using TOGA interpolation routine
343        CALL interp_amma_time(daytime,day1,annee_ref                        &
344     &       ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma       &
345     &       ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma                  &
346     &       ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma            &
347     &       ,sens_profamma)
348
349      print*,'apres interpolation temporelle AMMA'
350
351      do k=1,nlev_amma
352         th_profamma(k)=0.
353         q_profamma(k)=0.
354         u_profamma(k)=0.
355         v_profamma(k)=0.
356         vt_profamma(k)=0.
357         vq_profamma(k)=0.
358       enddo
359! vertical interpolation using TOGA interpolation routine:
360!      write(*,*)'avant interp vert', t_proftwp
361      CALL interp_toga_vertical(play,nlev_amma,plev_amma                      &
362     &         ,th_profamma,q_profamma,u_profamma,v_profamma                 &
363     &         ,vitw_profamma                                               &
364     &         ,ht_profamma,vt_profamma,hq_profamma,vq_profamma             &
365     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
366     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
367       write(*,*) 'Profil initial forcing AMMA interpole'
368
369
370!calcul de l'advection verticale a partir du omega
371!Calcul des gradients verticaux
372!initialisation
373      do l=1,llm
374      d_t_z(l)=0.
375      d_q_z(l)=0.
376      enddo
377
378      DO l=2,llm-1
379       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
380       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
381      ENDDO
382      d_t_z(1)=d_t_z(2)
383      d_q_z(1)=d_q_z(2)
384      d_t_z(llm)=d_t_z(llm-1)
385      d_q_z(llm)=d_q_z(llm-1)
386
387
388      do l = 1, llm
389       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
390       omega(l) = w_mod(l)*(-rg*rho(l))
391       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
392       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
393!calcul de l'advection totale
394!        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
395!attention: on impose dth
396        d_th_adv(l) = alpha*omega(l)/rcpd+                                  &
397     &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
398!        d_th_adv(l) = 0.
399!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
400        d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l)
401!        d_q_adv(l,1) = 0.
402!        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
403   
404       dt_cooling(l) = 0.0
405      enddo
406
407
408!     ok_flux_surf=.false.
409      fsens=-1.*sens_profamma
410      flat=-1.*lat_profamma
411
412      endif ! forcing_amma
413
414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415!---------------------------------------------------------------------
416! Interpolation forcing Rico
417!---------------------------------------------------------------------
418      if (forcing_rico) then
419!       call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,
420!     :  q,temp,u,v,play)
421       call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play)
422
423        do l=1,llm
424       d_th_adv(l) =  (dth_rico(l) +  dt_dyn(l))
425       d_q_adv(l,1) = (dqh_rico(l) +  dq_dyn(l,1))
426       d_q_adv(l,2) = 0.
427        enddo
428      endif  ! forcing_rico
429!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
430!---------------------------------------------------------------------
431! Interpolation forcing Arm_cu
432!---------------------------------------------------------------------
433      if (forcing_armcu) then
434
435        print*,                                                             &
436     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_armcu=',    &
437     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_armcu
438
439! time interpolation:
440! ATTENTION, cet appel ne convient pas pour TOGA !!
441! revoir 1DUTILS.h et les arguments
442      CALL interp_armcu_time(daytime,day1,annee_ref                         &
443     &            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu        &
444     &            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu          &
445     &            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof         &
446     &            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
447
448! vertical interpolation:
449! No vertical interpolation if nlev imposed to 19 or 40
450
451! For this case, fluxes are imposed
452       fsens=-1*sens_prof
453       flat=-1*flat_prof
454
455! Advective forcings are given in K or g/kg ... BY HOUR
456      do l = 1, llm
457       ug(l)= u_mod(l)
458       vg(l)= v_mod(l)
459       IF((phi(l)/RG).LT.1000) THEN
460         d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
461         d_q_adv(l,1) = adv_qt_prof/1000./3600.
462         d_q_adv(l,2) = 0.0
463!        print *,'INF1000: phi dth dq1 dq2',
464!    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
465       ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN
466         fact=((phi(l)/RG)-1000.)/2000.
467         fact=1-fact
468         d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
469         d_q_adv(l,1) = adv_qt_prof*fact/1000./3600.
470         d_q_adv(l,2) = 0.0
471!        print *,'SUP1000: phi fact dth dq1 dq2',
472!    :  phi(l)/RG,fact,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
473       ELSE
474         d_th_adv(l) = 0.0
475         d_q_adv(l,1) = 0.0
476         d_q_adv(l,2) = 0.0
477!        print *,'SUP3000: phi dth dq1 dq2',
478!    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
479       ENDIF
480      dt_cooling(l) = 0.0 
481!     print *,'Interp armcu: phi dth dq1 dq2',
482!    :  l,phi(l),d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
483      enddo
484      endif ! forcing_armcu
485!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
486!---------------------------------------------------------------------
487! Interpolation forcing in time and onto model levels
488!---------------------------------------------------------------------
489      if (forcing_sandu) then
490
491        print*,                                                             &
492     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=',    &
493     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu
494
495! time interpolation:
496! ATTENTION, cet appel ne convient pas pour TOGA !!
497! revoir 1DUTILS.h et les arguments
498      CALL interp_sandu_time(daytime,day1,annee_ref                         &
499     &             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu       &
500     &             ,nlev_sandu                                              &
501     &             ,ts_sandu,ts_prof)
502
503        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
504
505! vertical interpolation:
506      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs                 &
507     &         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs           &
508     &         ,omega_profs,o3mmr_profs                                     &
509     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                       &
510     &         ,omega_mod,o3mmr_mod,mxcalc)
511!calcul de l'advection verticale
512!Calcul des gradients verticaux
513!initialisation
514      d_t_z(:)=0.
515      d_q_z(:)=0.
516      d_t_dyn_z(:)=0.
517      d_q_dyn_z(:)=0.
518! schema centre
519!     DO l=2,llm-1
520!      d_t_z(l)=(temp(l+1)-temp(l-1))
521!    &          /(play(l+1)-play(l-1))
522!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
523!    &          /(play(l+1)-play(l-1))
524! schema amont
525      DO l=2,llm-1
526       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
527       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
528!     print *,'l temp2 temp0 play2 play0 omega_mod',
529!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
530      ENDDO
531      d_t_z(1)=d_t_z(2)
532      d_q_z(1)=d_q_z(2)
533      d_t_z(llm)=d_t_z(llm-1)
534      d_q_z(llm)=d_q_z(llm-1)
535
536!  calcul de l advection verticale
537! Confusion w (m/s) et omega (Pa/s) !!
538      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
539      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
540!     do l=1,llm
541!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
542!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
543!     enddo
544
545
546! large-scale forcing : pour le cas Sandu ces forcages sont la SST
547! et une divergence constante -> profil de omega
548      tsurf = ts_prof
549      write(*,*) 'SST suivante: ',tsurf
550      do l = 1, llm
551       omega(l) = omega_mod(l)
552       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
553
554       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
555!
556!      d_th_adv(l) = 0.0
557!      d_q_adv(l,1) = 0.0
558!CR:test advection=0
559!calcul de l'advection verticale
560        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
561!        print*,'temp adv',l,-d_t_dyn_z(l)
562        d_q_adv(l,1) = -d_q_dyn_z(l)
563!        print*,'q adv',l,-d_q_dyn_z(l)
564       dt_cooling(l) = 0.0
565      enddo
566      endif ! forcing_sandu
567!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
568!---------------------------------------------------------------------
569! Interpolation forcing in time and onto model levels
570!---------------------------------------------------------------------
571      if (forcing_astex) then
572
573        print*,                                                             &
574     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=',    &
575     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex
576
577! time interpolation:
578! ATTENTION, cet appel ne convient pas pour TOGA !!
579! revoir 1DUTILS.h et les arguments
580      CALL interp_astex_time(daytime,day1,annee_ref                         &
581     &             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex       &
582     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex         &
583     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
584     &             ,ufa_prof,vfa_prof)
585
586        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
587
588! vertical interpolation:
589      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
590     &         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa                &
591     &         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa               &
592     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod        &
593     &         ,tke_mod,o3mmr_mod,mxcalc)
594!calcul de l'advection verticale
595!Calcul des gradients verticaux
596!initialisation
597      d_t_z(:)=0.
598      d_q_z(:)=0.
599      d_t_dyn_z(:)=0.
600      d_q_dyn_z(:)=0.
601! schema centre
602!     DO l=2,llm-1
603!      d_t_z(l)=(temp(l+1)-temp(l-1))
604!    &          /(play(l+1)-play(l-1))
605!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
606!    &          /(play(l+1)-play(l-1))
607! schema amont
608      DO l=2,llm-1
609       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
610       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
611!     print *,'l temp2 temp0 play2 play0 omega_mod',
612!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
613      ENDDO
614      d_t_z(1)=d_t_z(2)
615      d_q_z(1)=d_q_z(2)
616      d_t_z(llm)=d_t_z(llm-1)
617      d_q_z(llm)=d_q_z(llm-1)
618
619!  calcul de l advection verticale
620! Confusion w (m/s) et omega (Pa/s) !!
621      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
622      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
623!     do l=1,llm
624!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
625!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
626!     enddo
627
628
629! large-scale forcing : pour le cas Astex ces forcages sont la SST
630! la divergence,ug,vg,ufa,vfa
631      tsurf = ts_prof
632      write(*,*) 'SST suivante: ',tsurf
633      do l = 1, llm
634       omega(l) = w_mod(l)
635       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
636
637       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
638!
639!      d_th_adv(l) = 0.0
640!      d_q_adv(l,1) = 0.0
641!CR:test advection=0
642!calcul de l'advection verticale
643        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
644!        print*,'temp adv',l,-d_t_dyn_z(l)
645        d_q_adv(l,1) = -d_q_dyn_z(l)
646!        print*,'q adv',l,-d_q_dyn_z(l)
647       dt_cooling(l) = 0.0
648      enddo
649      endif ! forcing_astex
650
651!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
652!---------------------------------------------------------------------
653! Interpolation forcing standard case
654!---------------------------------------------------------------------
655      if (forcing_case) then
656
657        print*,                                                             &
658     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
659     &    daytime,day1,(daytime-day1)*86400.,                               &
660     &    (daytime-day1)*86400/pdt_cas
661
662! time interpolation:
663        CALL interp_case_time(daytime,day1,annee_ref                                        &
664!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
665     &       ,nt_cas,nlev_cas                                                               &
666     &       ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas,ug_cas,vg_cas                         &
667     &       ,vitw_cas,du_cas,hu_cas,vu_cas                                                 &
668     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
669     &       ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas                               &
670     &       ,uw_cas,vw_cas,q1_cas,q2_cas                                                   &
671     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
672     &       ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
673     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
674     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
675     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
676
677             ts_cur = ts_prof_cas
678             psurf=plev_prof_cas(1)
679
680! vertical interpolation:
681      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
682     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
683     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
684     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
685     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
686     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
687     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
688
689
690!calcul de l'advection verticale a partir du omega
691!Calcul des gradients verticaux
692!initialisation
693      d_t_z(:)=0.
694      d_q_z(:)=0.
695      d_u_z(:)=0.
696      d_v_z(:)=0.
697      d_t_dyn_z(:)=0.
698      d_q_dyn_z(:)=0.
699      d_u_dyn_z(:)=0.
700      d_v_dyn_z(:)=0.
701      DO l=2,llm-1
702       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
703       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
704       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
705       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
706      ENDDO
707      d_t_z(1)=d_t_z(2)
708      d_q_z(1)=d_q_z(2)
709      d_u_z(1)=d_u_z(2)
710      d_v_z(1)=d_v_z(2)
711      d_t_z(llm)=d_t_z(llm-1)
712      d_q_z(llm)=d_q_z(llm-1)
713      d_u_z(llm)=d_u_z(llm-1)
714      d_v_z(llm)=d_v_z(llm-1)
715
716!Calcul de l advection verticale
717      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
718      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
719      d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
720      d_v_dyn_z(:)=w_mod_cas(:)*d_v_z(:)
721
722!wind nudging
723      if (nudge_u.gt.0.) then
724        do l=1,llm
725           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
726        enddo
727      else
728        do l=1,llm
729        u(l) = u_mod_cas(l) 
730        enddo
731      endif
732
733      if (nudge_v.gt.0.) then
734        do l=1,llm
735           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
736        enddo
737      else
738        do l=1,llm
739        v(l) = v_mod_cas(l) 
740        enddo
741      endif
742
743      if (nudge_w.gt.0.) then
744        do l=1,llm
745           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
746        enddo
747      else
748        do l=1,llm
749        w(l) = w_mod_cas(l) 
750        enddo
751      endif
752
753!nudging of q and temp
754      if (nudge_t.gt.0.) then
755        do l=1,llm
756           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
757        enddo
758      endif
759      if (nudge_q.gt.0.) then
760        do l=1,llm
761           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
762        enddo
763      endif
764
765      do l = 1, llm
766       omega(l) = w_mod_cas(l)
767       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
768       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
769
770!calcul advection
771        if ((tend_u.eq.1).and.(tend_w.eq.0)) then
772           d_u_adv(l)=du_mod_cas(l)
773        else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
774           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
775        endif
776
777        if ((tend_v.eq.1).and.(tend_w.eq.0)) then
778           d_v_adv(l)=dv_mod_cas(l)
779        else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
780           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
781        endif
782
783        if ((tend_t.eq.1).and.(tend_w.eq.0)) then
784!           d_th_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
785           d_th_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
786        else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
787!           d_th_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
788           d_th_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
789        endif
790
791        if ((tend_q.eq.1).and.(tend_w.eq.0)) then
792!           d_q_adv(l,1)=dq_mod_cas(l)
793           d_q_adv(l,1)=-1*dq_mod_cas(l)
794        else if ((tend_q.eq.1).and.(tend_w.eq.1)) then
795!           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
796           d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
797        endif
798         
799        if (tend_rayo.eq.1) then
800           dt_cooling(l) = dtrad_mod_cas(l)
801!          print *,'dt_cooling=',dt_cooling(l)
802        else
803           dt_cooling(l) = 0.0
804        endif
805      enddo
806
807! Faut-il multiplier par -1 ? (MPL 20160713)
808      IF(ok_flux_surf) THEN
809       fsens=sens_prof_cas
810       flat=lat_prof_cas
811      ENDIF
812!
813      IF (ok_prescr_ust) THEN
814       ust=ustar_prof_cas
815       print *,'ust=',ust
816      ENDIF
817      endif ! forcing_case
818
819
820!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
821!---------------------------------------------------------------------
822! Interpolation forcing standard case
823!---------------------------------------------------------------------
824      if (forcing_case2) then
825
826        print*,                                                             &
827     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
828     &    daytime,day1,(daytime-day1)*86400.,                               &
829     &    (daytime-day1)*86400/pdt_cas
830
831! time interpolation:
832        CALL interp2_case_time(daytime,day1,annee_ref                                       &
833!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
834     &       ,nt_cas,nlev_cas                                                               &
835     &       ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas      &
836     &       ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas             &
837     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
838     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
839     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
840!
841     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
842     &       ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas                              &
843     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas    &
844     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
845     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
846     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
847     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
848     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
849
850             ts_cur = ts_prof_cas
851!            psurf=plev_prof_cas(1)
852             psurf=ps_prof_cas
853
854! vertical interpolation:
855      CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas                                              &
856     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
857     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
858     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                      &
859     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
860     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
861     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                    &
862!
863     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas          &
864     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                         &
865     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
866     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
867     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
868
869
870      DO l=1,llm
871      teta(l)=temp(l)*(100000./play(l))**(rd/rcpd)
872      ENDDO
873!calcul de l'advection verticale a partir du omega
874!Calcul des gradients verticaux
875!initialisation
876      d_t_z(:)=0.
877      d_th_z(:)=0.
878      d_q_z(:)=0.
879      d_t_dyn_z(:)=0.
880      d_th_dyn_z(:)=0.
881      d_q_dyn_z(:)=0.
882      DO l=2,llm-1
883       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
884       d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
885       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
886      ENDDO
887      d_t_z(1)=d_t_z(2)
888      d_th_z(1)=d_th_z(2)
889      d_q_z(1)=d_q_z(2)
890      d_t_z(llm)=d_t_z(llm-1)
891      d_th_z(llm)=d_th_z(llm-1)
892      d_q_z(llm)=d_q_z(llm-1)
893
894!Calcul de l advection verticale
895      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
896      d_th_dyn_z(:)=w_mod_cas(:)*d_th_z(:)
897      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
898
899!wind nudging
900      if (nudging_u.gt.0.) then
901        do l=1,llm
902           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
903        enddo
904      else
905        do l=1,llm
906        ug(l) = u_mod_cas(l) 
907        enddo
908      endif
909
910      if (nudging_v.gt.0.) then
911        do l=1,llm
912           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
913        enddo
914      else
915        do l=1,llm
916        vg(l) = v_mod_cas(l) 
917        enddo
918      endif
919
920      if (nudging_w.gt.0.) then
921        do l=1,llm
922           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
923        enddo
924      else
925        do l=1,llm
926        w(l) = w_mod_cas(l) 
927        enddo
928      endif
929
930!nudging of q and temp
931      if (nudging_t.gt.0.) then
932        do l=1,llm
933           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
934        enddo
935      endif
936      if (nudging_q.gt.0.) then
937        do l=1,llm
938           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
939        enddo
940      endif
941
942      do l = 1, llm
943       omega(l) = w_mod_cas(l)
944       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
945       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
946
947!calcul advection
948!       if ((tend_u.eq.1).and.(tend_w.eq.0)) then
949!          d_u_adv(l)=du_mod_cas(l)
950!       else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
951!          d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
952!       endif
953!
954!       if ((tend_v.eq.1).and.(tend_w.eq.0)) then
955!          d_v_adv(l)=dv_mod_cas(l)
956!       else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
957!          d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
958!       endif
959!
960!-----------------------------------------------------
961        if (tadv.eq.1 .or. tadvh.eq.1) then
962           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
963        else if (tadvv.eq.1) then
964! ATTENTION d_t_dyn_z pas calcule (voir twpice)
965           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
966        endif
967        print *,'interp_case d_t_dyn_z=',d_t_dyn_z(l),d_q_dyn_z(l)
968
969! Verifier le signe !!
970        if (thadv.eq.1 .or. thadvh.eq.1) then
971           d_th_adv(l)=dth_mod_cas(l)
972           print *,'dthadv=',d_th_adv(l)*86400.
973        else if (thadvv.eq.1) then
974           d_th_adv(l)=hth_mod_cas(l)-d_th_dyn_z(l)
975        endif
976 
977! Verifier le signe !!
978        if ((qadv.eq.1).and.(forc_w.eq.0)) then
979           d_q_adv(l,1)=dq_mod_cas(l)
980        else if ((qadvh.eq.1).and.(forc_w.eq.1)) then
981           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
982        endif
983         
984        if (trad.eq.1) then
985           tend_rayo=1
986           dt_cooling(l) = dtrad_mod_cas(l)
987!          print *,'dt_cooling=',dt_cooling(l)
988        else
989           dt_cooling(l) = 0.0
990        endif
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       print *,'1D_interp: sens,flat',fsens,flat
998      ENDIF
999!
1000      IF (ok_prescr_ust) THEN
1001       ust=ustar_prof_cas
1002       print *,'ust=',ust
1003      ENDIF
1004      endif ! forcing_case2
1005!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006
Note: See TracBrowser for help on using the repository browser.