source: LMDZ5/trunk/libf/phy1d/1D_read_forc_cases.h @ 1953

Last change on this file since 1953 was 1953, checked in by fhourdin, 11 years ago

Correction in for readprofiles_fire

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 23.7 KB
Line 
1!----------------------------------------------------------------------
2! forcing_les = .T. : Impose a constant cooling
3! forcing_radconv = .T. : Pure radiative-convective equilibrium:
4!----------------------------------------------------------------------
5
6
7      nq1=0
8      nq2=0
9
10      if (forcing_les .or. forcing_radconv
11     :    .or. forcing_GCSSold .or. forcing_fire) then
12
13      if (forcing_fire) then
14!----------------------------------------------------------------------
15!read fire forcings from fire.nc
16!----------------------------------------------------------------------
17      fich_fire='fire.nc'
18      call read_fire(fich_fire,nlev_fire,nt_fire
19     :     ,height,tttprof,qtprof,uprof,vprof,e12prof
20     :     ,ugprof,vgprof,wfls,dqtdxls
21     :     ,dqtdyls,dqtdtls,thlpcar)
22      write(*,*) 'Forcing FIRE lu'
23      kmax=120            ! nombre de niveaux dans les profils et forcages
24      else 
25!----------------------------------------------------------------------
26! Read profiles from files: prof.inp.001 and lscale.inp.001
27! (repris de readlesfiles)
28!----------------------------------------------------------------------
29
30      call readprofiles(nlev_max,kmax,nqtot,height,
31     .           tttprof,qtprof,uprof,vprof,
32     .           e12prof,ugprof,vgprof,
33     .           wfls,dqtdxls,dqtdyls,dqtdtls,
34     .           thlpcar,qprof,nq1,nq2)
35      endif
36
37! compute altitudes of play levels.
38      zlay(1) =zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
39      do l = 2,llm
40        zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf)
41      enddo
42
43!----------------------------------------------------------------------
44! Interpolation of the profiles given on the input file to
45! model levels
46!----------------------------------------------------------------------
47      zlay(1) = zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
48      do l=1,llm
49        ! Above the max altutide of the input file
50
51        if (zlay(l)<height(kmax)) mxcalc=l
52
53        frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
54        ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
55       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
56          temp(l) = ttt*(play(l)/pzero)**rkappa
57          teta(l) = ttt
58       else
59          temp(l) = ttt
60          teta(l) = ttt*(pzero/play(l))**rkappa
61       endif
62          print *,' temp,teta ',l,temp(l),teta(l)
63        q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
64        u(l)    =  uprof(kmax)-frac*(  uprof(kmax)-  uprof(kmax-1))
65        v(l)    =  vprof(kmax)-frac*(  vprof(kmax)-  vprof(kmax-1))
66        ug(l)   = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1))
67        vg(l)   = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1))
68        IF (nq2>0) q(l,nq1:nq2)=qprof(kmax,nq1:nq2)
69     s               -frac*(qprof(kmax,nq1:nq2)-qprof(kmax-1,nq1:nq2))
70        omega(l)=   wfls(kmax)-frac*(   wfls(kmax)-   wfls(kmax-1))
71
72        dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1))
73        dt_cooling(l)
74     .          =thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
75        do k=2,kmax
76          frac = (height(k)-zlay(l))/(height(k)-height(k-1))
77          if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
78          if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
79            ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
80       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
81          temp(l) = ttt*(play(l)/pzero)**rkappa
82          teta(l) = ttt
83       else
84          temp(l) = ttt
85          teta(l) = ttt*(pzero/play(l))**rkappa
86       endif
87          print *,' temp,teta ',l,temp(l),teta(l)
88            q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
89            u(l)    =  uprof(k)-frac*(  uprof(k)-  uprof(k-1))
90            v(l)    =  vprof(k)-frac*(  vprof(k)-  vprof(k-1))
91            ug(l)   = ugprof(k)-frac*( ugprof(k)- ugprof(k-1))
92            vg(l)   = vgprof(k)-frac*( vgprof(k)- vgprof(k-1))
93            IF (nq2>0) q(l,nq1:nq2)=qprof(k,nq1:nq2)
94     s                   -frac*(qprof(k,nq1:nq2)-qprof(k-1,nq1:nq2))
95            omega(l)=   wfls(k)-frac*(   wfls(k)-   wfls(k-1))
96            dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1))
97            dt_cooling(l)
98     .              =thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
99          elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
100            ttt =tttprof(1)
101       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
102          temp(l) = ttt*(play(l)/pzero)**rkappa
103          teta(l) = ttt
104       else
105          temp(l) = ttt
106          teta(l) = ttt*(pzero/play(l))**rkappa
107       endif
108            q(l,1)  = qtprof(1)
109            u(l)    =  uprof(1)
110            v(l)    =  vprof(1)
111            ug(l)   = ugprof(1)
112            vg(l)   = vgprof(1)
113            omega(l)=   wfls(1)
114            IF (nq2>0) q(l,nq1:nq2)=qprof(1,nq1:nq2)
115            dq_dyn(l,1)  =dqtdtls(1)
116            dt_cooling(l)=thlpcar(1)
117          endif
118        enddo
119
120        temp(l)=max(min(temp(l),350.),150.)
121        rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
122        if (l .lt. llm) then
123          zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l))
124        endif
125        omega2(l)=-rho(l)*omega(l)
126        omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s
127        if (l>1) then
128        if(zlay(l-1)>height(kmax)) then
129           omega(l)=0.0
130           omega2(l)=0.0
131        endif   
132        endif
133        if(q(l,1)<0.) q(l,1)=0.0
134        q(l,2)  = 0.0
135      enddo
136
137      endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire
138!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139!---------------------------------------------------------------------
140! Forcing for GCSSold:
141!---------------------------------------------------------------------
142      if (forcing_GCSSold) then
143       fich_gcssold_ctl = './forcing.ctl'
144       fich_gcssold_dat = './forcing8.dat'
145       call copie(llm,play,psurf,fich_gcssold_ctl)
146       call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,
147     :               ht_gcssold,hq_gcssold,hw_gcssold,
148     :               hu_gcssold,hv_gcssold,
149     :               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,
150     :               imp_fcg_gcssold,ts_fcg_gcssold,
151     :               Tp_fcg_gcssold,Turb_fcg_gcssold)
152       print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold
153      endif ! forcing_GCSSold
154!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
155!---------------------------------------------------------------------
156! Forcing for RICO:
157!---------------------------------------------------------------------
158      if (forcing_rico) then
159
160!       call writefield_phy('omega', omega,llm+1)
161      fich_rico = 'rico.txt'
162       call read_rico(fich_rico,nlev_rico,ps_rico,play
163     :             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
164     :             ,dth_rico,dqh_rico)
165        print*, ' on a lu et prepare RICO'
166
167       mxcalc=llm
168       print *, airefi, ' airefi '
169       do l = 1, llm
170       rho(l)  = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l)))
171       temp(l) = t_rico(l)
172       q(l,1) = q_rico(l)
173       q(l,2) = 0.0
174       u(l) = u_rico(l)
175       v(l) = v_rico(l)
176       ug(l)=u_rico(l)
177       vg(l)=v_rico(l)
178       omega(l) = -w_rico(l)*rg
179       omega2(l) = omega(l)/rg*airefi
180       enddo
181      endif
182!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183!---------------------------------------------------------------------
184! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) :
185!---------------------------------------------------------------------
186
187      if (forcing_toga) then
188
189! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps):
190      fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt'
191      CALL read_togacoare(fich_toga,nlev_toga,nt_toga
192     :         ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
193     :         ,ht_toga,vt_toga,hq_toga,vq_toga)
194
195       write(*,*) 'Forcing TOGA lu'
196
197! time interpolation for initial conditions:
198      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
199      CALL interp_toga_time(daytime,day1,annee_ref
200     i             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga
201     i             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga
202     i             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga
203     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
204     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
205
206! vertical interpolation:
207      CALL interp_toga_vertical(play,nlev_toga,plev_prof
208     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
209     :         ,ht_prof,vt_prof,hq_prof,vq_prof
210     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
211     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
212       write(*,*) 'Profil initial forcing TOGA interpole'
213
214! initial and boundary conditions :
215      tsurf = ts_prof
216      write(*,*) 'SST initiale: ',tsurf
217      do l = 1, llm
218       temp(l) = t_mod(l)
219       q(l,1) = q_mod(l)
220       q(l,2) = 0.0
221       u(l) = u_mod(l)
222       v(l) = v_mod(l)
223       omega(l) = w_mod(l)
224       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
225!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
226!?       omega2(l)=-rho(l)*omega(l)
227       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
228       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
229       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
230       d_q_adv(l,2) = 0.0
231      enddo
232
233      endif ! forcing_toga
234!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
235!---------------------------------------------------------------------
236! Forcing from TWPICE experiment (Shaocheng et al. 2010) :
237!---------------------------------------------------------------------
238
239      if (forcing_twpice) then
240!read TWP-ICE forcings
241      fich_twpice=
242     : 'd_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
243      call read_twpice(fich_twpice,nlev_twpi,nt_twpi
244     :     ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
245     :     ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)
246
247      write(*,*) 'Forcing TWP-ICE lu'
248!Time interpolation for initial conditions using TOGA interpolation routine
249         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
250      CALL interp_toga_time(daytime,day1,annee_ref
251     i          ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi
252     i             ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
253     i             ,ht_twpi,vt_twpi,hq_twpi,vq_twpi
254     o             ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp
255     o             ,u_proftwp,v_proftwp,w_proftwp
256     o             ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
257
258! vertical interpolation using TOGA interpolation routine:
259!      write(*,*)'avant interp vert', t_proftwp
260      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp
261     :         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp
262     :         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp
263     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
264     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
265!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
266
267! initial and boundary conditions :
268!      tsurf = ts_proftwp
269      write(*,*) 'SST initiale: ',tsurf
270      do l = 1, llm
271       temp(l) = t_mod(l)
272       q(l,1) = q_mod(l)
273       q(l,2) = 0.0
274       u(l) = u_mod(l)
275       v(l) = v_mod(l)
276       omega(l) = w_mod(l)
277       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
278
279       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
280!on applique le forcage total au premier pas de temps
281!attention: signe different de toga
282       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
283       d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
284       d_q_adv(l,2) = 0.0
285      enddo     
286       
287      endif !forcing_twpice
288
289!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
290!---------------------------------------------------------------------
291! Forcing from AMMA experiment (Couvreux et al. 2010) :
292!---------------------------------------------------------------------
293
294      if (forcing_amma) then
295!read AMMA forcings
296      fich_amma='amma.nc'
297      call read_amma(fich_amma,nlev_amma,nt_amma
298     :     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma
299     :     ,ht_amma,hq_amma,sens_amma,lat_amma)
300
301      write(*,*) 'Forcing AMMA lu'
302
303!champs initiaux:
304      do k=1,nlev_amma
305         th_ammai(k)=th_amma(k)
306         q_ammai(k)=q_amma(k)
307         u_ammai(k)=u_amma(k)
308         v_ammai(k)=v_amma(k)
309         vitw_ammai(k)=vitw_amma(k,12)
310         ht_ammai(k)=ht_amma(k,12)
311         hq_ammai(k)=hq_amma(k,12)
312         vt_ammai(k)=0.
313         vq_ammai(k)=0.
314      enddo   
315      omega(:)=0.     
316      omega2(:)=0.
317      rho(:)=0.
318! vertical interpolation using TOGA interpolation routine:
319!      write(*,*)'avant interp vert', t_proftwp
320      CALL interp_toga_vertical(play,nlev_amma,plev_amma
321     :         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai
322     :         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai
323     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
324     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
325!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
326
327! initial and boundary conditions :
328!      tsurf = ts_proftwp
329      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
330      do l = 1, llm
331! Ligne du dessous à decommenter si on lit theta au lieu de temp
332!      temp(l) = t_mod(l)*(play(l)/pzero)**rkappa
333       temp(l) = t_mod(l) 
334       q(l,1) = q_mod(l)
335       q(l,2) = 0.0
336!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
337       u(l) = u_mod(l)
338       v(l) = v_mod(l)
339       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
340       omega(l) = w_mod(l)*(-rg*rho(l))
341       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
342
343       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
344!on applique le forcage total au premier pas de temps
345!attention: signe different de toga
346       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
347!forcage en th
348!       d_th_adv(l) = ht_mod(l)
349       d_q_adv(l,1) = hq_mod(l)
350       d_q_adv(l,2) = 0.0
351       dt_cooling(l)=0.
352      enddo     
353       write(*,*) 'Profil initial forcing AMMA interpole temp39',
354     &           temp(39)
355     
356
357!     ok_flux_surf=.false.
358      fsens=-1.*sens_amma(12)
359      flat=-1.*lat_amma(12)
360       
361      endif !forcing_amma
362
363
364!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
365!---------------------------------------------------------------------
366! Forcing from Arm_Cu case                   
367! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
368! large scale advective forcing,radiative forcing
369! and advective tendency of theta and qt to be applied
370!---------------------------------------------------------------------
371
372      if (forcing_armcu) then
373! read armcu forcing :
374       write(*,*) 'Avant lecture Forcing Arm_Cu'
375      fich_armcu = './ifa_armcu.txt'
376      CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu,
377     : sens_armcu,flat_armcu,adv_theta_armcu,
378     : rad_theta_armcu,adv_qt_armcu)
379       write(*,*) 'Forcing Arm_Cu lu'
380
381!----------------------------------------------------------------------
382! Read profiles from file: prof.inp.19 or prof.inp.40 
383! For this case, profiles are given for two vertical resolution
384! 19 or 40 levels
385!
386! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html
387! Note that the initial profiles contain no liquid water! 
388! (so potential temperature can be interpreted as liquid water
389! potential temperature and water vapor as total water)
390! profiles are given at full levels
391!----------------------------------------------------------------------
392
393      call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod,
394     .           v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)
395
396! time interpolation for initial conditions:
397      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
398
399      print *,'Avant interp_armcu_time'
400      print *,'daytime=',daytime
401      print *,'day1=',day1
402      print *,'annee_ref=',annee_ref
403      print *,'year_ini_armcu=',year_ini_armcu
404      print *,'day_ju_ini_armcu=',day_ju_ini_armcu
405      print *,'nt_armcu=',nt_armcu
406      print *,'dt_armcu=',dt_armcu
407      print *,'nlev_armcu=',nlev_armcu
408      CALL interp_armcu_time(daytime,day1,annee_ref
409     i            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu
410     i            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu
411     i            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof
412     i            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
413       write(*,*) 'Forcages interpoles dans temps'
414
415! No vertical interpolation if nlev imposed to 19 or 40
416! The vertical grid stops at 4000m # 600hPa
417      mxcalc=llm
418
419! initial and boundary conditions :
420!     tsurf = ts_prof
421! tsurf read in lmdz1d.def
422      write(*,*) 'Tsurf initiale: ',tsurf
423      do l = 1, llm
424       play(l)=play_mod(l)*100.
425       presnivs(l)=play(l)
426       zlay(l)=height(l)
427       temp(l) = t_mod(l)
428       teta(l)=theta_mod(l)
429       q(l,1) = qv_mod(l)/1000.
430! No liquid water in the initial profil
431       q(l,2) = 0.
432       u(l) = u_mod(l)
433       ug(l)= u_mod(l)
434       v(l) = v_mod(l)
435       vg(l)= v_mod(l)
436! Advective forcings are given in K or g/kg ... per HOUR
437!      IF(height(l).LT.1000) THEN
438!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
439!        d_q_adv(l,1) = adv_qt_prof/1000./3600.
440!        d_q_adv(l,2) = 0.0
441!      ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
442!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*
443!    :               (1-(height(l)-1000.)/2000.)
444!        d_th_adv(l) = d_th_adv(l)/3600.
445!        d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
446!        d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
447!        d_q_adv(l,2) = 0.0
448!      ELSE
449!        d_th_adv(l) = 0.0
450!        d_q_adv(l,1) = 0.0
451!        d_q_adv(l,2) = 0.0
452!      ENDIF
453      enddo
454! plev at half levels is given in proh.inp.19 or proh.inp.40 files
455      plev(1)= ap(llm+1)+bp(llm+1)*psurf
456      do l = 1, llm
457      plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf
458      print *,'Read_forc: l height play plev zlay temp',
459     :   l,height(l),play(l),plev(l),zlay(l),temp(l)
460      enddo
461! For this case, fluxes are imposed
462       fsens=-1*sens_prof
463       flat=-1*flat_prof
464
465      endif ! forcing_armcu
466!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
467!---------------------------------------------------------------------
468! Forcing from transition case of Irina Sandu                 
469!---------------------------------------------------------------------
470
471      if (forcing_sandu) then
472       write(*,*) 'Avant lecture Forcing SANDU'
473
474! read sanduref forcing :
475      fich_sandu = './ifa_sanduref.txt'
476      CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
477
478       write(*,*) 'Forcing SANDU lu'
479
480!----------------------------------------------------------------------
481! Read profiles from file: prof.inp.001 
482!----------------------------------------------------------------------
483
484      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,
485     .           thl_profs,q_profs,u_profs,v_profs,
486     .           w_profs,omega_profs,o3mmr_profs)
487
488! time interpolation for initial conditions:
489      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
490! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
491! revoir 1DUTILS.h et les arguments
492
493      print *,'Avant interp_sandu_time'
494      print *,'daytime=',daytime
495      print *,'day1=',day1
496      print *,'annee_ref=',annee_ref
497      print *,'year_ini_sandu=',year_ini_sandu
498      print *,'day_ju_ini_sandu=',day_ju_ini_sandu
499      print *,'nt_sandu=',nt_sandu
500      print *,'dt_sandu=',dt_sandu
501      print *,'nlev_sandu=',nlev_sandu
502      CALL interp_sandu_time(daytime,day1,annee_ref
503     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
504     i             ,nlev_sandu
505     i             ,ts_sandu,ts_prof)
506
507! vertical interpolation:
508      print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
509      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
510     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
511     :         ,omega_profs,o3mmr_profs
512     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
513     :         ,omega_mod,o3mmr_mod,mxcalc)
514       write(*,*) 'Profil initial forcing SANDU interpole'
515
516! initial and boundary conditions :
517      tsurf = ts_prof
518      write(*,*) 'SST initiale: ',tsurf
519      do l = 1, llm
520       temp(l) = t_mod(l)
521       tetal(l)=thl_mod(l)
522       q(l,1) = q_mod(l)
523       q(l,2) = 0.0
524       u(l) = u_mod(l)
525       v(l) = v_mod(l)
526       w(l) = w_mod(l)
527       omega(l) = omega_mod(l)
528       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
529!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
530!?       omega2(l)=-rho(l)*omega(l)
531       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
532!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
533!      d_q_adv(l,1) = vq_mod(l)
534       d_th_adv(l) = alpha*omega(l)/rcpd
535       d_q_adv(l,1) = 0.0
536       d_q_adv(l,2) = 0.0
537      enddo
538
539      endif ! forcing_sandu
540!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
541!---------------------------------------------------------------------
542! Forcing from Astex case
543!---------------------------------------------------------------------
544
545      if (forcing_astex) then
546       write(*,*) 'Avant lecture Forcing Astex'
547
548! read astex forcing :
549      fich_astex = './ifa_astex.txt'
550      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,
551     :  ug_astex,vg_astex,ufa_astex,vfa_astex)
552
553       write(*,*) 'Forcing Astex lu'
554
555!----------------------------------------------------------------------
556! Read profiles from file: prof.inp.001
557!----------------------------------------------------------------------
558
559      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,
560     .           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,
561     .           w_profa,tke_profa,o3mmr_profa)
562
563! time interpolation for initial conditions:
564      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
565! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
566! revoir 1DUTILS.h et les arguments
567
568      print *,'Avant interp_astex_time'
569      print *,'daytime=',daytime
570      print *,'day1=',day1
571      print *,'annee_ref=',annee_ref
572      print *,'year_ini_astex=',year_ini_astex
573      print *,'day_ju_ini_astex=',day_ju_ini_astex
574      print *,'nt_astex=',nt_astex
575      print *,'dt_astex=',dt_astex
576      print *,'nlev_astex=',nlev_astex
577      CALL interp_astex_time(daytime,day1,annee_ref
578     i             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex
579     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
580     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
581     i             ,ufa_prof,vfa_prof)
582
583! vertical interpolation:
584      print *,'Avant interp_vertical: nlev_astex=',nlev_astex
585      CALL interp_astex_vertical(play,nlev_astex,plev_profa
586     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
587     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
588     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
589     :         ,tke_mod,o3mmr_mod,mxcalc)
590       write(*,*) 'Profil initial forcing Astex interpole'
591
592! initial and boundary conditions :
593      tsurf = ts_prof
594      write(*,*) 'SST initiale: ',tsurf
595      do l = 1, llm
596       temp(l) = t_mod(l)
597       tetal(l)=thl_mod(l)
598       q(l,1) = qv_mod(l)
599       q(l,2) = ql_mod(l)
600       u(l) = u_mod(l)
601       v(l) = v_mod(l)
602       w(l) = w_mod(l)
603       omega(l) = w_mod(l)
604!      omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
605!      rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
606!      omega2(l)=-rho(l)*omega(l)
607       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
608!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
609!      d_q_adv(l,1) = vq_mod(l)
610       d_th_adv(l) = alpha*omega(l)/rcpd
611       d_q_adv(l,1) = 0.0
612       d_q_adv(l,2) = 0.0
613      enddo
614
615      endif ! forcing_astex
616
Note: See TracBrowser for help on using the repository browser.