Ignore:
Timestamp:
Jul 5, 2013, 10:38:13 AM (11 years ago)
Author:
Laurent Fairhead
Message:

Inclusion des cas AMMA et Sandu dans la configuration 1D
MPL

AMMA and Sandu cases included in 1D configuration
MPL

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phy1d/1D_read_forc_cases.h

    r1607 r1780  
    44!----------------------------------------------------------------------
    55
    6       if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then
    7 
     6      if (forcing_les .or. forcing_radconv
     7     :    .or. forcing_GCSSold .or. forcing_fire) then
     8
     9      if (forcing_fire) then
     10!----------------------------------------------------------------------
     11!read fire forcings from fire.nc
     12!----------------------------------------------------------------------
     13      fich_fire='fire.nc'
     14      call read_fire(fich_fire,nlev_fire,nt_fire
     15     :     ,height,tttprof,qtprof,uprof,vprof,e12prof
     16     :     ,ugprof,vgprof,wfls,dqtdxls
     17     :     ,dqtdyls,dqtdtls,thlpcar)
     18      write(*,*) 'Forcing FIRE lu'
     19      kmax=120            ! nombre de niveaux dans les profils et forcages
     20      else
    821!----------------------------------------------------------------------
    922! Read profiles from files: prof.inp.001 and lscale.inp.001
    … …  
    1629     .           wfls,dqtdxls,dqtdyls,dqtdtls,
    1730     .           thlpcar)
     31      endif
    1832
    1933! compute altitudes of play levels.
    … …  
    3549        frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
    3650        ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
    37         if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
     51       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
    3852          temp(l) = ttt*(play(l)/pzero)**rkappa
    3953          teta(l) = ttt
    40         else
     54       else
    4155          temp(l) = ttt
    4256          teta(l) = ttt*(pzero/play(l))**rkappa
    43         endif
     57       endif
    4458          print *,' temp,teta ',l,temp(l),teta(l)
    4559        q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
    … …  
    5872          if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
    5973            ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
    60         if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
     74       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
    6175          temp(l) = ttt*(play(l)/pzero)**rkappa
    6276          teta(l) = ttt
    63         else
     77       else
    6478          temp(l) = ttt
    6579          teta(l) = ttt*(pzero/play(l))**rkappa
    66         endif
     80       endif
    6781          print *,' temp,teta ',l,temp(l),teta(l)
    6882            q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
    … …  
    7791          elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
    7892            ttt =tttprof(1)
    79         if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
     93       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
    8094          temp(l) = ttt*(play(l)/pzero)**rkappa
    8195          teta(l) = ttt
    82         else
     96       else
    8397          temp(l) = ttt
    8498          teta(l) = ttt*(pzero/play(l))**rkappa
    85         endif
     99       endif
    86100            q(l,1)  = qtprof(1)
    87101            u(l)    =  uprof(1)
    … …  
    112126      enddo
    113127
    114       endif ! forcing_les .or. forcing_GCSSold
     128      endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire
    115129!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    116130!---------------------------------------------------------------------
    … …  
    263277       
    264278      endif !forcing_twpice
     279
     280!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     281!---------------------------------------------------------------------
     282! Forcing from AMMA experiment (Couvreux et al. 2010) :
     283!---------------------------------------------------------------------
     284
     285      if (forcing_amma) then
     286!read AMMA forcings
     287      fich_amma='amma.nc'
     288      call read_amma(fich_amma,nlev_amma,nt_amma
     289     :     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma
     290     :     ,ht_amma,hq_amma,sens_amma,lat_amma)
     291
     292      write(*,*) 'Forcing AMMA lu'
     293
     294!champs initiaux:
     295      do k=1,nlev_amma
     296         th_ammai(k)=th_amma(k)
     297         q_ammai(k)=q_amma(k)
     298         u_ammai(k)=u_amma(k)
     299         v_ammai(k)=v_amma(k)
     300         vitw_ammai(k)=vitw_amma(k,12)
     301         ht_ammai(k)=ht_amma(k,12)
     302         hq_ammai(k)=hq_amma(k,12)
     303         vt_ammai(k)=0.
     304         vq_ammai(k)=0.
     305      enddo   
     306      omega(:)=0.     
     307      omega2(:)=0.
     308      rho(:)=0.
     309! vertical interpolation using TOGA interpolation routine:
     310!      write(*,*)'avant interp vert', t_proftwp
     311      CALL interp_toga_vertical(play,nlev_amma,plev_amma
     312     :         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai
     313     :         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai
     314     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
     315     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
     316!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod
     317
     318! initial and boundary conditions :
     319!      tsurf = ts_proftwp
     320      write(*,*) 'SST initiale mxcalc: ',tsurf,mxcalc
     321      do l = 1, llm
     322! Ligne du dessous à decommenter si on lit theta au lieu de temp
     323!      temp(l) = t_mod(l)*(play(l)/pzero)**rkappa
     324       temp(l) = t_mod(l)
     325       q(l,1) = q_mod(l)
     326       q(l,2) = 0.0
     327!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
     328       u(l) = u_mod(l)
     329       v(l) = v_mod(l)
     330       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
     331       omega(l) = w_mod(l)*(-rg*rho(l))
     332       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     333
     334       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     335!on applique le forcage total au premier pas de temps
     336!attention: signe different de toga
     337       d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
     338!forcage en th
     339!       d_th_adv(l) = ht_mod(l)
     340       d_q_adv(l,1) = hq_mod(l)
     341       d_q_adv(l,2) = 0.0
     342       dt_cooling(l)=0.
     343      enddo     
     344       write(*,*) 'Profil initial forcing AMMA interpole temp39',
     345     &           temp(39)
     346     
     347
     348!     ok_flux_surf=.false.
     349      fsens=-1.*sens_amma(12)
     350      flat=-1.*lat_amma(12)
     351       
     352      endif !forcing_amma
     353
     354
    265355!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    266356!---------------------------------------------------------------------
    … …  
    366456      endif ! forcing_armcu
    367457!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    368 
     458!---------------------------------------------------------------------
     459! Forcing from transition case of Irina Sandu                 
     460!---------------------------------------------------------------------
     461
     462      if (forcing_sandu) then
     463       write(*,*) 'Avant lecture Forcing SANDU'
     464
     465! read sanduref forcing :
     466      fich_sandu = './ifa_sanduref.txt'
     467      CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
     468
     469       write(*,*) 'Forcing SANDU lu'
     470
     471!----------------------------------------------------------------------
     472! Read profiles from file: prof.inp.001
     473!----------------------------------------------------------------------
     474
     475      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,
     476     .           thl_profs,q_profs,u_profs,v_profs,
     477     .           w_profs,omega_profs,o3mmr_profs)
     478
     479! time interpolation for initial conditions:
     480      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
     481! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
     482! revoir 1DUTILS.h et les arguments
     483
     484      print *,'Avant interp_sandu_time'
     485      print *,'daytime=',daytime
     486      print *,'day1=',day1
     487      print *,'annee_ref=',annee_ref
     488      print *,'year_ini_sandu=',year_ini_sandu
     489      print *,'day_ju_ini_sandu=',day_ju_ini_sandu
     490      print *,'nt_sandu=',nt_sandu
     491      print *,'dt_sandu=',dt_sandu
     492      print *,'nlev_sandu=',nlev_sandu
     493      CALL interp_sandu_time(daytime,day1,annee_ref
     494     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
     495     i             ,nlev_sandu
     496     i             ,ts_sandu,ts_prof)
     497
     498! vertical interpolation:
     499      print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
     500      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
     501     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
     502     :         ,omega_profs,o3mmr_profs
     503     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
     504     :         ,omega_mod,o3mmr_mod,mxcalc)
     505       write(*,*) 'Profil initial forcing SANDU interpole'
     506
     507! initial and boundary conditions :
     508      tsurf = ts_prof
     509      write(*,*) 'SST initiale: ',tsurf
     510      do l = 1, llm
     511       temp(l) = t_mod(l)
     512       tetal(l)=thl_mod(l)
     513       q(l,1) = q_mod(l)
     514       q(l,2) = 0.0
     515       u(l) = u_mod(l)
     516       v(l) = v_mod(l)
     517       w(l) = w_mod(l)
     518       omega(l) = omega_mod(l)
     519       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     520!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
     521!?       omega2(l)=-rho(l)*omega(l)
     522       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     523!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
     524!      d_q_adv(l,1) = vq_mod(l)
     525       d_th_adv(l) = alpha*omega(l)/rcpd
     526       d_q_adv(l,1) = 0.0
     527       d_q_adv(l,2) = 0.0
     528      enddo
     529
     530      endif ! forcing_sandu
     531!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     532!---------------------------------------------------------------------
     533! Forcing from Astex case
     534!---------------------------------------------------------------------
     535
     536      if (forcing_astex) then
     537       write(*,*) 'Avant lecture Forcing Astex'
     538
     539! read astex forcing :
     540      fich_astex = './ifa_astex.txt'
     541      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,
     542     :  ug_astex,vg_astex,ufa_astex,vfa_astex)
     543
     544       write(*,*) 'Forcing Astex lu'
     545
     546!----------------------------------------------------------------------
     547! Read profiles from file: prof.inp.001
     548!----------------------------------------------------------------------
     549
     550      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,
     551     .           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,
     552     .           w_profa,tke_profa,o3mmr_profa)
     553
     554! time interpolation for initial conditions:
     555      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
     556! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
     557! revoir 1DUTILS.h et les arguments
     558
     559      print *,'Avant interp_astex_time'
     560      print *,'daytime=',daytime
     561      print *,'day1=',day1
     562      print *,'annee_ref=',annee_ref
     563      print *,'year_ini_astex=',year_ini_astex
     564      print *,'day_ju_ini_astex=',day_ju_ini_astex
     565      print *,'nt_astex=',nt_astex
     566      print *,'dt_astex=',dt_astex
     567      print *,'nlev_astex=',nlev_astex
     568      CALL interp_astex_time(daytime,day1,annee_ref
     569     i             ,year_ini_astex,day_ju_ini_astex,nt_astex,dt_astex
     570     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
     571     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
     572     i             ,ufa_prof,vfa_prof)
     573
     574! vertical interpolation:
     575      print *,'Avant interp_vertical: nlev_astex=',nlev_astex
     576      CALL interp_astex_vertical(play,nlev_astex,plev_profa
     577     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
     578     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
     579     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
     580     :         ,tke_mod,o3mmr_mod,mxcalc)
     581       write(*,*) 'Profil initial forcing Astex interpole'
     582
     583! initial and boundary conditions :
     584      tsurf = ts_prof
     585      write(*,*) 'SST initiale: ',tsurf
     586      do l = 1, llm
     587       temp(l) = t_mod(l)
     588       tetal(l)=thl_mod(l)
     589       q(l,1) = qv_mod(l)
     590       q(l,2) = ql_mod(l)
     591       u(l) = u_mod(l)
     592       v(l) = v_mod(l)
     593       w(l) = w_mod(l)
     594       omega(l) = w_mod(l)
     595!      omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     596!      rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
     597!      omega2(l)=-rho(l)*omega(l)
     598       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     599!      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
     600!      d_q_adv(l,1) = vq_mod(l)
     601       d_th_adv(l) = alpha*omega(l)/rcpd
     602       d_q_adv(l,1) = 0.0
     603       d_q_adv(l,2) = 0.0
     604      enddo
     605
     606      endif ! forcing_astex
     607
Note: See TracChangeset for help on using the changeset viewer.