source: LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h @ 5437

Last change on this file since 5437 was 2920, checked in by fhourdin, 7 years ago

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

  • Property svn:keywords set to Id
File size: 42.0 KB
Line 
1!
2! $Id: 1D_read_forc_cases.h 2920 2017-06-29 09:58:07Z evignon $
3!
4!----------------------------------------------------------------------
5! forcing_les = .T. : Impose a constant cooling
6! forcing_radconv = .T. : Pure radiative-convective equilibrium:
7!----------------------------------------------------------------------
8
9
10      nq1=0
11      nq2=0
12
13      if (forcing_les .or. forcing_radconv                                      &
14     &    .or. forcing_GCSSold .or. forcing_fire) then
15
16      if (forcing_fire) then
17!----------------------------------------------------------------------
18!read fire forcings from fire.nc
19!----------------------------------------------------------------------
20      fich_fire='fire.nc'
21      call read_fire(fich_fire,nlev_fire,nt_fire                                &
22     &     ,height,tttprof,qtprof,uprof,vprof,e12prof                           &
23     &     ,ugprof,vgprof,wfls,dqtdxls                                          &
24     &     ,dqtdyls,dqtdtls,thlpcar)
25      write(*,*) 'Forcing FIRE lu'
26      kmax=120            ! nombre de niveaux dans les profils et forcages
27      else 
28!----------------------------------------------------------------------
29! Read profiles from files: prof.inp.001 and lscale.inp.001
30! (repris de readlesfiles)
31!----------------------------------------------------------------------
32
33      call readprofiles(nlev_max,kmax,nqtot,height,                             &
34     &           tttprof,qtprof,uprof,vprof,                                    &
35     &           e12prof,ugprof,vgprof,                                         &
36     &           wfls,dqtdxls,dqtdyls,dqtdtls,                                  &
37     &           thlpcar,qprof,nq1,nq2)
38      endif
39
40! compute altitudes of play levels.
41      zlay(1) =zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
42      do l = 2,llm
43        zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf)
44      enddo
45
46!----------------------------------------------------------------------
47! Interpolation of the profiles given on the input file to
48! model levels
49!----------------------------------------------------------------------
50      zlay(1) = zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
51      do l=1,llm
52        ! Above the max altutide of the input file
53
54        if (zlay(l)<height(kmax)) mxcalc=l
55
56        frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
57        ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
58       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
59          temp(l) = ttt*(play(l)/pzero)**rkappa
60          teta(l) = ttt
61       else
62          temp(l) = ttt
63          teta(l) = ttt*(pzero/play(l))**rkappa
64       endif
65          print *,' temp,teta ',l,temp(l),teta(l)
66        q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
67        u(l)    =  uprof(kmax)-frac*(  uprof(kmax)-  uprof(kmax-1))
68        v(l)    =  vprof(kmax)-frac*(  vprof(kmax)-  vprof(kmax-1))
69        ug(l)   = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1))
70        vg(l)   = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1))
71        IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2)                         &
72     &               -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2))
73        omega(l)=   wfls(kmax)-frac*(   wfls(kmax)-   wfls(kmax-1))
74
75        dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1))
76        dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
77        do k=2,kmax
78          print *,'k l height(k) height(k-1) zlay(l) frac=',k,l,height(k),height(k-1),zlay(l),frac
79          frac = (height(k)-zlay(l))/(height(k)-height(k-1))
80          if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
81          if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
82            ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
83       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
84          temp(l) = ttt*(play(l)/pzero)**rkappa
85          teta(l) = ttt
86       else
87          temp(l) = ttt
88          teta(l) = ttt*(pzero/play(l))**rkappa
89       endif
90          print *,' temp,teta ',l,temp(l),teta(l)
91            q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
92            u(l)    =  uprof(k)-frac*(  uprof(k)-  uprof(k-1))
93            v(l)    =  vprof(k)-frac*(  vprof(k)-  vprof(k-1))
94            ug(l)   = ugprof(k)-frac*( ugprof(k)- ugprof(k-1))
95            vg(l)   = vgprof(k)-frac*( vgprof(k)- vgprof(k-1))
96            IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2)                        &
97     &                   -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2))
98            omega(l)=   wfls(k)-frac*(   wfls(k)-   wfls(k-1))
99            dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1))
100            dt_cooling(l)=thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
101          elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
102            ttt =tttprof(1)
103       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
104          temp(l) = ttt*(play(l)/pzero)**rkappa
105          teta(l) = ttt
106       else
107          temp(l) = ttt
108          teta(l) = ttt*(pzero/play(l))**rkappa
109       endif
110            q(l,1)  = qtprof(1)
111            u(l)    =  uprof(1)
112            v(l)    =  vprof(1)
113            ug(l)   = ugprof(1)
114            vg(l)   = vgprof(1)
115            omega(l)=   wfls(1)
116            IF (nq2>0) q(l,nq1:nq2)=qprof(1,nq1:nq2)
117            dq_dyn(l,1)  =dqtdtls(1)
118            dt_cooling(l)=thlpcar(1)
119          endif
120        enddo
121
122        temp(l)=max(min(temp(l),350.),150.)
123        rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
124        if (l .lt. llm) then
125          zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l))
126        endif
127        omega2(l)=-rho(l)*omega(l)
128        omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s
129        if (l>1) then
130        if(zlay(l-1)>height(kmax)) then
131           omega(l)=0.0
132           omega2(l)=0.0
133        endif   
134        endif
135        if(q(l,1)<0.) q(l,1)=0.0
136        q(l,2)  = 0.0
137      enddo
138
139      endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire
140!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141!---------------------------------------------------------------------
142! Forcing for GCSSold:
143!---------------------------------------------------------------------
144      if (forcing_GCSSold) then
145       fich_gcssold_ctl = './forcing.ctl'
146       fich_gcssold_dat = './forcing8.dat'
147       call copie(llm,play,psurf,fich_gcssold_ctl)
148       call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,         &
149     &               ht_gcssold,hq_gcssold,hw_gcssold,                      &
150     &               hu_gcssold,hv_gcssold,                                 &
151     &               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,             &
152     &               imp_fcg_gcssold,ts_fcg_gcssold,                        &
153     &               Tp_fcg_gcssold,Turb_fcg_gcssold)
154       print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold
155      endif ! forcing_GCSSold
156!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
157!---------------------------------------------------------------------
158! Forcing for RICO:
159!---------------------------------------------------------------------
160      if (forcing_rico) then
161
162!       call writefield_phy('omega', omega,llm+1)
163      fich_rico = 'rico.txt'
164       call read_rico(fich_rico,nlev_rico,ps_rico,play                      &
165     &             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico              &
166     &             ,dth_rico,dqh_rico)
167        print*, ' on a lu et prepare RICO'
168
169       mxcalc=llm
170       print *, airefi, ' airefi '
171       do l = 1, llm
172       rho(l)  = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l)))
173       temp(l) = t_rico(l)
174       q(l,1) = q_rico(l)
175       q(l,2) = 0.0
176       u(l) = u_rico(l)
177       v(l) = v_rico(l)
178       ug(l)=u_rico(l)
179       vg(l)=v_rico(l)
180       omega(l) = -w_rico(l)*rg
181       omega2(l) = omega(l)/rg*airefi
182       enddo
183      endif
184!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185!---------------------------------------------------------------------
186! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) :
187!---------------------------------------------------------------------
188
189      if (forcing_toga) then
190
191! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps):
192      fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt'
193      CALL read_togacoare(fich_toga,nlev_toga,nt_toga                       &
194     &         ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga        &
195     &         ,ht_toga,vt_toga,hq_toga,vq_toga)
196
197       write(*,*) 'Forcing TOGA lu'
198
199! time interpolation for initial conditions:
200      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
201      CALL interp_toga_time(daytime,day1,annee_ref                          &
202     &             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga           &
203     &             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga        &
204     &             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga           &
205     &             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof    &
206     &             ,ht_prof,vt_prof,hq_prof,vq_prof)
207
208! vertical interpolation:
209      CALL interp_toga_vertical(play,nlev_toga,plev_prof                    &
210     &         ,t_prof,q_prof,u_prof,v_prof,w_prof                          &
211     &         ,ht_prof,vt_prof,hq_prof,vq_prof                             &
212     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
213     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
214       write(*,*) 'Profil initial forcing TOGA interpole'
215
216! initial and boundary conditions :
217      tsurf = ts_prof
218      write(*,*) 'SST initiale: ',tsurf
219      do l = 1, llm
220       temp(l) = t_mod(l)
221       q(l,1) = q_mod(l)
222       q(l,2) = 0.0
223       u(l) = u_mod(l)
224       v(l) = v_mod(l)
225       omega(l) = w_mod(l)
226       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
227!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
228!?       omega2(l)=-rho(l)*omega(l)
229       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
230       d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
231       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
232       d_q_adv(l,2) = 0.0
233      enddo
234
235      endif ! forcing_toga
236!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
237!---------------------------------------------------------------------
238! Forcing from TWPICE experiment (Shaocheng et al. 2010) :
239!---------------------------------------------------------------------
240
241      if (forcing_twpice) then
242!read TWP-ICE forcings
243     fich_twpice='d_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
244      call read_twpice(fich_twpice,nlev_twpi,nt_twpi                        &
245     &     ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi            &
246     &     ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)
247
248      write(*,*) 'Forcing TWP-ICE lu'
249!Time interpolation for initial conditions using TOGA interpolation routine
250         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
251      CALL interp_toga_time(daytime,day1,annee_ref                          &
252     &          ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi    &
253     &             ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi    &
254     &             ,ht_twpi,vt_twpi,hq_twpi,vq_twpi                         &
255     &             ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp             &
256     &             ,u_proftwp,v_proftwp,w_proftwp                           &
257     &             ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
258
259! vertical interpolation using TOGA interpolation routine:
260!      write(*,*)'avant interp vert', t_proftwp
261      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp                 &
262     &         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp           &
263     &         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp                 &
264     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
265     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
266!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
267
268! initial and boundary conditions :
269!      tsurf = ts_proftwp
270      write(*,*) 'SST initiale: ',tsurf
271      do l = 1, llm
272       temp(l) = t_mod(l)
273       q(l,1) = q_mod(l)
274       q(l,2) = 0.0
275       u(l) = u_mod(l)
276       v(l) = v_mod(l)
277       omega(l) = w_mod(l)
278       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
279
280       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
281!on applique le forcage total au premier pas de temps
282!attention: signe different de toga
283       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
284       d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
285       d_q_adv(l,2) = 0.0
286      enddo     
287       
288      endif !forcing_twpice
289
290!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291!---------------------------------------------------------------------
292! Forcing from AMMA experiment (Couvreux et al. 2010) :
293!---------------------------------------------------------------------
294
295      if (forcing_amma) then
296
297      call read_1D_cases
298
299      write(*,*) 'Forcing AMMA lu'
300
301!champs initiaux:
302      do k=1,nlev_amma
303         th_ammai(k)=th_amma(k)
304         q_ammai(k)=q_amma(k)
305         u_ammai(k)=u_amma(k)
306         v_ammai(k)=v_amma(k)
307         vitw_ammai(k)=vitw_amma(k,12)
308         ht_ammai(k)=ht_amma(k,12)
309         hq_ammai(k)=hq_amma(k,12)
310         vt_ammai(k)=0.
311         vq_ammai(k)=0.
312      enddo   
313      omega(:)=0.     
314      omega2(:)=0.
315      rho(:)=0.
316! vertical interpolation using TOGA interpolation routine:
317!      write(*,*)'avant interp vert', t_proftwp
318      CALL interp_toga_vertical(play,nlev_amma,plev_amma                    &
319     &         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai                 &
320     &         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai                         &
321     &         ,t_mod,q_mod,u_mod,v_mod,w_mod                               &
322     &         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
323!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
324
325! initial and boundary conditions :
326!      tsurf = ts_proftwp
327      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
328      do l = 1, llm
329! Ligne du dessous ?? decommenter si on lit theta au lieu de temp
330!      temp(l) = t_mod(l)*(play(l)/pzero)**rkappa
331       temp(l) = t_mod(l) 
332       q(l,1) = q_mod(l)
333       q(l,2) = 0.0
334!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
335       u(l) = u_mod(l)
336       v(l) = v_mod(l)
337       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
338       omega(l) = w_mod(l)*(-rg*rho(l))
339       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
340
341       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
342!on applique le forcage total au premier pas de temps
343!attention: signe different de toga
344       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
345!forcage en th
346!       d_t_adv(l) = ht_mod(l)
347       d_q_adv(l,1) = hq_mod(l)
348       d_q_adv(l,2) = 0.0
349       dt_cooling(l)=0.
350      enddo     
351       write(*,*) 'Prof initeforcing AMMA interpole temp39',temp(39)
352     
353
354!     ok_flux_surf=.false.
355      fsens=-1.*sens_amma(12)
356      flat=-1.*lat_amma(12)
357       
358      endif !forcing_amma
359
360
361!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
362!---------------------------------------------------------------------
363! Forcing from DICE experiment (see file DICE_protocol_vn2-3.pdf)
364!---------------------------------------------------------------------
365
366      if (forcing_dice) then
367!read DICE forcings
368      fich_dice='dice_driver.nc'
369      call read_dice(fich_dice,nlev_dice,nt_dice                    &
370     &     ,zz_dice,plev_dice,t_dice,qv_dice,u_dice,v_dice,o3_dice &
371     &     ,shf_dice,lhf_dice,lwup_dice,swup_dice,tg_dice,ustar_dice& 
372     &     ,psurf_dice,ug_dice,vg_dice,ht_dice,hq_dice              &
373     &     ,hu_dice,hv_dice,w_dice,omega_dice)
374
375      write(*,*) 'Forcing DICE lu'
376
377!champs initiaux:
378      do k=1,nlev_dice
379         t_dicei(k)=t_dice(k)
380         qv_dicei(k)=qv_dice(k)
381         u_dicei(k)=u_dice(k)
382         v_dicei(k)=v_dice(k)
383         o3_dicei(k)=o3_dice(k)
384         ht_dicei(k)=ht_dice(k,1)
385         hq_dicei(k)=hq_dice(k,1)
386         hu_dicei(k)=hu_dice(k,1)
387         hv_dicei(k)=hv_dice(k,1)
388         w_dicei(k)=w_dice(k,1)
389         omega_dicei(k)=omega_dice(k,1)
390      enddo   
391      omega(:)=0.     
392      omega2(:)=0.
393      rho(:)=0.
394! vertical interpolation using TOGA interpolation routine:
395!      write(*,*)'avant interp vert', t_proftwp
396!
397!     CALL interp_dice_time(daytime,day1,annee_ref
398!    i             ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice
399!    i             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice
400!    i             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice
401!    i             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice
402!    o             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
403!    o             ,ustar_prof,psurf_prof,ug_profd,vg_profd
404!    o             ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd
405!    o             ,omega_profd)
406
407      CALL interp_dice_vertical(play,nlev_dice,nt_dice,plev_dice       &
408     &         ,t_dicei,qv_dicei,u_dicei,v_dicei,o3_dicei             &
409     &         ,ht_dicei,hq_dicei,hu_dicei,hv_dicei,w_dicei,omega_dicei&
410     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                       &
411     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
412
413! Pour tester les advections horizontales de T et Q, on met w_mod et omega_mod ?? zero (MPL 20131108)
414!     w_mod(:,:)=0.
415!     omega_mod(:,:)=0.
416
417!       write(*,*) 'Profil initial forcing DICE interpole',t_mod
418! Les forcages DICE sont donnes /jour et non /seconde !
419      ht_mod(:)=ht_mod(:)/86400.
420      hq_mod(:)=hq_mod(:)/86400.
421      hu_mod(:)=hu_mod(:)/86400.
422      hv_mod(:)=hv_mod(:)/86400.
423
424! initial and boundary conditions :
425      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
426      do l = 1, llm
427! Ligne du dessous ?? decommenter si on lit theta au lieu de temp
428!      temp(l) = th_mod(l)*(play(l)/pzero)**rkappa
429       temp(l) = t_mod(l) 
430       q(l,1) = qv_mod(l)
431       q(l,2) = 0.0
432!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
433       u(l) = u_mod(l)
434       v(l) = v_mod(l)
435       ug(l)=ug_dice(1)
436       vg(l)=vg_dice(1)
437       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
438!      omega(l) = w_mod(l)*(-rg*rho(l))
439       omega(l) = omega_mod(l)
440       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
441
442       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
443!on applique le forcage total au premier pas de temps
444!attention: signe different de toga
445       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
446!forcage en th
447!       d_t_adv(l) = ht_mod(l)
448       d_q_adv(l,1) = hq_mod(l)
449       d_q_adv(l,2) = 0.0
450       dt_cooling(l)=0.
451      enddo     
452       write(*,*) 'Profil initial forcing DICE interpole temp39',temp(39)
453     
454
455!     ok_flux_surf=.false.
456      fsens=-1.*shf_dice(1)
457      flat=-1.*lhf_dice(1)
458! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par
459! le coefficient de trainee en surface cd**2=ustar*vent(k=1)
460! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface
461! MPL 05082013
462      ust=ustar_dice(1)
463      tg=tg_dice(1)
464      print *,'ust= ',ust
465      IF (tsurf .LE. 0.) THEN
466       tsurf= tg_dice(1)
467      ENDIF
468      psurf= psurf_dice(1)
469      solsw_in = (1.-albedo)/albedo*swup_dice(1)
470      sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1)
471      PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in
472      endif !forcing_dice
473
474!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
475!---------------------------------------------------------------------
476! Forcing from GABLS4 experiment
477!---------------------------------------------------------------------
478
479!!!! Si la temperature de surface n'est pas impos??e:
480 
481      if (forcing_gabls4) then
482!read GABLS4 forcings
483     
484      fich_gabls4='gabls4_driver.nc'
485     
486       
487      call read_gabls4(fich_gabls4,nlev_gabls4,nt_gabls4,nsol_gabls4,zz_gabls4,depth_sn_gabls4,ug_gabls4,vg_gabls4 &
488     & ,plev_gabls4,th_gabls4,t_gabls4,qv_gabls4,u_gabls4,v_gabls4,ht_gabls4,hq_gabls4,tg_gabls4,tsnow_gabls4,snow_dens_gabls4)
489
490      write(*,*) 'Forcing GABLS4 lu'
491
492!champs initiaux:
493      do k=1,nlev_gabls4
494         t_gabi(k)=t_gabls4(k)
495         qv_gabi(k)=qv_gabls4(k)
496         u_gabi(k)=u_gabls4(k)
497         v_gabi(k)=v_gabls4(k)
498         poub(k)=0.
499         ht_gabi(k)=ht_gabls4(k,1)
500         hq_gabi(k)=hq_gabls4(k,1)
501         ug_gabi(k)=ug_gabls4(k,1)
502         vg_gabi(k)=vg_gabls4(k,1)
503      enddo
504 
505      omega(:)=0.     
506      omega2(:)=0.
507      rho(:)=0.
508! vertical interpolation using TOGA interpolation routine:
509!      write(*,*)'avant interp vert', t_proftwp
510!
511!     CALL interp_dice_time(daytime,day1,annee_ref
512!    i             ,year_ini_dice,day_ju_ini_dice,nt_dice,dt_dice
513!    i             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice
514!    i             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice
515!    i             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice
516!    o             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof
517!    o             ,ustar_prof,psurf_prof,ug_profd,vg_profd
518!    o             ,ht_profd,hq_profd,hu_profd,hv_profd,w_profd
519!    o             ,omega_profd)
520
521      CALL interp_dice_vertical(play,nlev_gabls4,nt_gabls4,plev_gabls4       &
522     &         ,t_gabi,qv_gabi,u_gabi,v_gabi,poub                  &
523     &         ,ht_gabi,hq_gabi,ug_gabi,vg_gabi,poub,poub          &
524     &         ,t_mod,qv_mod,u_mod,v_mod,o3_mod                    &
525     &         ,ht_mod,hq_mod,ug_mod,vg_mod,w_mod,omega_mod,mxcalc)
526
527! Les forcages GABLS4 ont l air d etre en K/S quoiqu en dise le fichier gabls4_driver.nc !? MPL 20141024
528!     ht_mod(:)=ht_mod(:)/86400.
529!     hq_mod(:)=hq_mod(:)/86400.
530
531! initial and boundary conditions :
532      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
533      do l = 1, llm
534! Ligne du dessous ?? decommenter si on lit theta au lieu de temp
535!      temp(l) = th_mod(l)*(play(l)/pzero)**rkappa
536       temp(l) = t_mod(l) 
537       q(l,1) = qv_mod(l)
538       q(l,2) = 0.0
539!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
540       u(l) = u_mod(l)
541       v(l) = v_mod(l)
542       ug(l)=ug_mod(l)
543       vg(l)=vg_mod(l)
544       
545!
546!       tg=tsurf
547!       
548
549       print *,'***** tsurf=',tsurf
550       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
551!      omega(l) = w_mod(l)*(-rg*rho(l))
552       omega(l) = omega_mod(l)
553       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
554       
555   
556
557       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
558!on applique le forcage total au premier pas de temps
559!attention: signe different de toga
560!      d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
561!forcage en th
562       d_t_adv(l) = ht_mod(l)
563       d_q_adv(l,1) = hq_mod(l)
564       d_q_adv(l,2) = 0.0
565       dt_cooling(l)=0.
566      enddo     
567
568!--------------- Residus forcages du cas Dice (a supprimer) MPL 20141024---------------
569! Le cas Dice doit etre force avec ustar mais on peut simplifier en forcant par
570! le coefficient de trainee en surface cd**2=ustar*vent(k=1)
571! On commence ici a stocker ustar dans cdrag puis on terminera le calcul dans pbl_surface
572! MPL 05082013
573!     ust=ustar_dice(1)
574!     tg=tg_dice(1)
575!     print *,'ust= ',ust
576!     IF (tsurf .LE. 0.) THEN
577!      tsurf= tg_dice(1)
578!     ENDIF
579!     psurf= psurf_dice(1)
580!     solsw_in = (1.-albedo)/albedo*swup_dice(1)
581!     sollw_in = (0.7*RSIGMA*temp(1)**4)-lwup_dice(1)
582!     PRINT *,'1D_READ_FORC : solsw, sollw',solsw_in,sollw_in
583!--------------------------------------------------------------------------------------
584      endif !forcing_gabls4
585
586
587
588! Forcing from Arm_Cu case                   
589! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
590! large scale advective forcing,radiative forcing
591! and advective tendency of theta and qt to be applied
592!---------------------------------------------------------------------
593
594      if (forcing_armcu) then
595! read armcu forcing :
596       write(*,*) 'Avant lecture Forcing Arm_Cu'
597      fich_armcu = './ifa_armcu.txt'
598      CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu,                       &
599     & sens_armcu,flat_armcu,adv_theta_armcu,                               &
600     & rad_theta_armcu,adv_qt_armcu)
601       write(*,*) 'Forcing Arm_Cu lu'
602
603!----------------------------------------------------------------------
604! Read profiles from file: prof.inp.19 or prof.inp.40 
605! For this case, profiles are given for two vertical resolution
606! 19 or 40 levels
607!
608! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html
609! Note that the initial profiles contain no liquid water! 
610! (so potential temperature can be interpreted as liquid water
611! potential temperature and water vapor as total water)
612! profiles are given at full levels
613!----------------------------------------------------------------------
614
615      call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod,           &
616     &           v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)
617
618! time interpolation for initial conditions:
619      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
620
621      print *,'Avant interp_armcu_time'
622      print *,'daytime=',daytime
623      print *,'day1=',day1
624      print *,'annee_ref=',annee_ref
625      print *,'year_ini_armcu=',year_ini_armcu
626      print *,'day_ju_ini_armcu=',day_ju_ini_armcu
627      print *,'nt_armcu=',nt_armcu
628      print *,'dt_armcu=',dt_armcu
629      print *,'nlev_armcu=',nlev_armcu
630      CALL interp_armcu_time(daytime,day1,annee_ref                         &
631     &            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu        &
632     &            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu         &
633     &            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof         &
634     &            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
635       write(*,*) 'Forcages interpoles dans temps'
636
637! No vertical interpolation if nlev imposed to 19 or 40
638! The vertical grid stops at 4000m # 600hPa
639      mxcalc=llm
640
641! initial and boundary conditions :
642!     tsurf = ts_prof
643! tsurf read in lmdz1d.def
644      write(*,*) 'Tsurf initiale: ',tsurf
645      do l = 1, llm
646       play(l)=play_mod(l)*100.
647       presnivs(l)=play(l)
648       zlay(l)=height(l)
649       temp(l) = t_mod(l)
650       teta(l)=theta_mod(l)
651       q(l,1) = qv_mod(l)/1000.
652! No liquid water in the initial profil
653       q(l,2) = 0.
654       u(l) = u_mod(l)
655       ug(l)= u_mod(l)
656       v(l) = v_mod(l)
657       vg(l)= v_mod(l)
658! Advective forcings are given in K or g/kg ... per HOUR
659!      IF(height(l).LT.1000) THEN
660!        d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
661!        d_q_adv(l,1) = adv_qt_prof/1000./3600.
662!        d_q_adv(l,2) = 0.0
663!      ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
664!        d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*
665!    :               (1-(height(l)-1000.)/2000.)
666!        d_t_adv(l) = d_t_adv(l)/3600.
667!        d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
668!        d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
669!        d_q_adv(l,2) = 0.0
670!      ELSE
671!        d_t_adv(l) = 0.0
672!        d_q_adv(l,1) = 0.0
673!        d_q_adv(l,2) = 0.0
674!      ENDIF
675      enddo
676! plev at half levels is given in proh.inp.19 or proh.inp.40 files
677      plev(1)= ap(llm+1)+bp(llm+1)*psurf
678      do l = 1, llm
679      plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf
680      print *,'Read_forc: l height play plev zlay temp',                    &
681     &   l,height(l),play(l),plev(l),zlay(l),temp(l)
682      enddo
683! For this case, fluxes are imposed
684       fsens=-1*sens_prof
685       flat=-1*flat_prof
686
687      endif ! forcing_armcu
688!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
689!---------------------------------------------------------------------
690! Forcing from transition case of Irina Sandu                 
691!---------------------------------------------------------------------
692
693      if (forcing_sandu) then
694       write(*,*) 'Avant lecture Forcing SANDU'
695
696! read sanduref forcing :
697      fich_sandu = './ifa_sanduref.txt'
698      CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
699
700       write(*,*) 'Forcing SANDU lu'
701
702!----------------------------------------------------------------------
703! Read profiles from file: prof.inp.001 
704!----------------------------------------------------------------------
705
706      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,       &
707     &           thl_profs,q_profs,u_profs,v_profs,                         &
708     &           w_profs,omega_profs,o3mmr_profs)
709
710! time interpolation for initial conditions:
711      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
712! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
713! revoir 1DUTILS.h et les arguments
714
715      print *,'Avant interp_sandu_time'
716      print *,'daytime=',daytime
717      print *,'day1=',day1
718      print *,'annee_ref=',annee_ref
719      print *,'year_ini_sandu=',year_ini_sandu
720      print *,'day_ju_ini_sandu=',day_ju_ini_sandu
721      print *,'nt_sandu=',nt_sandu
722      print *,'dt_sandu=',dt_sandu
723      print *,'nlev_sandu=',nlev_sandu
724      CALL interp_sandu_time(daytime,day1,annee_ref                         &
725     &             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu       &
726     &             ,nlev_sandu                                              &
727     &             ,ts_sandu,ts_prof)
728
729! vertical interpolation:
730      print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
731      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs                 &
732     &         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs           &
733     &         ,omega_profs,o3mmr_profs                                     &
734     &         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod                       &
735     &         ,omega_mod,o3mmr_mod,mxcalc)
736       write(*,*) 'Profil initial forcing SANDU interpole'
737
738! initial and boundary conditions :
739      tsurf = ts_prof
740      write(*,*) 'SST initiale: ',tsurf
741      do l = 1, llm
742       temp(l) = t_mod(l)
743       tetal(l)=thl_mod(l)
744       q(l,1) = q_mod(l)
745       q(l,2) = 0.0
746       u(l) = u_mod(l)
747       v(l) = v_mod(l)
748       w(l) = w_mod(l)
749       omega(l) = omega_mod(l)
750       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
751!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
752!?       omega2(l)=-rho(l)*omega(l)
753       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
754!      d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
755!      d_q_adv(l,1) = vq_mod(l)
756       d_t_adv(l) = alpha*omega(l)/rcpd
757       d_q_adv(l,1) = 0.0
758       d_q_adv(l,2) = 0.0
759      enddo
760
761      endif ! forcing_sandu
762!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
763!---------------------------------------------------------------------
764! Forcing from Astex case
765!---------------------------------------------------------------------
766
767      if (forcing_astex) then
768       write(*,*) 'Avant lecture Forcing Astex'
769
770! read astex forcing :
771      fich_astex = './ifa_astex.txt'
772      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,    &
773     &  ug_astex,vg_astex,ufa_astex,vfa_astex)
774
775       write(*,*) 'Forcing Astex lu'
776
777!----------------------------------------------------------------------
778! Read profiles from file: prof.inp.001
779!----------------------------------------------------------------------
780
781      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,       &
782     &           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,      &
783     &           w_profa,tke_profa,o3mmr_profa)
784
785! time interpolation for initial conditions:
786      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
787! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
788! revoir 1DUTILS.h et les arguments
789
790      print *,'Avant interp_astex_time'
791      print *,'daytime=',daytime
792      print *,'day1=',day1
793      print *,'annee_ref=',annee_ref
794      print *,'year_ini_astex=',year_ini_astex
795      print *,'day_ju_ini_astex=',day_ju_ini_astex
796      print *,'nt_astex=',nt_astex
797      print *,'dt_astex=',dt_astex
798      print *,'nlev_astex=',nlev_astex
799      CALL interp_astex_time(daytime,day1,annee_ref                         &
800     &             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex       &
801     &             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex         &
802     &             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof    &
803     &             ,ufa_prof,vfa_prof)
804
805! vertical interpolation:
806      print *,'Avant interp_vertical: nlev_astex=',nlev_astex
807      CALL interp_astex_vertical(play,nlev_astex,plev_profa                 &
808     &         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa                &
809     &         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa               &
810     &         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod        &
811     &         ,tke_mod,o3mmr_mod,mxcalc)
812       write(*,*) 'Profil initial forcing Astex interpole'
813
814! initial and boundary conditions :
815      tsurf = ts_prof
816      write(*,*) 'SST initiale: ',tsurf
817      do l = 1, llm
818       temp(l) = t_mod(l)
819       tetal(l)=thl_mod(l)
820       q(l,1) = qv_mod(l)
821       q(l,2) = ql_mod(l)
822       u(l) = u_mod(l)
823       v(l) = v_mod(l)
824       w(l) = w_mod(l)
825       omega(l) = w_mod(l)
826!      omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
827!      rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
828!      omega2(l)=-rho(l)*omega(l)
829       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
830!      d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
831!      d_q_adv(l,1) = vq_mod(l)
832       d_t_adv(l) = alpha*omega(l)/rcpd
833       d_q_adv(l,1) = 0.0
834       d_q_adv(l,2) = 0.0
835      enddo
836
837      endif ! forcing_astex
838!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
839!---------------------------------------------------------------------
840! Forcing from standard case :
841!---------------------------------------------------------------------
842
843      if (forcing_case) then
844
845         write(*,*),'avant call read_1D_cas'
846         call read_1D_cas
847         write(*,*) 'Forcing read'
848
849!Time interpolation for initial conditions using TOGA interpolation routine
850         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
851      CALL interp_case_time(day,day1,annee_ref                                                              &
852!    &         ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                                         &
853     &         ,nt_cas,nlev_cas                                                                             &
854     &         ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas                                                     &
855     &         ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas                                                 &
856     &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                                         &
857     &         ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas,ustar_cas                                             &
858     &         ,uw_cas,vw_cas,q1_cas,q2_cas                                                                 &
859     &         ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas                       &
860     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas                   &
861     &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas      &
862     &         ,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas,ustar_prof_cas               &
863     &         ,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas)
864
865! vertical interpolation using TOGA interpolation routine:
866!      write(*,*)'avant interp vert', t_prof
867      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
868     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas    &
869     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
870     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
871     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas           &
872     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
873     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
874!       write(*,*) 'Profil initial forcing case interpole',t_mod
875
876! initial and boundary conditions :
877!      tsurf = ts_prof_cas
878      ts_cur = ts_prof_cas
879      psurf=plev_prof_cas(1)
880      write(*,*) 'SST initiale: ',tsurf
881      do l = 1, llm
882       temp(l) = t_mod_cas(l)
883       q(l,1) = q_mod_cas(l)
884       q(l,2) = 0.0
885       u(l) = u_mod_cas(l)
886       v(l) = v_mod_cas(l)
887       omega(l) = w_mod_cas(l)
888       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
889
890       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
891!on applique le forcage total au premier pas de temps
892!attention: signe different de toga
893       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
894       d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
895       d_q_adv(l,2) = 0.0
896       d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
897! correction bug d_u -> d_v (MM+MPL 20170310)
898       d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
899      enddo     
900
901! In case fluxes are imposed
902       IF (ok_flux_surf) THEN
903       fsens=sens_prof_cas
904       flat=lat_prof_cas
905       ENDIF
906       IF (ok_prescr_ust) THEN
907       ust=ustar_prof_cas
908       print *,'ust=',ust
909       ENDIF
910
911      endif !forcing_case
912!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
913!---------------------------------------------------------------------
914! Forcing from standard case :
915!---------------------------------------------------------------------
916
917      if (forcing_case2) then
918
919         write(*,*),'avant call read2_1D_cas'
920         call read2_1D_cas
921         write(*,*) 'Forcing read'
922
923!Time interpolation for initial conditions using interpolation routine
924         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
925        CALL interp2_case_time(daytime,day1,annee_ref                                       &
926!    &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas                           &
927     &       ,nt_cas,nlev_cas                                                               &
928     &       ,ts_cas,ps_cas,plev_cas,t_cas,th_cas,thv_cas,thl_cas,qv_cas,ql_cas,qi_cas      &
929     &       ,u_cas,v_cas,ug_cas,vg_cas,vitw_cas,omega_cas,du_cas,hu_cas,vu_cas             &
930     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas                           &
931     &       ,dq_cas,hq_cas,vq_cas,dth_cas,hth_cas,vth_cas,lat_cas,sens_cas,ustar_cas       &
932     &       ,uw_cas,vw_cas,q1_cas,q2_cas,tke_cas                                           &
933!
934     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,theta_prof_cas,thv_prof_cas  &
935     &       ,thl_prof_cas,qv_prof_cas,ql_prof_cas,qi_prof_cas                              &
936     &       ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas    &
937     &       ,du_prof_cas,hu_prof_cas,vu_prof_cas                                           &
938     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
939     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
940     &       ,dth_prof_cas,hth_prof_cas,vth_prof_cas,lat_prof_cas                           &
941     &       ,sens_prof_cas,ustar_prof_cas,uw_prof_cas,vw_prof_cas,q1_prof_cas,q2_prof_cas,tke_prof_cas)
942
943      do l = 1, nlev_cas
944      print *,'apres 1ere interp: plev_cas, plev_prof_cas=',l,plev_cas(l,1),plev_prof_cas(l)
945      enddo
946
947! vertical interpolation using interpolation routine:
948!      write(*,*)'avant interp vert', t_prof
949      CALL interp2_case_vertical(play,nlev_cas,plev_prof_cas                                              &
950     &         ,t_prof_cas,theta_prof_cas,thv_prof_cas,thl_prof_cas                                          &
951     &         ,qv_prof_cas,ql_prof_cas,qi_prof_cas,u_prof_cas,v_prof_cas                                 &
952     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,omega_prof_cas                                      &
953     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas                   &
954     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas    &
955     &         ,dth_prof_cas,hth_prof_cas,vth_prof_cas                                                    &
956!
957     &         ,t_mod_cas,theta_mod_cas,thv_mod_cas,thl_mod_cas,qv_mod_cas,ql_mod_cas,qi_mod_cas          &
958     &         ,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas,omega_mod_cas                         &
959     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas                         &
960     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas           &
961     &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
962
963!       write(*,*) 'Profil initial forcing case interpole',t_mod
964
965! initial and boundary conditions :
966!      tsurf = ts_prof_cas
967      ts_cur = ts_prof_cas
968      psurf=plev_prof_cas(1)
969      write(*,*) 'SST initiale: ',tsurf
970      do l = 1, llm
971       temp(l) = t_mod_cas(l)
972       q(l,1) = qv_mod_cas(l)
973       q(l,2) = ql_mod_cas(l)
974       u(l) = u_mod_cas(l)
975       ug(l)= ug_mod_cas(l)
976       v(l) = v_mod_cas(l)
977       vg(l)= vg_mod_cas(l)
978! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
979       omega(l) = omega_mod_cas(l)
980       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
981
982       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
983!on applique le forcage total au premier pas de temps
984!attention: signe different de toga
985       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
986       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
987!      d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
988       d_q_adv(l,1) = dq_mod_cas(l)
989       d_q_adv(l,2) = 0.0
990!      d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
991       d_u_adv(l) = du_mod_cas(l)
992!      d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
993! correction bug d_u -> d_v (MM+MPL 20170310)
994       d_v_adv(l) = dv_mod_cas(l)
995      enddo     
996
997! Faut-il multiplier par -1 ? (MPL 20160713)
998       IF (ok_flux_surf) THEN
999       fsens=-1.*sens_prof_cas
1000       flat=-1.*lat_prof_cas
1001       ENDIF
1002!
1003       IF (ok_prescr_ust) THEN
1004       ust=ustar_prof_cas
1005       print *,'ust=',ust
1006       ENDIF
1007
1008      endif !forcing_case2
1009!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1010
Note: See TracBrowser for help on using the repository browser.