source: LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h @ 5006

Last change on this file since 5006 was 2920, checked in by fhourdin, 8 years ago

Modification 1D pour le format standard (Marie-Pierre)

  • Property svn:keywords set to Id
File size: 39.2 KB
Line 
1!
2! $Id: 1D_interp_cases.h 2920 2017-06-29 09:58:07Z abarral $
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_t_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_t_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_t_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_t_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_t_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_t_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_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
395!attention: on impose dth
396        d_t_adv(l) = alpha*omega(l)/rcpd+                                  &
397     &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
398!        d_t_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,q,temp,u,v,play)
420       call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play)
421
422        do l=1,llm
423       d_t_adv(l) =  (dth_rico(l) +  dt_dyn(l))
424       d_q_adv(l,1) = (dqh_rico(l) +  dq_dyn(l,1))
425       d_q_adv(l,2) = 0.
426        enddo
427      endif  ! forcing_rico
428!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
429!---------------------------------------------------------------------
430! Interpolation forcing Arm_cu
431!---------------------------------------------------------------------
432      if (forcing_armcu) then
433
434        print*,                                                             &
435     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_armcu=',    &
436     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_armcu
437
438! time interpolation:
439! ATTENTION, cet appel ne convient pas pour TOGA !!
440! revoir 1DUTILS.h et les arguments
441      CALL interp_armcu_time(daytime,day1,annee_ref                         &
442     &            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu        &
443     &            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu          &
444     &            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof         &
445     &            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
446
447! vertical interpolation:
448! No vertical interpolation if nlev imposed to 19 or 40
449
450! For this case, fluxes are imposed
451       fsens=-1*sens_prof
452       flat=-1*flat_prof
453
454! Advective forcings are given in K or g/kg ... BY HOUR
455      do l = 1, llm
456       ug(l)= u_mod(l)
457       vg(l)= v_mod(l)
458       IF((phi(l)/RG).LT.1000) THEN
459         d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
460         d_q_adv(l,1) = adv_qt_prof/1000./3600.
461         d_q_adv(l,2) = 0.0
462!        print *,'INF1000: phi dth dq1 dq2',
463!    :  phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
464       ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN
465         fact=((phi(l)/RG)-1000.)/2000.
466         fact=1-fact
467         d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
468         d_q_adv(l,1) = adv_qt_prof*fact/1000./3600.
469         d_q_adv(l,2) = 0.0
470!        print *,'SUP1000: phi fact dth dq1 dq2',
471!    :  phi(l)/RG,fact,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
472       ELSE
473         d_t_adv(l) = 0.0
474         d_q_adv(l,1) = 0.0
475         d_q_adv(l,2) = 0.0
476!        print *,'SUP3000: phi dth dq1 dq2',
477!    :  phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
478       ENDIF
479      dt_cooling(l) = 0.0 
480!     print *,'Interp armcu: phi dth dq1 dq2',
481!    :  l,phi(l),d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
482      enddo
483      endif ! forcing_armcu
484!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485!---------------------------------------------------------------------
486! Interpolation forcing in time and onto model levels
487!---------------------------------------------------------------------
488      if (forcing_sandu) then
489
490        print*,                                                             &
491     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=',    &
492     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu
493
494! time interpolation:
495! ATTENTION, cet appel ne convient pas pour TOGA !!
496! revoir 1DUTILS.h et les arguments
497      CALL interp_sandu_time(daytime,day1,annee_ref                         &
498     &             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu       &
499     &             ,nlev_sandu                                              &
500     &             ,ts_sandu,ts_prof)
501
502        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
503
504! vertical interpolation:
505      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs                 &
506     &         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs           &
507     &         ,omega_profs,o3mmr_profs                                     &
508     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                       &
509     &         ,omega_mod,o3mmr_mod,mxcalc)
510!calcul de l'advection verticale
511!Calcul des gradients verticaux
512!initialisation
513      d_t_z(:)=0.
514      d_q_z(:)=0.
515      d_t_dyn_z(:)=0.
516      d_q_dyn_z(:)=0.
517! schema centre
518!     DO l=2,llm-1
519!      d_t_z(l)=(temp(l+1)-temp(l-1))
520!    &          /(play(l+1)-play(l-1))
521!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
522!    &          /(play(l+1)-play(l-1))
523! schema amont
524      DO l=2,llm-1
525       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
526       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
527!     print *,'l temp2 temp0 play2 play0 omega_mod',
528!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
529      ENDDO
530      d_t_z(1)=d_t_z(2)
531      d_q_z(1)=d_q_z(2)
532      d_t_z(llm)=d_t_z(llm-1)
533      d_q_z(llm)=d_q_z(llm-1)
534
535!  calcul de l advection verticale
536! Confusion w (m/s) et omega (Pa/s) !!
537      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
538      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
539!     do l=1,llm
540!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
541!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
542!     enddo
543
544
545! large-scale forcing : pour le cas Sandu ces forcages sont la SST
546! et une divergence constante -> profil de omega
547      tsurf = ts_prof
548      write(*,*) 'SST suivante: ',tsurf
549      do l = 1, llm
550       omega(l) = omega_mod(l)
551       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
552
553       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
554!
555!      d_t_adv(l) = 0.0
556!      d_q_adv(l,1) = 0.0
557!CR:test advection=0
558!calcul de l'advection verticale
559        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
560!        print*,'temp adv',l,-d_t_dyn_z(l)
561        d_q_adv(l,1) = -d_q_dyn_z(l)
562!        print*,'q adv',l,-d_q_dyn_z(l)
563       dt_cooling(l) = 0.0
564      enddo
565      endif ! forcing_sandu
566!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
567!---------------------------------------------------------------------
568! Interpolation forcing in time and onto model levels
569!---------------------------------------------------------------------
570      if (forcing_astex) then
571
572        print*,                                                             &
573     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=',    &
574     &    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex
575
576! time interpolation:
577! ATTENTION, cet appel ne convient pas pour TOGA !!
578! revoir 1DUTILS.h et les arguments
579      CALL interp_astex_time(daytime,day1,annee_ref                         &
580     &             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex       &
581     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex         &
582     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
583     &             ,ufa_prof,vfa_prof)
584
585        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
586
587! vertical interpolation:
588      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
589     &         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa                &
590     &         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa               &
591     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod        &
592     &         ,tke_mod,o3mmr_mod,mxcalc)
593!calcul de l'advection verticale
594!Calcul des gradients verticaux
595!initialisation
596      d_t_z(:)=0.
597      d_q_z(:)=0.
598      d_t_dyn_z(:)=0.
599      d_q_dyn_z(:)=0.
600! schema centre
601!     DO l=2,llm-1
602!      d_t_z(l)=(temp(l+1)-temp(l-1))
603!    &          /(play(l+1)-play(l-1))
604!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
605!    &          /(play(l+1)-play(l-1))
606! schema amont
607      DO l=2,llm-1
608       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
609       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
610!     print *,'l temp2 temp0 play2 play0 omega_mod',
611!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
612      ENDDO
613      d_t_z(1)=d_t_z(2)
614      d_q_z(1)=d_q_z(2)
615      d_t_z(llm)=d_t_z(llm-1)
616      d_q_z(llm)=d_q_z(llm-1)
617
618!  calcul de l advection verticale
619! Confusion w (m/s) et omega (Pa/s) !!
620      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
621      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
622!     do l=1,llm
623!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
624!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
625!     enddo
626
627
628! large-scale forcing : pour le cas Astex ces forcages sont la SST
629! la divergence,ug,vg,ufa,vfa
630      tsurf = ts_prof
631      write(*,*) 'SST suivante: ',tsurf
632      do l = 1, llm
633       omega(l) = w_mod(l)
634       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
635
636       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
637!
638!      d_t_adv(l) = 0.0
639!      d_q_adv(l,1) = 0.0
640!CR:test advection=0
641!calcul de l'advection verticale
642        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
643!        print*,'temp adv',l,-d_t_dyn_z(l)
644        d_q_adv(l,1) = -d_q_dyn_z(l)
645!        print*,'q adv',l,-d_q_dyn_z(l)
646       dt_cooling(l) = 0.0
647      enddo
648      endif ! forcing_astex
649
650!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
651!---------------------------------------------------------------------
652! Interpolation forcing standard case
653!---------------------------------------------------------------------
654      if (forcing_case) then
655
656        print*,                                                             &
657     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
658     &    daytime,day1,(daytime-day1)*86400.,                               &
659     &    (daytime-day1)*86400/pdt_cas
660
661! time interpolation:
662        CALL interp_case_time(daytime,day1,annee_ref                                        &
663!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
664     &       ,nt_cas,nlev_cas                                                               &
665     &       ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas,ug_cas,vg_cas                         &
666     &       ,vitw_cas,du_cas,hu_cas,vu_cas                                                 &
667     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
668     &       ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas                               &
669     &       ,uw_cas,vw_cas,q1_cas,q2_cas                                                   &
670     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
671     &       ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
672     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
673     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
674     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
675
676             ts_cur = ts_prof_cas
677             psurf=plev_prof_cas(1)
678
679! vertical interpolation:
680      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
681     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
682     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
683     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
684     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
685     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
686     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
687
688
689!calcul de l'advection verticale a partir du omega
690!Calcul des gradients verticaux
691!initialisation
692      d_t_z(:)=0.
693      d_q_z(:)=0.
694      d_u_z(:)=0.
695      d_v_z(:)=0.
696      d_t_dyn_z(:)=0.
697      d_q_dyn_z(:)=0.
698      d_u_dyn_z(:)=0.
699      d_v_dyn_z(:)=0.
700      DO l=2,llm-1
701       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
702       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
703       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
704       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
705      ENDDO
706      d_t_z(1)=d_t_z(2)
707      d_q_z(1)=d_q_z(2)
708      d_u_z(1)=d_u_z(2)
709      d_v_z(1)=d_v_z(2)
710      d_t_z(llm)=d_t_z(llm-1)
711      d_q_z(llm)=d_q_z(llm-1)
712      d_u_z(llm)=d_u_z(llm-1)
713      d_v_z(llm)=d_v_z(llm-1)
714
715!Calcul de l advection verticale
716
717      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
718
719      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
720      d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
721      d_v_dyn_z(:)=w_mod_cas(:)*d_v_z(:)
722
723!wind nudging
724      if (nudge_u.gt.0.) then
725        do l=1,llm
726           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
727        enddo
728      else
729        do l=1,llm
730        u(l) = u_mod_cas(l) 
731        enddo
732      endif
733
734      if (nudge_v.gt.0.) then
735        do l=1,llm
736           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
737        enddo
738      else
739        do l=1,llm
740        v(l) = v_mod_cas(l) 
741        enddo
742      endif
743
744      if (nudge_w.gt.0.) then
745        do l=1,llm
746           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
747        enddo
748      else
749        do l=1,llm
750        w(l) = w_mod_cas(l) 
751        enddo
752      endif
753
754!nudging of q and temp
755      if (nudge_t.gt.0.) then
756        do l=1,llm
757           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
758        enddo
759      endif
760      if (nudge_q.gt.0.) then
761        do l=1,llm
762           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
763        enddo
764      endif
765
766      do l = 1, llm
767       omega(l) = w_mod_cas(l)  ! juste car w_mod_cas en Pa/s (MPL 20170310)
768       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
769       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
770
771!calcul advection
772        if ((tend_u.eq.1).and.(tend_w.eq.0)) then
773           d_u_adv(l)=du_mod_cas(l)
774        else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
775           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
776        endif
777
778        if ((tend_v.eq.1).and.(tend_w.eq.0)) then
779           d_v_adv(l)=dv_mod_cas(l)
780        else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
781           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
782        endif
783
784        if ((tend_t.eq.1).and.(tend_w.eq.0)) then
785!           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
786           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
787        else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
788!           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
789           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
790        endif
791
792        if ((tend_q.eq.1).and.(tend_w.eq.0)) then
793!           d_q_adv(l,1)=dq_mod_cas(l)
794           d_q_adv(l,1)=-1*dq_mod_cas(l)
795        else if ((tend_q.eq.1).and.(tend_w.eq.1)) then
796!           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
797           d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
798        endif
799         
800        if (tend_rayo.eq.1) then
801           dt_cooling(l) = dtrad_mod_cas(l)
802!          print *,'dt_cooling=',dt_cooling(l)
803        else
804           dt_cooling(l) = 0.0
805        endif
806      enddo
807
808! Faut-il multiplier par -1 ? (MPL 20160713)
809      IF(ok_flux_surf) THEN
810       fsens=sens_prof_cas
811       flat=lat_prof_cas
812      ENDIF
813!
814      IF (ok_prescr_ust) THEN
815       ust=ustar_prof_cas
816       print *,'ust=',ust
817      ENDIF
818      endif ! forcing_case
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_u_z(:)=0.
880      d_v_z(:)=0.
881      d_t_dyn_z(:)=0.
882      d_th_dyn_z(:)=0.
883      d_q_dyn_z(:)=0.
884      d_u_dyn_z(:)=0.
885      d_v_dyn_z(:)=0.
886      DO l=2,llm-1
887       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
888       d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
889       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
890       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
891       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
892      ENDDO
893      d_t_z(1)=d_t_z(2)
894      d_th_z(1)=d_th_z(2)
895      d_q_z(1)=d_q_z(2)
896      d_u_z(1)=d_u_z(2)
897      d_v_z(1)=d_v_z(2)
898      d_t_z(llm)=d_t_z(llm-1)
899      d_th_z(llm)=d_th_z(llm-1)
900      d_q_z(llm)=d_q_z(llm-1)
901      d_u_z(llm)=d_u_z(llm-1)
902      d_v_z(llm)=d_v_z(llm-1)
903
904!Calcul de l advection verticale
905! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
906      d_t_dyn_z(:)=omega_mod_cas(:)*d_t_z(:)
907      d_th_dyn_z(:)=omega_mod_cas(:)*d_th_z(:)
908      d_q_dyn_z(:)=omega_mod_cas(:)*d_q_z(:)
909      d_u_dyn_z(:)=omega_mod_cas(:)*d_u_z(:)
910      d_v_dyn_z(:)=omega_mod_cas(:)*d_v_z(:)
911
912!geostrophic wind
913      if (forc_geo.eq.1) then
914        do l=1,llm
915        ug(l) = ug_mod_cas(l) 
916        vg(l) = vg_mod_cas(l) 
917        enddo
918      endif
919!wind nudging
920      if (nudging_u.gt.0.) then
921        do l=1,llm
922           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
923        enddo
924!     else
925!       do l=1,llm
926!          u(l) = u_mod_cas(l) 
927!       enddo
928      endif
929
930      if (nudging_v.gt.0.) then
931        do l=1,llm
932           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
933        enddo
934!     else
935!       do l=1,llm
936!          v(l) = v_mod_cas(l) 
937!       enddo
938      endif
939
940      if (nudging_w.gt.0.) then
941        do l=1,llm
942           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
943        enddo
944 !    else
945 !      do l=1,llm
946 !         w(l) = w_mod_cas(l) 
947 !      enddo
948      endif
949
950!nudging of q and temp
951      if (nudging_t.gt.0.) then
952        do l=1,llm
953           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
954        enddo
955      endif
956      if (nudging_q.gt.0.) then
957        do l=1,llm
958           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
959        enddo
960      endif
961
962      do l = 1, llm
963! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
964       omega(l) = omega_mod_cas(l)
965       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
966       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
967
968!calcul advections
969        if ((forc_u.eq.1).and.(forc_w.eq.0)) then
970           d_u_adv(l)=du_mod_cas(l)
971        else if ((forc_u.eq.1).and.(forc_w.eq.1)) then
972           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
973        endif
974
975        if ((forc_v.eq.1).and.(forc_w.eq.0)) then
976           d_v_adv(l)=dv_mod_cas(l)
977        else if ((forc_v.eq.1).and.(forc_w.eq.1)) then
978           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
979        endif
980
981! Puisque dth a ete converti en dt, on traite de la meme facon
982! les flags tadv et thadv
983        if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.0)) then
984!          d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
985           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
986        else if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.1)) then
987!          d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
988           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
989        endif
990
991!       if ((thadv.eq.1) .and. (forc_w.eq.0)) then
992!          d_t_adv(l)=alpha*omega(l)/rcpd-dth_mod_cas(l)
993!          d_t_adv(l)=alpha*omega(l)/rcpd+dth_mod_cas(l)
994!       else if ((thadv.eq.1) .and. (forc_w.eq.1)) then
995!          d_t_adv(l)=alpha*omega(l)/rcpd-hth_mod_cas(l)-d_t_dyn_z(l)
996!          d_t_adv(l)=alpha*omega(l)/rcpd+hth_mod_cas(l)-d_t_dyn_z(l)
997!       endif
998
999        if ((qadv.eq.1) .and. (forc_w.eq.0)) then
1000           d_q_adv(l,1)=dq_mod_cas(l)
1001!          d_q_adv(l,1)=-1*dq_mod_cas(l)
1002        else if ((qadv.eq.1) .and. (forc_w.eq.1)) then
1003           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
1004!          d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
1005        endif
1006         
1007        if (trad.eq.1) then
1008           tend_rayo=1
1009           dt_cooling(l) = dtrad_mod_cas(l)
1010!          print *,'dt_cooling=',dt_cooling(l)
1011        else
1012           dt_cooling(l) = 0.0
1013        endif
1014      enddo
1015
1016! Faut-il multiplier par -1 ? (MPL 20160713)
1017      IF(ok_flux_surf) THEN
1018       fsens=-1.*sens_prof_cas
1019       flat=-1.*lat_prof_cas
1020       print *,'1D_interp: sens,flat',fsens,flat
1021      ENDIF
1022!
1023      IF (ok_prescr_ust) THEN
1024       ust=ustar_prof_cas
1025       print *,'ust=',ust
1026      ENDIF
1027      endif ! forcing_case2
1028!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1029
Note: See TracBrowser for help on using the repository browser.