Ignore:
Timestamp:
Jun 29, 2017, 11:58:07 AM (7 years ago)
Author:
fhourdin
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h

    r2716 r2920  
    3737
    3838       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    39        d_th_adv(l) = ht_gcssold(l)
     39       d_t_adv(l) = ht_gcssold(l)
    4040       d_q_adv(l,1) = hq_gcssold(l)
    4141       dt_cooling(l) = 0.0
    … …  
    8484
    8585       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    86        d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
     86       d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
    8787       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
    8888       dt_cooling(l) = 0.0
    … …  
    183183
    184184       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    185        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
     185       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
    186186       d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
    187187       d_u_adv(l) = hu_mod(l)-d_u_dyn_z(l)
    … …  
    224224       ug(l)= ug_mod(l)
    225225       vg(l)= vg_mod(l)
    226        d_th_adv(l)=ht_mod(l)
     226       d_t_adv(l)=ht_mod(l)
    227227       d_q_adv(l,1)=hq_mod(l)
    228228      enddo
    … …  
    315315!calcul de l'advection totale
    316316        if (cptadvw) then
    317         d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
     317        d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-d_t_dyn_z(l)
    318318!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
    319319        d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l)
    320320!        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
    321321        else
    322         d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
     322        d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
    323323        d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
    324324        endif
    … …  
    392392       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    393393!calcul de l'advection totale
    394 !        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
     394!        d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
    395395!attention: on impose dth
    396         d_th_adv(l) = alpha*omega(l)/rcpd+                                  &
     396        d_t_adv(l) = alpha*omega(l)/rcpd+                                  &
    397397     &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
    398 !        d_th_adv(l) = 0.
     398!        d_t_adv(l) = 0.
    399399!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
    400400        d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l)
    … …  
    417417!---------------------------------------------------------------------
    418418      if (forcing_rico) then
    419 !       call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,
    420 !     :  q,temp,u,v,play)
     419!      call lstendH(llm,omega,dt_dyn,dq_dyn,du_dyn, dv_dyn,q,temp,u,v,play)
    421420       call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play)
    422421
    423422        do l=1,llm
    424        d_th_adv(l) =  (dth_rico(l) +  dt_dyn(l))
     423       d_t_adv(l) =  (dth_rico(l) +  dt_dyn(l))
    425424       d_q_adv(l,1) = (dqh_rico(l) +  dq_dyn(l,1))
    426425       d_q_adv(l,2) = 0.
    … …  
    458457       vg(l)= v_mod(l)
    459458       IF((phi(l)/RG).LT.1000) THEN
    460          d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
     459         d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
    461460         d_q_adv(l,1) = adv_qt_prof/1000./3600.
    462461         d_q_adv(l,2) = 0.0
    463462!        print *,'INF1000: phi dth dq1 dq2',
    464 !    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
     463!    :  phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
    465464       ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN
    466465         fact=((phi(l)/RG)-1000.)/2000.
    467466         fact=1-fact
    468          d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
     467         d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.
    469468         d_q_adv(l,1) = adv_qt_prof*fact/1000./3600.
    470469         d_q_adv(l,2) = 0.0
    471470!        print *,'SUP1000: phi fact dth dq1 dq2',
    472 !    :  phi(l)/RG,fact,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
     471!    :  phi(l)/RG,fact,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
    473472       ELSE
    474          d_th_adv(l) = 0.0
     473         d_t_adv(l) = 0.0
    475474         d_q_adv(l,1) = 0.0
    476475         d_q_adv(l,2) = 0.0
    477476!        print *,'SUP3000: phi dth dq1 dq2',
    478 !    :  phi(l)/RG,d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
     477!    :  phi(l)/RG,d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
    479478       ENDIF
    480479      dt_cooling(l) = 0.0 
    481480!     print *,'Interp armcu: phi dth dq1 dq2',
    482 !    :  l,phi(l),d_th_adv(l),d_q_adv(l,1),d_q_adv(l,2)
     481!    :  l,phi(l),d_t_adv(l),d_q_adv(l,1),d_q_adv(l,2)
    483482      enddo
    484483      endif ! forcing_armcu
    … …  
    554553       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    555554!
    556 !      d_th_adv(l) = 0.0
     555!      d_t_adv(l) = 0.0
    557556!      d_q_adv(l,1) = 0.0
    558557!CR:test advection=0
    559558!calcul de l'advection verticale
    560         d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
     559        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
    561560!        print*,'temp adv',l,-d_t_dyn_z(l)
    562561        d_q_adv(l,1) = -d_q_dyn_z(l)
    … …  
    637636       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    638637!
    639 !      d_th_adv(l) = 0.0
     638!      d_t_adv(l) = 0.0
    640639!      d_q_adv(l,1) = 0.0
    641640!CR:test advection=0
    642641!calcul de l'advection verticale
    643         d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
     642        d_t_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
    644643!        print*,'temp adv',l,-d_t_dyn_z(l)
    645644        d_q_adv(l,1) = -d_q_dyn_z(l)
    … …  
    715714
    716715!Calcul de l advection verticale
     716
    717717      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
     718
    718719      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
    719720      d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
    … …  
    764765
    765766      do l = 1, llm
    766        omega(l) = w_mod_cas(l)
     767       omega(l) = w_mod_cas(l)  ! juste car w_mod_cas en Pa/s (MPL 20170310)
    767768       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
    768769       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    … …  
    782783
    783784        if ((tend_t.eq.1).and.(tend_w.eq.0)) then
    784 !           d_th_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
    785            d_th_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
     785!           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
     786           d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
    786787        else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
    787 !           d_th_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
    788            d_th_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
     788!           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
     789           d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
    789790        endif
    790791
    … …  
    816817      ENDIF
    817818      endif ! forcing_case
    818 
    819819
    820820!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    877877      d_th_z(:)=0.
    878878      d_q_z(:)=0.
     879      d_u_z(:)=0.
     880      d_v_z(:)=0.
    879881      d_t_dyn_z(:)=0.
    880882      d_th_dyn_z(:)=0.
    881883      d_q_dyn_z(:)=0.
     884      d_u_dyn_z(:)=0.
     885      d_v_dyn_z(:)=0.
    882886      DO l=2,llm-1
    883887       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
    884888       d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1))
    885889       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
     890       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
     891       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
    886892      ENDDO
    887893      d_t_z(1)=d_t_z(2)
    888894      d_th_z(1)=d_th_z(2)
    889895      d_q_z(1)=d_q_z(2)
     896      d_u_z(1)=d_u_z(2)
     897      d_v_z(1)=d_v_z(2)
    890898      d_t_z(llm)=d_t_z(llm-1)
    891899      d_th_z(llm)=d_th_z(llm-1)
    892900      d_q_z(llm)=d_q_z(llm-1)
     901      d_u_z(llm)=d_u_z(llm-1)
     902      d_v_z(llm)=d_v_z(llm-1)
    893903
    894904!Calcul de l advection verticale
    895       d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
    896       d_th_dyn_z(:)=w_mod_cas(:)*d_th_z(:)
    897       d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
    898 
     905! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170310)
     906      d_t_dyn_z(:)=omega_mod_cas(:)*d_t_z(:)
     907      d_th_dyn_z(:)=omega_mod_cas(:)*d_th_z(:)
     908      d_q_dyn_z(:)=omega_mod_cas(:)*d_q_z(:)
     909      d_u_dyn_z(:)=omega_mod_cas(:)*d_u_z(:)
     910      d_v_dyn_z(:)=omega_mod_cas(:)*d_v_z(:)
     911
     912!geostrophic wind
     913      if (forc_geo.eq.1) then
     914        do l=1,llm
     915        ug(l) = ug_mod_cas(l)
     916        vg(l) = vg_mod_cas(l)
     917        enddo
     918      endif
    899919!wind nudging
    900920      if (nudging_u.gt.0.) then
    … …  
    902922           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
    903923        enddo
    904       else
    905         do l=1,llm
    906         ug(l) = u_mod_cas(l)
    907         enddo
     924!     else
     925!       do l=1,llm
     926!          u(l) = u_mod_cas(l)
     927!       enddo
    908928      endif
    909929
    … …  
    912932           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
    913933        enddo
    914       else
    915         do l=1,llm
    916         vg(l) = v_mod_cas(l)
    917         enddo
     934!     else
     935!       do l=1,llm
     936!          v(l) = v_mod_cas(l)
     937!       enddo
    918938      endif
    919939
    … …  
    922942           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
    923943        enddo
    924       else
    925         do l=1,llm
    926         w(l) = w_mod_cas(l)
    927         enddo
     944 !    else
     945 !      do l=1,llm
     946 !         w(l) = w_mod_cas(l)
     947 !      enddo
    928948      endif
    929949
    … …  
    941961
    942962      do l = 1, llm
    943        omega(l) = w_mod_cas(l)
     963! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
     964       omega(l) = omega_mod_cas(l)
    944965       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
    945966       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    946967
    947 !calcul advection
    948 !       if ((tend_u.eq.1).and.(tend_w.eq.0)) then
    949 !          d_u_adv(l)=du_mod_cas(l)
    950 !       else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
    951 !          d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
     968!calcul advections
     969        if ((forc_u.eq.1).and.(forc_w.eq.0)) then
     970           d_u_adv(l)=du_mod_cas(l)
     971        else if ((forc_u.eq.1).and.(forc_w.eq.1)) then
     972           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
     973        endif
     974
     975        if ((forc_v.eq.1).and.(forc_w.eq.0)) then
     976           d_v_adv(l)=dv_mod_cas(l)
     977        else if ((forc_v.eq.1).and.(forc_w.eq.1)) then
     978           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
     979        endif
     980
     981! Puisque dth a ete converti en dt, on traite de la meme facon
     982! les flags tadv et thadv
     983        if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.0)) then
     984!          d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
     985           d_t_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
     986        else if ((tadv.eq.1.or.thadv.eq.1) .and. (forc_w.eq.1)) then
     987!          d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
     988           d_t_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
     989        endif
     990
     991!       if ((thadv.eq.1) .and. (forc_w.eq.0)) then
     992!          d_t_adv(l)=alpha*omega(l)/rcpd-dth_mod_cas(l)
     993!          d_t_adv(l)=alpha*omega(l)/rcpd+dth_mod_cas(l)
     994!       else if ((thadv.eq.1) .and. (forc_w.eq.1)) then
     995!          d_t_adv(l)=alpha*omega(l)/rcpd-hth_mod_cas(l)-d_t_dyn_z(l)
     996!          d_t_adv(l)=alpha*omega(l)/rcpd+hth_mod_cas(l)-d_t_dyn_z(l)
    952997!       endif
    953 !
    954 !       if ((tend_v.eq.1).and.(tend_w.eq.0)) then
    955 !          d_v_adv(l)=dv_mod_cas(l)
    956 !       else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
    957 !          d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
    958 !       endif
    959 !
    960 !-----------------------------------------------------
    961         if (tadv.eq.1 .or. tadvh.eq.1) then
    962            d_t_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
    963         else if (tadvv.eq.1) then
    964 ! ATTENTION d_t_dyn_z pas calcule (voir twpice)
    965            d_t_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
    966         endif
    967         print *,'interp_case d_t_dyn_z=',d_t_dyn_z(l),d_q_dyn_z(l)
    968 
    969 ! Verifier le signe !!
    970         if (thadv.eq.1 .or. thadvh.eq.1) then
    971            d_th_adv(l)=dth_mod_cas(l)
    972            print *,'dthadv=',d_th_adv(l)*86400.
    973         else if (thadvv.eq.1) then
    974            d_th_adv(l)=hth_mod_cas(l)-d_th_dyn_z(l)
    975         endif
    976  
    977 ! Verifier le signe !!
    978         if ((qadv.eq.1).and.(forc_w.eq.0)) then
     998
     999        if ((qadv.eq.1) .and. (forc_w.eq.0)) then
    9791000           d_q_adv(l,1)=dq_mod_cas(l)
    980         else if ((qadvh.eq.1).and.(forc_w.eq.1)) then
     1001!          d_q_adv(l,1)=-1*dq_mod_cas(l)
     1002        else if ((qadv.eq.1) .and. (forc_w.eq.1)) then
    9811003           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
     1004!          d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
    9821005        endif
    9831006         
Note: See TracChangeset for help on using the changeset viewer.