Changeset 1795 for LMDZ5/branches/testing/libf/phy1d/1D_read_forc_cases.h

Timestamp:

Jul 18, 2013, 10:20:28 AM (11 years ago)

Author:

Ehouarn Millour

Message:

Version testing basee sur la r1794

---

Testing release based on r1794

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phy1d/1D_read_forc_cases.h (modified) (8 diffs)

LMDZ5/branches/testing/libf/phy1d/1D_read_forc_cases.h

@@ -4 +4 @@
 !----------------------------------------------------------------------
 
-      if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then
-
+      if (forcing_les .or. forcing_radconv
+     :    .or. forcing_GCSSold .or. forcing_fire) then
+
+      if (forcing_fire) then
+!----------------------------------------------------------------------
+!read fire forcings from fire.nc
+!----------------------------------------------------------------------
+      fich_fire='fire.nc'
+      call read_fire(fich_fire,nlev_fire,nt_fire
+     :     ,height,tttprof,qtprof,uprof,vprof,e12prof
+     :     ,ugprof,vgprof,wfls,dqtdxls
+     :     ,dqtdyls,dqtdtls,thlpcar)
+      write(*,*) 'Forcing FIRE lu'
+      kmax=120            ! nombre de niveaux dans les profils et forcages
+      else 
 !----------------------------------------------------------------------
 ! Read profiles from files: prof.inp.001 and lscale.inp.001
@@ -16 +29 @@
      .           wfls,dqtdxls,dqtdyls,dqtdtls,
      .           thlpcar) 
+      endif
 
 ! compute altitudes of play levels.
@@ -35 +49 @@
         frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
         ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
-        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
+       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
           temp(l) = ttt*(play(l)/pzero)**rkappa
           teta(l) = ttt
-        else
+       else
           temp(l) = ttt
           teta(l) = ttt*(pzero/play(l))**rkappa
-        endif
+       endif
           print *,' temp,teta ',l,temp(l),teta(l)
         q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
@@ -58 +72 @@
           if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
             ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
-        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
+       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
           temp(l) = ttt*(play(l)/pzero)**rkappa
           teta(l) = ttt
-        else
+       else
           temp(l) = ttt
           teta(l) = ttt*(pzero/play(l))**rkappa
-        endif
+       endif
           print *,' temp,teta ',l,temp(l),teta(l)
             q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
@@ -77 +91 @@
           elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
             ttt =tttprof(1)
-        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
+       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
           temp(l) = ttt*(play(l)/pzero)**rkappa
           teta(l) = ttt
-        else
+       else
           temp(l) = ttt
           teta(l) = ttt*(pzero/play(l))**rkappa
-        endif
+       endif
             q(l,1)  = qtprof(1)
             u(l)    =  uprof(1)
@@ -112 +126 @@
       enddo
 
