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

Last change on this file since 1615 was 1607, checked in by lguez, 13 years ago

Import 1D files. Files added to directory "phy1d" were in directories:

lmdz1d_source_20111207/phy1d_source
lmdz1d_source_20111207/phy1d_source_upd

extracted from:

http://www.lmd.jussieu.fr/~jyg/lmdz1d_source_20111207.tar.gz

File size: 14.6 KB
Line 
1!----------------------------------------------------------------------
2! forcing_les = .T. : Impose a constant cooling
3! forcing_radconv = .T. : Pure radiative-convective equilibrium:
4!----------------------------------------------------------------------
5
6      if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then
7
8!----------------------------------------------------------------------
9! Read profiles from files: prof.inp.001 and lscale.inp.001
10! (repris de readlesfiles)
11!----------------------------------------------------------------------
12
13      call readprofiles(nlev_max,kmax,height,
14     .           tttprof,qtprof,uprof,vprof,
15     .           e12prof,ugprof,vgprof,
16     .           wfls,dqtdxls,dqtdyls,dqtdtls,
17     .           thlpcar) 
18
19! compute altitudes of play levels.
20      zlay(1) =zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
21      do l = 2,llm
22        zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf)
23      enddo
24
25!----------------------------------------------------------------------
26! Interpolation of the profiles given on the input file to
27! model levels
28!----------------------------------------------------------------------
29      zlay(1) = zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
30      do l=1,llm
31        ! Above the max altutide of the input file
32
33        if (zlay(l)<height(kmax)) mxcalc=l
34
35        frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
36        ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
37        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
38          temp(l) = ttt*(play(l)/pzero)**rkappa
39          teta(l) = ttt
40        else
41          temp(l) = ttt
42          teta(l) = ttt*(pzero/play(l))**rkappa
43        endif
44          print *,' temp,teta ',l,temp(l),teta(l)
45        q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
46        u(l)    =  uprof(kmax)-frac*(  uprof(kmax)-  uprof(kmax-1))
47        v(l)    =  vprof(kmax)-frac*(  vprof(kmax)-  vprof(kmax-1))
48        ug(l)   = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1))
49        vg(l)   = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1))
50        omega(l)=   wfls(kmax)-frac*(   wfls(kmax)-   wfls(kmax-1))
51
52        dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1))
53        dt_cooling(l)
54     .          =thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
55        do k=2,kmax
56          frac = (height(k)-zlay(l))/(height(k)-height(k-1))
57          if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
58          if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
59            ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
60        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
61          temp(l) = ttt*(play(l)/pzero)**rkappa
62          teta(l) = ttt
63        else
64          temp(l) = ttt
65          teta(l) = ttt*(pzero/play(l))**rkappa
66        endif
67          print *,' temp,teta ',l,temp(l),teta(l)
68            q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
69            u(l)    =  uprof(k)-frac*(  uprof(k)-  uprof(k-1))
70            v(l)    =  vprof(k)-frac*(  vprof(k)-  vprof(k-1))
71            ug(l)   = ugprof(k)-frac*( ugprof(k)- ugprof(k-1))
72            vg(l)   = vgprof(k)-frac*( vgprof(k)- vgprof(k-1))
73            omega(l)=   wfls(k)-frac*(   wfls(k)-   wfls(k-1))
74            dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1))
75            dt_cooling(l)
76     .              =thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
77          elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
78            ttt =tttprof(1)
79        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
80          temp(l) = ttt*(play(l)/pzero)**rkappa
81          teta(l) = ttt
82        else
83          temp(l) = ttt
84          teta(l) = ttt*(pzero/play(l))**rkappa
85        endif
86            q(l,1)  = qtprof(1)
87            u(l)    =  uprof(1)
88            v(l)    =  vprof(1)
89            ug(l)   = ugprof(1)
90            vg(l)   = vgprof(1)
91            omega(l)=   wfls(1)
92            dq_dyn(l,1)  =dqtdtls(1)
93            dt_cooling(l)=thlpcar(1)
94          endif
95        enddo
96
97        temp(l)=max(min(temp(l),350.),150.)
98        rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
99        if (l .lt. llm) then
100          zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l))
101        endif
102        omega2(l)=-rho(l)*omega(l)
103        omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s
104        if (l>1) then
105        if(zlay(l-1)>height(kmax)) then
106           omega(l)=0.0
107           omega2(l)=0.0
108        endif   
109        endif
110        if(q(l,1)<0.) q(l,1)=0.0
111        q(l,2)  = 0.0
112      enddo
113
114      endif ! forcing_les .or. forcing_GCSSold
115!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116!---------------------------------------------------------------------
117! Forcing for GCSSold:
118!---------------------------------------------------------------------
119      if (forcing_GCSSold) then
120       fich_gcssold_ctl = './forcing.ctl'
121       fich_gcssold_dat = './forcing8.dat'
122       call copie(llm,play,psurf,fich_gcssold_ctl)
123       call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,
124     :               ht_gcssold,hq_gcssold,hw_gcssold,
125     :               hu_gcssold,hv_gcssold,
126     :               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,
127     :               imp_fcg_gcssold,ts_fcg_gcssold,
128     :               Tp_fcg_gcssold,Turb_fcg_gcssold)
129       print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold
130      endif ! forcing_GCSSold
131!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
132!---------------------------------------------------------------------
133! Forcing for RICO:
134!---------------------------------------------------------------------
135      if (forcing_rico) then
136
137!       call writefield_phy('omega', omega,llm+1)
138      fich_rico = 'rico.txt'
139       call read_rico(fich_rico,nlev_rico,ps_rico,play
140     :             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
141     :             ,dth_rico,dqh_rico)
142        print*, ' on a lu et prepare RICO'
143
144       mxcalc=llm
145       print *, airefi, ' airefi '
146       do l = 1, llm
147       rho(l)  = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l)))
148       temp(l) = t_rico(l)
149       q(l,1) = q_rico(l)
150       q(l,2) = 0.0
151       u(l) = u_rico(l)
152       v(l) = v_rico(l)
153       ug(l)=u_rico(l)
154       vg(l)=v_rico(l)
155       omega(l) = -w_rico(l)*rg
156       omega2(l) = omega(l)/rg*airefi
157       enddo
158      endif
159!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
160!---------------------------------------------------------------------
161! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) :
162!---------------------------------------------------------------------
163
164      if (forcing_toga) then
165
166! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps):
167      fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt'
168      CALL read_togacoare(fich_toga,nlev_toga,nt_toga
169     :         ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
170     :         ,ht_toga,vt_toga,hq_toga,vq_toga)
171
172       write(*,*) 'Forcing TOGA lu'
173
174! time interpolation for initial conditions:
175      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
176      CALL interp_toga_time(daytime,day1,annee_ref
177     i             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga
178     i             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga
179     i             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga
180     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
181     o             ,ht_prof,vt_prof,hq_prof,vq_prof)
182
183! vertical interpolation:
184      CALL interp_toga_vertical(play,nlev_toga,plev_prof
185     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
186     :         ,ht_prof,vt_prof,hq_prof,vq_prof
187     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
188     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
189       write(*,*) 'Profil initial forcing TOGA interpole'
190
191! initial and boundary conditions :
192      tsurf = ts_prof
193      write(*,*) 'SST initiale: ',tsurf
194      do l = 1, llm
195       temp(l) = t_mod(l)
196       q(l,1) = q_mod(l)
197       q(l,2) = 0.0
198       u(l) = u_mod(l)
199       v(l) = v_mod(l)
200       omega(l) = w_mod(l)
201       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
202!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
203!?       omega2(l)=-rho(l)*omega(l)
204       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
205       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
206       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
207       d_q_adv(l,2) = 0.0
208      enddo
209
210      endif ! forcing_toga
211!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
212!---------------------------------------------------------------------
213! Forcing from TWPICE experiment (Shaocheng et al. 2010) :
214!---------------------------------------------------------------------
215
216      if (forcing_twpice) then
217!read TWP-ICE forcings
218      fich_twpice=
219     : 'd_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
220      call read_twpice(fich_twpice,nlev_twpi,nt_twpi
221     :     ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
222     :     ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)
223
224      write(*,*) 'Forcing TWP-ICE lu'
225!Time interpolation for initial conditions using TOGA interpolation routine
226         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
227      CALL interp_toga_time(daytime,day1,annee_ref
228     i          ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi
229     i             ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
230     i             ,ht_twpi,vt_twpi,hq_twpi,vq_twpi
231     o             ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp
232     o             ,u_proftwp,v_proftwp,w_proftwp
233     o             ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)
234
235! vertical interpolation using TOGA interpolation routine:
236!      write(*,*)'avant interp vert', t_proftwp
237      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp
238     :         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp
239     :         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp
240     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
241     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
242!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
243
244! initial and boundary conditions :
245!      tsurf = ts_proftwp
246      write(*,*) 'SST initiale: ',tsurf
247      do l = 1, llm
248       temp(l) = t_mod(l)
249       q(l,1) = q_mod(l)
250       q(l,2) = 0.0
251       u(l) = u_mod(l)
252       v(l) = v_mod(l)
253       omega(l) = w_mod(l)
254       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
255
256       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
257!on applique le forcage total au premier pas de temps
258!attention: signe different de toga
259       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
260       d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
261       d_q_adv(l,2) = 0.0
262      enddo     
263       
264      endif !forcing_twpice
265!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
266!---------------------------------------------------------------------
267! Forcing from Arm_Cu case                   
268! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
269! large scale advective forcing,radiative forcing
270! and advective tendency of theta and qt to be applied
271!---------------------------------------------------------------------
272
273      if (forcing_armcu) then
274! read armcu forcing :
275       write(*,*) 'Avant lecture Forcing Arm_Cu'
276      fich_armcu = './ifa_armcu.txt'
277      CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu,
278     : sens_armcu,flat_armcu,adv_theta_armcu,
279     : rad_theta_armcu,adv_qt_armcu)
280       write(*,*) 'Forcing Arm_Cu lu'
281
282!----------------------------------------------------------------------
283! Read profiles from file: prof.inp.19 or prof.inp.40 
284! For this case, profiles are given for two vertical resolution
285! 19 or 40 levels
286!
287! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html
288! Note that the initial profiles contain no liquid water! 
289! (so potential temperature can be interpreted as liquid water
290! potential temperature and water vapor as total water)
291! profiles are given at full levels
292!----------------------------------------------------------------------
293
294      call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod,
295     .           v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)
296
297! time interpolation for initial conditions:
298      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
299
300      print *,'Avant interp_armcu_time'
301      print *,'daytime=',daytime
302      print *,'day1=',day1
303      print *,'annee_ref=',annee_ref
304      print *,'year_ini_armcu=',year_ini_armcu
305      print *,'day_ju_ini_armcu=',day_ju_ini_armcu
306      print *,'nt_armcu=',nt_armcu
307      print *,'dt_armcu=',dt_armcu
308      print *,'nlev_armcu=',nlev_armcu
309      CALL interp_armcu_time(daytime,day1,annee_ref
310     i            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu
311     i            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu
312     i            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof
313     i            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
314       write(*,*) 'Forcages interpoles dans temps'
315
316! No vertical interpolation if nlev imposed to 19 or 40
317! The vertical grid stops at 4000m # 600hPa
318      mxcalc=llm
319
320! initial and boundary conditions :
321!     tsurf = ts_prof
322! tsurf read in lmdz1d.def
323      write(*,*) 'Tsurf initiale: ',tsurf
324      do l = 1, llm
325       play(l)=play_mod(l)*100.
326       presnivs(l)=play(l)
327       zlay(l)=height(l)
328       temp(l) = t_mod(l)
329       teta(l)=theta_mod(l)
330       q(l,1) = qv_mod(l)/1000.
331! No liquid water in the initial profil
332       q(l,2) = 0.
333       u(l) = u_mod(l)
334       ug(l)= u_mod(l)
335       v(l) = v_mod(l)
336       vg(l)= v_mod(l)
337! Advective forcings are given in K or g/kg ... per HOUR
338!      IF(height(l).LT.1000) THEN
339!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
340!        d_q_adv(l,1) = adv_qt_prof/1000./3600.
341!        d_q_adv(l,2) = 0.0
342!      ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
343!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*
344!    :               (1-(height(l)-1000.)/2000.)
345!        d_th_adv(l) = d_th_adv(l)/3600.
346!        d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
347!        d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
348!        d_q_adv(l,2) = 0.0
349!      ELSE
350!        d_th_adv(l) = 0.0
351!        d_q_adv(l,1) = 0.0
352!        d_q_adv(l,2) = 0.0
353!      ENDIF
354      enddo
355! plev at half levels is given in proh.inp.19 or proh.inp.40 files
356      plev(1)= ap(llm+1)+bp(llm+1)*psurf
357      do l = 1, llm
358      plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf
359      print *,'Read_forc: l height play plev zlay temp',
360     :   l,height(l),play(l),plev(l),zlay(l),temp(l)
361      enddo
362! For this case, fluxes are imposed
363       fsens=-1*sens_prof
364       flat=-1*flat_prof
365
366      endif ! forcing_armcu
367!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
368
Note: See TracBrowser for help on using the repository browser.