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

Last change on this file since 5450 was 3780, checked in by evignon, 4 years ago

Premiere comission Etienne: changements pour le 1D (forcage en Ts au dessus des continents) et inclusion drag arbres dans yamada4_num=6

File size: 39.5 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! EV: tg instead of ts_cur
65        if (type_ts_forcing.eq.1) tg = ts_prof ! 
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! EV tg instead of ts_cur
116        if (type_ts_forcing.eq.1) tg = tg_prof ! SST used
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!EV tg instead of ts_cur
209        if (type_ts_forcing.eq.1) tg = tg_prof ! SST used
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! EV tg instead of ts_cur
502        if (type_ts_forcing.eq.1) tg = 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! EV tg instead of ts_cur
585        if (type_ts_forcing.eq.1) tg = ts_prof ! SST used
586! vertical interpolation:
587      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
588     &         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa                &
589     &         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa               &
590     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod        &
591     &         ,tke_mod,o3mmr_mod,mxcalc)
592!calcul de l'advection verticale
593!Calcul des gradients verticaux
594!initialisation
595      d_t_z(:)=0.
596      d_q_z(:)=0.
597      d_t_dyn_z(:)=0.
598      d_q_dyn_z(:)=0.
599! schema centre
600!     DO l=2,llm-1
601!      d_t_z(l)=(temp(l+1)-temp(l-1))
602!    &          /(play(l+1)-play(l-1))
603!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
604!    &          /(play(l+1)-play(l-1))
605! schema amont
606      DO l=2,llm-1
607       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
608       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
609!     print *,'l temp2 temp0 play2 play0 omega_mod',
610!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
611      ENDDO
612      d_t_z(1)=d_t_z(2)
613      d_q_z(1)=d_q_z(2)
614      d_t_z(llm)=d_t_z(llm-1)
615      d_q_z(llm)=d_q_z(llm-1)
616
617!  calcul de l advection verticale
618! Confusion w (m/s) et omega (Pa/s) !!
619      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
620      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
621!     do l=1,llm
622!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
623!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
624!     enddo
625
626
627! large-scale forcing : pour le cas Astex ces forcages sont la SST
628! la divergence,ug,vg,ufa,vfa
629      tsurf = ts_prof
630      write(*,*) 'SST suivante: ',tsurf
631      do l = 1, llm
632       omega(l) = w_mod(l)
633       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
634
635       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
636!
637!      d_t_adv(l) = 0.0
638!      d_q_adv(l,1) = 0.0
639!CR:test advection=0
640!calcul de l'advection verticale
641        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
642!        print*,'temp adv',l,-d_t_dyn_z(l)
643        d_q_adv(l,1) = -d_q_dyn_z(l)
644!        print*,'q adv',l,-d_q_dyn_z(l)
645       dt_cooling(l) = 0.0
646      enddo
647      endif ! forcing_astex
648
649!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
650!---------------------------------------------------------------------
651! Interpolation forcing standard case
652!---------------------------------------------------------------------
653      if (forcing_case) then
654
655         print*,'FORCING CASE forcing_case'
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! EV tg instead of ts_cur
677
678             tg = 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! EV tg instead of ts_cur
853
854             tg = ts_prof_cas
855!            psurf=plev_prof_cas(1)
856             psurf=ps_prof_cas
857
858! vertical interpolation:
859      CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas                                              &
860     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
861     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
862     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                      &
863     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
864     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
865     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                    &
866!
867     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas          &
868     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                         &
869     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
870     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
871     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
872
873
874      DO l=1,llm
875      teta(l)=temp(l)*(100000./play(l))**(rd/rcpd)
876      ENDDO
877!calcul de l'advection verticale a partir du omega
878!Calcul des gradients verticaux
879!initialisation
880      d_t_z(:)=0.
881      d_th_z(:)=0.
882      d_q_z(:)=0.
883      d_u_z(:)=0.
884      d_v_z(:)=0.
885      d_t_dyn_z(:)=0.
886      d_th_dyn_z(:)=0.
887      d_q_dyn_z(:)=0.
888      d_u_dyn_z(:)=0.
889      d_v_dyn_z(:)=0.
890      DO l=2,llm-1
891       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
892       d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
893       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
894       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
895       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
896      ENDDO
897      d_t_z(1)=d_t_z(2)
898      d_th_z(1)=d_th_z(2)
899      d_q_z(1)=d_q_z(2)
900      d_u_z(1)=d_u_z(2)
901      d_v_z(1)=d_v_z(2)
902      d_t_z(llm)=d_t_z(llm-1)
903      d_th_z(llm)=d_th_z(llm-1)
904      d_q_z(llm)=d_q_z(llm-1)
905      d_u_z(llm)=d_u_z(llm-1)
906      d_v_z(llm)=d_v_z(llm-1)
907
908!Calcul de l advection verticale
909! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
910      d_t_dyn_z(:)=omega_mod_cas(:)*d_t_z(:)
911      d_th_dyn_z(:)=omega_mod_cas(:)*d_th_z(:)
912      d_q_dyn_z(:)=omega_mod_cas(:)*d_q_z(:)
913      d_u_dyn_z(:)=omega_mod_cas(:)*d_u_z(:)
914      d_v_dyn_z(:)=omega_mod_cas(:)*d_v_z(:)
915
916!geostrophic wind
917      if (forc_geo.eq.1) then
918        do l=1,llm
919        ug(l) = ug_mod_cas(l) 
920        vg(l) = vg_mod_cas(l) 
921        enddo
922      endif
923!wind nudging
924      if (nudging_u.gt.0.) then
925        do l=1,llm
926           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
927        enddo
928!     else
929!       do l=1,llm
930!          u(l) = u_mod_cas(l) 
931!       enddo
932      endif
933
934      if (nudging_v.gt.0.) then
935        do l=1,llm
936           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
937        enddo
938!     else
939!       do l=1,llm
940!          v(l) = v_mod_cas(l) 
941!       enddo
942      endif
943
944      if (nudging_w.gt.0.) then
945        do l=1,llm
946           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
947        enddo
948 !    else
949 !      do l=1,llm
950 !         w(l) = w_mod_cas(l) 
951 !      enddo
952      endif
953
954!nudging of q and temp
955      if (nudging_t.gt.0.) then
956        do l=1,llm
957           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
958        enddo
959      endif
960      if (nudging_qv.gt.0.) then
961        do l=1,llm
962           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
963        enddo
964      endif
965
966      do l = 1, llm
967! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
968       omega(l) = omega_mod_cas(l)
969       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
970       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
971
972!calcul advections
973        if ((forc_u.eq.1).and.(forc_w.eq.0)) then
974           d_u_adv(l)=du_mod_cas(l)
975        else if ((forc_u.eq.1).and.(forc_w.eq.1)) then
976           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
977        endif
978
979        if ((forc_v.eq.1).and.(forc_w.eq.0)) then
980           d_v_adv(l)=dv_mod_cas(l)
981        else if ((forc_v.eq.1).and.(forc_w.eq.1)) then
982           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
983        endif
984
985! Puisque dth a ete converti en dt, on traite de la meme facon
986! les flags tadv et thadv
987        if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.0)) then
988!          d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
989           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
990        else if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.1)) then
991!          d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
992           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
993        endif
994
995!       if ((thadv.eq.1) .and. (forc_w.eq.0)) then
996!          d_t_adv(l)=alpha*omega(l)/rcpd-dth_mod_cas(l)
997!          d_t_adv(l)=alpha*omega(l)/rcpd+dth_mod_cas(l)
998!       else if ((thadv.eq.1) .and. (forc_w.eq.1)) then
999!          d_t_adv(l)=alpha*omega(l)/rcpd-hth_mod_cas(l)-d_t_dyn_z(l)
1000!          d_t_adv(l)=alpha*omega(l)/rcpd+hth_mod_cas(l)-d_t_dyn_z(l)
1001!       endif
1002
1003        if ((qadv.eq.1) .and. (forc_w.eq.0)) then
1004           d_q_adv(l,1)=dq_mod_cas(l)
1005!          d_q_adv(l,1)=-1*dq_mod_cas(l)
1006        else if ((qadv.eq.1) .and. (forc_w.eq.1)) then
1007           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
1008!          d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
1009        endif
1010         
1011        if (trad.eq.1) then
1012           tend_rayo=1
1013           dt_cooling(l) = dtrad_mod_cas(l)
1014!          print *,'dt_cooling=',dt_cooling(l)
1015        else
1016           dt_cooling(l) = 0.0
1017        endif
1018      enddo
1019
1020! Faut-il multiplier par -1 ? (MPL 20160713)
1021      IF(ok_flux_surf) THEN
1022       fsens=-1.*sens_prof_cas
1023       flat=-1.*lat_prof_cas
1024       print *,'1D_interp: sens,flat',fsens,flat
1025      ENDIF
1026!
1027      IF (ok_prescr_ust) THEN
1028       ust=ustar_prof_cas
1029       print *,'ust=',ust
1030      ENDIF
1031      endif ! forcing_case2
1032!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1033
Note: See TracBrowser for help on using the repository browser.