source: LMDZ6/trunk/libf/phylmd/dyn1d/old_1D_interp_cases.h @ 3541

Last change on this file since 3541 was 3541, checked in by fhourdin, 5 years ago

Gros nettoyage en cours sur le 1D.
Le nouveau lmdz1d.F90 est une coquille vide qui choisit entre
old_lmdz1d.F90 (l'ancien lmdz1d.F90) et scm.F90 (le nouveau au format standard).
Plusieur fichiers sont renommés old_truc, le truc étant au format standard,
nettoyé des anciens formats.
Le 1DUTILS.h est lui même séparé entre des routines génériques venant remplacer
notamment des routines de dyn3d (la vocation d'origine de 1DUTILS.h) et
les routiles de lecture spécifiques mises dans old_1DUTILS.h
On perdra un peu de l'utilister de svn au moment de cette grosse bascule.
Mais les old_ sont faits pour ne plus bouger, et les versions standard
sont en pleine évolution.
Fredho

File size: 39.4 KB
Line 
1!
2! $Id: 1D_interp_cases.h 3537 2019-06-19 08:29:16Z 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_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*,'FORCING CASE forcing_case'
657
658        print*,                                                             &
659     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
660     &    daytime,day1,(daytime-day1)*86400.,                               &
661     &    (daytime-day1)*86400/pdt_cas
662
663! time interpolation:
664        CALL interp_case_time(daytime,day1,annee_ref                                        &
665!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
666     &       ,nt_cas,nlev_cas                                                               &
667     &       ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas,ug_cas,vg_cas                         &
668     &       ,vitw_cas,du_cas,hu_cas,vu_cas                                                 &
669     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
670     &       ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas                               &
671     &       ,uw_cas,vw_cas,q1_cas,q2_cas                                                   &
672     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
673     &       ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
674     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
675     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
676     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
677
678             ts_cur = ts_prof_cas
679             psurf=plev_prof_cas(1)
680
681! vertical interpolation:
682      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
683     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
684     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
685     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
686     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
687     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
688     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
689
690
691!calcul de l'advection verticale a partir du omega
692!Calcul des gradients verticaux
693!initialisation
694      d_t_z(:)=0.
695      d_q_z(:)=0.
696      d_u_z(:)=0.
697      d_v_z(:)=0.
698      d_t_dyn_z(:)=0.
699      d_q_dyn_z(:)=0.
700      d_u_dyn_z(:)=0.
701      d_v_dyn_z(:)=0.
702      DO l=2,llm-1
703       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
704       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
705       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
706       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
707      ENDDO
708      d_t_z(1)=d_t_z(2)
709      d_q_z(1)=d_q_z(2)
710      d_u_z(1)=d_u_z(2)
711      d_v_z(1)=d_v_z(2)
712      d_t_z(llm)=d_t_z(llm-1)
713      d_q_z(llm)=d_q_z(llm-1)
714      d_u_z(llm)=d_u_z(llm-1)
715      d_v_z(llm)=d_v_z(llm-1)
716
717!Calcul de l advection verticale
718
719      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
720
721      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
722      d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
723      d_v_dyn_z(:)=w_mod_cas(:)*d_v_z(:)
724
725!wind nudging
726      if (nudge_u.gt.0.) then
727        do l=1,llm
728           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
729        enddo
730      else
731        do l=1,llm
732        u(l) = u_mod_cas(l) 
733        enddo
734      endif
735
736      if (nudge_v.gt.0.) then
737        do l=1,llm
738           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
739        enddo
740      else
741        do l=1,llm
742        v(l) = v_mod_cas(l) 
743        enddo
744      endif
745
746      if (nudge_w.gt.0.) then
747        do l=1,llm
748           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
749        enddo
750      else
751        do l=1,llm
752        w(l) = w_mod_cas(l) 
753        enddo
754      endif
755
756!nudging of q and temp
757      if (nudge_t.gt.0.) then
758        do l=1,llm
759           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
760        enddo
761      endif
762      if (nudge_q.gt.0.) then
763        do l=1,llm
764           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
765        enddo
766      endif
767
768      do l = 1, llm
769       omega(l) = w_mod_cas(l)  ! juste car w_mod_cas en Pa/s (MPL 20170310)
770       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
771       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
772
773!calcul advection
774        if ((tend_u.eq.1).and.(tend_w.eq.0)) then
775           d_u_adv(l)=du_mod_cas(l)
776        else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
777           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
778        endif
779
780        if ((tend_v.eq.1).and.(tend_w.eq.0)) then
781           d_v_adv(l)=dv_mod_cas(l)
782        else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
783           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
784        endif
785
786        if ((tend_t.eq.1).and.(tend_w.eq.0)) then
787!           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
788           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
789        else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
790!           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
791           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
792        endif
793
794        if ((tend_q.eq.1).and.(tend_w.eq.0)) then
795!           d_q_adv(l,1)=dq_mod_cas(l)
796           d_q_adv(l,1)=-1*dq_mod_cas(l)
797        else if ((tend_q.eq.1).and.(tend_w.eq.1)) then
798!           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
799           d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
800        endif
801         
802        if (tend_rayo.eq.1) then
803           dt_cooling(l) = dtrad_mod_cas(l)
804!          print *,'dt_cooling=',dt_cooling(l)
805        else
806           dt_cooling(l) = 0.0
807        endif
808      enddo
809
810! Faut-il multiplier par -1 ? (MPL 20160713)
811      IF(ok_flux_surf) THEN
812       fsens=sens_prof_cas
813       flat=lat_prof_cas
814      ENDIF
815!
816      IF (ok_prescr_ust) THEN
817       ust=ustar_prof_cas
818       print *,'ust=',ust
819      ENDIF
820      endif ! forcing_case
821
822!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
823!---------------------------------------------------------------------
824! Interpolation forcing standard case
825!---------------------------------------------------------------------
826      if (forcing_case2 .OR. forcing_SCM) then
827
828         print*,'FORCING CASE forcing_case2'
829        print*,                                                             &
830     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
831     &    daytime,day1,(daytime-day1)*86400.,                               &
832     &    (daytime-day1)*86400/pdt_cas
833
834! time interpolation:
835        CALL interp2_case_time(daytime,day1,annee_ref                                       &
836!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
837     &       ,nt_cas,nlev_cas                                                               &
838     &       ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas      &
839     &       ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas             &
840     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
841     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
842     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
843!
844     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
845     &       ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas                              &
846     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas    &
847     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
848     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
849     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
850     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
851     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
852
853             ts_cur = ts_prof_cas
854!            psurf=plev_prof_cas(1)
855             psurf=ps_prof_cas
856
857! vertical interpolation:
858      CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas                                              &
859     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
860     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
861     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                      &
862     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
863     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
864     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                    &
865!
866     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas          &
867     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                         &
868     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
869     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
870     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
871
872
873      DO l=1,llm
874      teta(l)=temp(l)*(100000./play(l))**(rd/rcpd)
875      ENDDO
876!calcul de l'advection verticale a partir du omega
877!Calcul des gradients verticaux
878!initialisation
879      d_t_z(:)=0.
880      d_th_z(:)=0.
881      d_q_z(:)=0.
882      d_u_z(:)=0.
883      d_v_z(:)=0.
884      d_t_dyn_z(:)=0.
885      d_th_dyn_z(:)=0.
886      d_q_dyn_z(:)=0.
887      d_u_dyn_z(:)=0.
888      d_v_dyn_z(:)=0.
889      DO l=2,llm-1
890       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
891       d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
892       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
893       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
894       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
895      ENDDO
896      d_t_z(1)=d_t_z(2)
897      d_th_z(1)=d_th_z(2)
898      d_q_z(1)=d_q_z(2)
899      d_u_z(1)=d_u_z(2)
900      d_v_z(1)=d_v_z(2)
901      d_t_z(llm)=d_t_z(llm-1)
902      d_th_z(llm)=d_th_z(llm-1)
903      d_q_z(llm)=d_q_z(llm-1)
904      d_u_z(llm)=d_u_z(llm-1)
905      d_v_z(llm)=d_v_z(llm-1)
906
907!Calcul de l advection verticale
908! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
909      d_t_dyn_z(:)=omega_mod_cas(:)*d_t_z(:)
910      d_th_dyn_z(:)=omega_mod_cas(:)*d_th_z(:)
911      d_q_dyn_z(:)=omega_mod_cas(:)*d_q_z(:)
912      d_u_dyn_z(:)=omega_mod_cas(:)*d_u_z(:)
913      d_v_dyn_z(:)=omega_mod_cas(:)*d_v_z(:)
914
915!geostrophic wind
916      if (forc_geo.eq.1) then
917        do l=1,llm
918        ug(l) = ug_mod_cas(l) 
919        vg(l) = vg_mod_cas(l) 
920        enddo
921      endif
922!wind nudging
923      if (nudging_u.gt.0.) then
924        do l=1,llm
925           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
926        enddo
927!     else
928!       do l=1,llm
929!          u(l) = u_mod_cas(l) 
930!       enddo
931      endif
932
933      if (nudging_v.gt.0.) then
934        do l=1,llm
935           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
936        enddo
937!     else
938!       do l=1,llm
939!          v(l) = v_mod_cas(l) 
940!       enddo
941      endif
942
943      if (nudging_w.gt.0.) then
944        do l=1,llm
945           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
946        enddo
947 !    else
948 !      do l=1,llm
949 !         w(l) = w_mod_cas(l) 
950 !      enddo
951      endif
952
953!nudging of q and temp
954      if (nudging_t.gt.0.) then
955        do l=1,llm
956           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
957        enddo
958      endif
959      if (nudging_q.gt.0.) then
960        do l=1,llm
961           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
962        enddo
963      endif
964
965      do l = 1, llm
966! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
967       omega(l) = omega_mod_cas(l)
968       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
969       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
970
971!calcul advections
972        if ((forc_u.eq.1).and.(forc_w.eq.0)) then
973           d_u_adv(l)=du_mod_cas(l)
974        else if ((forc_u.eq.1).and.(forc_w.eq.1)) then
975           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
976        endif
977
978        if ((forc_v.eq.1).and.(forc_w.eq.0)) then
979           d_v_adv(l)=dv_mod_cas(l)
980        else if ((forc_v.eq.1).and.(forc_w.eq.1)) then
981           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
982        endif
983
984! Puisque dth a ete converti en dt, on traite de la meme facon
985! les flags tadv et thadv
986        if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.0)) then
987!          d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
988           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
989        else if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.1)) then
990!          d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
991           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
992        endif
993
994!       if ((thadv.eq.1) .and. (forc_w.eq.0)) then
995!          d_t_adv(l)=alpha*omega(l)/rcpd-dth_mod_cas(l)
996!          d_t_adv(l)=alpha*omega(l)/rcpd+dth_mod_cas(l)
997!       else if ((thadv.eq.1) .and. (forc_w.eq.1)) then
998!          d_t_adv(l)=alpha*omega(l)/rcpd-hth_mod_cas(l)-d_t_dyn_z(l)
999!          d_t_adv(l)=alpha*omega(l)/rcpd+hth_mod_cas(l)-d_t_dyn_z(l)
1000!       endif
1001
1002        if ((qadv.eq.1) .and. (forc_w.eq.0)) then
1003           d_q_adv(l,1)=dq_mod_cas(l)
1004!          d_q_adv(l,1)=-1*dq_mod_cas(l)
1005        else if ((qadv.eq.1) .and. (forc_w.eq.1)) then
1006           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
1007!          d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
1008        endif
1009         
1010        if (trad.eq.1) then
1011           tend_rayo=1
1012           dt_cooling(l) = dtrad_mod_cas(l)
1013!          print *,'dt_cooling=',dt_cooling(l)
1014        else
1015           dt_cooling(l) = 0.0
1016        endif
1017      enddo
1018
1019! Faut-il multiplier par -1 ? (MPL 20160713)
1020      IF(ok_flux_surf) THEN
1021       fsens=-1.*sens_prof_cas
1022       flat=-1.*lat_prof_cas
1023       print *,'1D_interp: sens,flat',fsens,flat
1024      ENDIF
1025!
1026      IF (ok_prescr_ust) THEN
1027       ust=ustar_prof_cas
1028       print *,'ust=',ust
1029      ENDIF
1030      endif ! forcing_case2
1031!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1032
Note: See TracBrowser for help on using the repository browser.