-      endif ! forcing_les .or. forcing_GCSSold
+      endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
@@ -263 +277 @@
        
       endif !forcing_twpice
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Forcing from AMMA experiment (Couvreux et al. 2010) :
+!---------------------------------------------------------------------
+
+      if (forcing_amma) then
+!read AMMA forcings
+      fich_amma='amma.nc'
+      call read_amma(fich_amma,nlev_amma,nt_amma
+     :     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma
+     :     ,ht_amma,hq_amma,sens_amma,lat_amma)
+
+      write(*,*) 'Forcing AMMA lu'
+
+!champs initiaux:
+      do k=1,nlev_amma
+         th_ammai(k)=th_amma(k)
+         q_ammai(k)=q_amma(k)
+         u_ammai(k)=u_amma(k)
+         v_ammai(k)=v_amma(k)
+         vitw_ammai(k)=vitw_amma(k,12)
+         ht_ammai(k)=ht_amma(k,12)
+         hq_ammai(k)=hq_amma(k,12)
+         vt_ammai(k)=0.
+         vq_ammai(k)=0.
+      enddo   
+      omega(:)=0.      
+      omega2(:)=0.
+      rho(:)=0.
+! vertical interpolation using TOGA interpolation routine:
+!      write(*,*)'avant interp vert', t_proftwp
+      CALL interp_toga_vertical(play,nlev_amma,plev_amma
+     :         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai
+     :         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai
+     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
+     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
+
+! initial and boundary conditions :
+!      tsurf = ts_proftwp
+      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
+      do l = 1, llm
+! Ligne du dessous à decommenter si on lit theta au lieu de temp
+!      temp(l) = t_mod(l)*(play(l)/pzero)**rkappa 
+       temp(l) = t_mod(l) 
+       q(l,1) = q_mod(l)
+       q(l,2) = 0.0
+!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
+       omega(l) = w_mod(l)*(-rg*rho(l))
+       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+!on applique le forcage total au premier pas de temps
+!attention: signe different de toga
+       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
+!forcage en th
+!       d_th_adv(l) = ht_mod(l)
+       d_q_adv(l,1) = hq_mod(l)
+       d_q_adv(l,2) = 0.0
+       dt_cooling(l)=0.
+      enddo     
+       write(*,*) 'Profil initial forcing AMMA interpole temp39',
+     &           temp(39)
+      
+
+!     ok_flux_surf=.false.
+      fsens=-1.*sens_amma(12)
+      flat=-1.*lat_amma(12)
+       
+      endif !forcing_amma
+
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
@@ -366 +456 @@
       endif ! forcing_armcu
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!---------------------------------------------------------------------
+! Forcing from transition case of Irina Sandu                  
+!---------------------------------------------------------------------
+
+      if (forcing_sandu) then
+       write(*,*) 'Avant lecture Forcing SANDU'
+
+! read sanduref forcing :
+      fich_sandu = './ifa_sanduref.txt'
+      CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
+
+       write(*,*) 'Forcing SANDU lu'
+
+!----------------------------------------------------------------------
+! Read profiles from file: prof.inp.001 
+!----------------------------------------------------------------------
+
+      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,
+     .           thl_profs,q_profs,u_profs,v_profs,
+     .           w_profs,omega_profs,o3mmr_profs)
+
+! time interpolation for initial conditions:
+      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
+! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
+! revoir 1DUTILS.h et les arguments
+
+      print *,'Avant interp_sandu_time'
+      print *,'daytime=',daytime
+      print *,'day1=',day1
+      print *,'annee_ref=',annee_ref
+      print *,'year_ini_sandu=',year_ini_sandu
+      print *,'day_ju_ini_sandu=',day_ju_ini_sandu
+      print *,'nt_sandu=',nt_sandu
+      print *,'dt_sandu=',dt_sandu
+      print *,'nlev_sandu=',nlev_sandu
+      CALL interp_sandu_time(daytime,day1,annee_ref
+     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
+     i             ,nlev_sandu
+     i             ,ts_sandu,ts_prof)
+
+! vertical interpolation:
+      print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
+      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
+     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
+     :         ,omega_profs,o3mmr_profs
+     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
+     :         ,omega_mod,o3mmr_mod,mxcalc)
+       write(*,*) 'Profil initial forcing SANDU interpole'
+
+! initial and boundary conditions :
+      tsurf = ts_prof
+      write(*,*) 'SST initiale: ',tsurf
+      do l = 1, llm
+       temp(l) = t_mod(l)
+       tetal(l)=thl_mod(l)
+       q(l,1) = q_mod(l)
+       q(l,2) = 0.0
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       w(l) = w_mod(l)
+       omega(l) = omega_mod(l)
+       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
+!?       omega2(l)=-rho(l)*omega(l)
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
+!      d_q_adv(l,1) = vq_mod(l)
+       d_th_adv(l) = alpha*omega(l)/rcpd
+       d_q_adv(l,1) = 0.0
+       d_q_adv(l,2) = 0.0
+      enddo
+
+      endif ! forcing_sandu
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Forcing from Astex case
+!---------------------------------------------------------------------
+
+      if (forcing_astex) then
+       write(*,*) 'Avant lecture Forcing Astex'
+
+! read astex forcing :
+      fich_astex = './ifa_astex.txt'
+      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,
+     :  ug_astex,vg_astex,ufa_astex,vfa_astex)
+
+       write(*,*) 'Forcing Astex lu'
+
+!----------------------------------------------------------------------
+! Read profiles from file: prof.inp.001
+!----------------------------------------------------------------------
+
+      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,
+     .           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,
+     .           w_profa,tke_profa,o3mmr_profa)
+
+! time interpolation for initial conditions:
+      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
+! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
+! revoir 1DUTILS.h et les arguments
+
+      print *,'Avant interp_astex_time'
+      print *,'daytime=',daytime
+      print *,'day1=',day1
+      print *,'annee_ref=',annee_ref
+      print *,'year_ini_astex=',year_ini_astex
+      print *,'day_ju_ini_astex=',day_ju_ini_astex
+      print *,'nt_astex=',nt_astex
+      print *,'dt_astex=',dt_astex
+      print *,'nlev_astex=',nlev_astex
+      CALL interp_astex_time(daytime,day1,annee_ref
+     i             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex
+     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
+     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
+     i             ,ufa_prof,vfa_prof)
+
+! vertical interpolation:
+      print *,'Avant interp_vertical: nlev_astex=',nlev_astex
+      CALL interp_astex_vertical(play,nlev_astex,plev_profa
+     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
+     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
+     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
+     :         ,tke_mod,o3mmr_mod,mxcalc)
+       write(*,*) 'Profil initial forcing Astex interpole'
+
+! initial and boundary conditions :
+      tsurf = ts_prof
+      write(*,*) 'SST initiale: ',tsurf
+      do l = 1, llm
+       temp(l) = t_mod(l)
+       tetal(l)=thl_mod(l)
+       q(l,1) = qv_mod(l)
+       q(l,2) = ql_mod(l)
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       w(l) = w_mod(l)
+       omega(l) = w_mod(l)
+!      omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+!      rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
+!      omega2(l)=-rho(l)*omega(l)
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
+!      d_q_adv(l,1) = vq_mod(l)
+       d_th_adv(l) = alpha*omega(l)/rcpd
+       d_q_adv(l,1) = 0.0
+       d_q_adv(l,2) = 0.0
+      enddo
+
+      endif ! forcing_astex
+