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

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

Location:
LMDZ5/trunk/libf/phylmd/dyn1d
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r2904 r2920  
    482482       forc_omega =0
    483483       CALL getin('forc_omega',forc_omega)
     484
     485!Config  Key  = forc_u
     486!Config  Desc = forcage ou non par u
     487!Config  Def  = false
     488!Config  Help = forcage ou non par u
     489       forc_u =0
     490       CALL getin('forc_u',forc_u)
     491
     492!Config  Key  = forc_v
     493!Config  Desc = forcage ou non par v
     494!Config  Def  = false
     495!Config  Help = forcage ou non par v
     496       forc_v =0
     497       CALL getin('forc_v',forc_v)
    484498
    485499!Config  Key  = forc_w
    … …  
    653667
    654668      modname = 'dyn1deta0 : '
    655 !!      nmq(1)="vap"
    656 !!      nmq(2)="cond"
    657 !!      do iq=3,nqtot
    658 !!        write(nmq(iq),'("tra",i1)') iq-2
    659 !!      enddo
    660       DO iq = 1,nqtot
    661         nmq(iq) = trim(tname(iq))
    662       ENDDO
     669      nmq(1)="vap"
     670      nmq(2)="cond"
     671      do iq=3,nqtot
     672        write(nmq(iq),'("tra",i1)') iq-2
     673      enddo
    663674      print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot
    664675      CALL open_startphy(fichnom)
    … …  
    807818      CALL open_restartphy(fichnom)
    808819      print*,'redm1 ',fichnom,klon,klev,nqtot
    809 !!      nmq(1)="vap"
    810 !!      nmq(2)="cond"
    811 !!      nmq(3)="tra1"
    812 !!      nmq(4)="tra2"
    813       DO iq = 1,nqtot
    814         nmq(iq) = trim(tname(iq))
    815       ENDDO
     820      nmq(1)="vap"
     821      nmq(2)="cond"
     822      nmq(3)="tra1"
     823      nmq(4)="tra2"
    816824
    817825      modname = 'dyn1dredem'
    … …  
    13961404      cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1)))
    13971405
    1398 
    13991406      do k=2,llm-1
    14001407
    … …  
    27892796         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
    27902797         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
    2791          dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
    27922798     
    27932799         else !play>plev_prof_cas(1)
    … …  
    28162822         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
    28172823         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
    2818          dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
    28192824
    28202825         endif ! play.le.plev_prof_cas(1)
    … …  
    28452850         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                            !jyg
    28462851         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                            !jyg
    2847          dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact                      !jyg
    28482852 
    28492853        endif ! play
    … …  
    50255029      REAL tau
    50265030!c      DATA tau /3600./
    5027 !      DATA tau /5400./
    5028        DATA tau /43200./
     5031      DATA tau /5400./
    50295032!
    50305033   INTEGER k,i
    … …  
    50385041        DO k = 1,klev
    50395042         DO i = 1,klon
    5040 !CR: nudging everywhere
    5041 !           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
     5043           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
    50425044!
    50435045            d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
    … …  
    50465048            print *,' k,u,d_u,v,d_v ',    &
    50475049                      k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
    5048 !           ENDIF
     5050           ENDIF
    50495051!
    50505052         ENDDO
    … …  
    51735175         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
    51745176         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
    5175          dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))
    51765177     
    51775178         else !play>plev_prof_cas(1)
    … …  
    52105211         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
    52115212         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
    5212          dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)
    52135213
    52145214         endif ! play.le.plev_prof_cas(1)
    … …  
    52475247         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                     !jyg
    52485248         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                     !jyg
    5249          dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact               !jyg
    52505249 
    52515250        endif ! play

  • 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         

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

    r2716 r2920  
    7676        dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
    7777        do k=2,kmax
     78          print *,'k l height(k) height(k-1) zlay(l) frac=',k,l,height(k),height(k-1),zlay(l),frac
    7879          frac = (height(k)-zlay(l))/(height(k)-height(k-1))
    7980          if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
    … …  
    227228!?       omega2(l)=-rho(l)*omega(l)
    228229       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    229        d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
     230       d_t_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
    230231       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
    231232       d_q_adv(l,2) = 0.0
    … …  
    280281!on applique le forcage total au premier pas de temps
    281282!attention: signe different de toga
    282        d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
     283       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
    283284       d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
    284285       d_q_adv(l,2) = 0.0
    … …  
    341342!on applique le forcage total au premier pas de temps
    342343!attention: signe different de toga
    343        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
     344       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
    344345!forcage en th
    345 !       d_th_adv(l) = ht_mod(l)
     346!       d_t_adv(l) = ht_mod(l)
    346347       d_q_adv(l,1) = hq_mod(l)
    347348       d_q_adv(l,2) = 0.0
    … …  
    442443!on applique le forcage total au premier pas de temps
    443444!attention: signe different de toga
    444        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
     445       d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
    445446!forcage en th
    446 !       d_th_adv(l) = ht_mod(l)
     447!       d_t_adv(l) = ht_mod(l)
    447448       d_q_adv(l,1) = hq_mod(l)
    448449       d_q_adv(l,2) = 0.0
    … …  
    557558!on applique le forcage total au premier pas de temps
    558559!attention: signe different de toga
    559 !      d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
     560!      d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)
    560561!forcage en th
    561        d_th_adv(l) = ht_mod(l)
     562       d_t_adv(l) = ht_mod(l)
    562563       d_q_adv(l,1) = hq_mod(l)
    563564       d_q_adv(l,2) = 0.0
    … …  
    657658! Advective forcings are given in K or g/kg ... per HOUR
    658659!      IF(height(l).LT.1000) THEN
    659 !        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
     660!        d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
    660661!        d_q_adv(l,1) = adv_qt_prof/1000./3600.
    661662!        d_q_adv(l,2) = 0.0
    662663!      ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
    663 !        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*
     664!        d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*
    664665!    :               (1-(height(l)-1000.)/2000.)
    665 !        d_th_adv(l) = d_th_adv(l)/3600.
     666!        d_t_adv(l) = d_t_adv(l)/3600.
    666667!        d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
    667668!        d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
    668669!        d_q_adv(l,2) = 0.0
    669670!      ELSE
    670 !        d_th_adv(l) = 0.0
     671!        d_t_adv(l) = 0.0
    671672!        d_q_adv(l,1) = 0.0
    672673!        d_q_adv(l,2) = 0.0
    … …  
    751752!?       omega2(l)=-rho(l)*omega(l)
    752753       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    753 !      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
     754!      d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
    754755!      d_q_adv(l,1) = vq_mod(l)
    755        d_th_adv(l) = alpha*omega(l)/rcpd
     756       d_t_adv(l) = alpha*omega(l)/rcpd
    756757       d_q_adv(l,1) = 0.0
    757758       d_q_adv(l,2) = 0.0
    … …  
    827828!      omega2(l)=-rho(l)*omega(l)
    828829       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
    829 !      d_th_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
     830!      d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)
    830831!      d_q_adv(l,1) = vq_mod(l)
    831        d_th_adv(l) = alpha*omega(l)/rcpd
     832       d_t_adv(l) = alpha*omega(l)/rcpd
    832833       d_q_adv(l,1) = 0.0
    833834       d_q_adv(l,2) = 0.0
    … …  
    890891!on applique le forcage total au premier pas de temps
    891892!attention: signe different de toga
    892        d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
     893       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
    893894       d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
    894895       d_q_adv(l,2) = 0.0
    895896       d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
    896        d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
     897! correction bug d_u -> d_v (MM+MPL 20170310)
     898       d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
    897899      enddo     
    898900
    … …  
    971973       q(l,2) = ql_mod_cas(l)
    972974       u(l) = u_mod_cas(l)
    973        ug(l)= u_mod_cas(l)
     975       ug(l)= ug_mod_cas(l)
    974976       v(l) = v_mod_cas(l)
    975        vg(l)= v_mod_cas(l)
    976        omega(l) = w_mod_cas(l)
     977       vg(l)= vg_mod_cas(l)
     978! Modif w_mod_cas -> omega_mod_cas (MM+MPL 20170309)
     979       omega(l) = omega_mod_cas(l)
    977980       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
    978981
    … …  
    980983!on applique le forcage total au premier pas de temps
    981984!attention: signe different de toga
    982        d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
     985       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
    983986       d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
    984987!      d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
    … …  
    987990!      d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
    988991       d_u_adv(l) = du_mod_cas(l)
    989 !      d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
    990        d_u_adv(l) = dv_mod_cas(l)
     992!      d_v_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
     993! correction bug d_u -> d_v (MM+MPL 20170310)
     994       d_v_adv(l) = dv_mod_cas(l)
    991995      enddo     
    992996

  • LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90

    r2904 r2920  
    2020       wake_deltaq, wake_deltat, wake_s, wake_dens, &
    2121       zgam, zmax0, zmea, zpic, zsig, &
    22        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, ql_ancien, qs_ancien, &
    23        prlw_ancien, prsw_ancien, prw_ancien
     22       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
    2423 
    2524   USE dimphy
    … …  
    195194      real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
    196195      real :: dt_dyn(llm)
    197       real :: dt_cooling(llm),d_t_adv(llm),d_th_adv(llm),d_t_nudge(llm)
     196      real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
    198197      real :: d_u_nudge(llm),d_v_nudge(llm)
    199198      real :: du_adv(llm),dv_adv(llm)
    … …  
    207206      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
    208207      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
    209 !      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
     208!      REAL, ALLOCATABLE, DIMENSION(:):: d_t_adv
    210209
    211210!---------------------------------------------------------------------
    … …  
    272271      dt_dyn(:)=0.
    273272      dt_cooling(:)=0.
    274       d_th_adv(:)=0.
     273      d_t_adv(:)=0.
    275274      d_t_nudge(:)=0.
    276275      d_u_nudge(:)=0.
    … …  
    405404       heure_ini_cas=0.
    406405       pdt_cas=1800.         ! forcing frequency
     406      elseif (forcing_type .eq.105) THEN ! bomex starts 16-12-2004 0h
     407       forcing_case2 = .true.
     408       year_ini_cas=1969
     409       mth_ini_cas=6
     410       day_deb=24
     411       heure_ini_cas=0.
     412       pdt_cas=1800.         ! forcing frequency
    407413      elseif (forcing_type .eq.40) THEN
    408414       forcing_GCSSold = .true.
    … …  
    574580      allocate(d_q_adv(llm,nqtot))
    575581      allocate(d_q_nudge(llm,nqtot))
    576 !      allocate(d_th_adv(llm))
     582!      allocate(d_t_adv(llm))
    577583
    578584      q(:,:) = 0.
    … …  
    807813        t_ancien(1,:)=temp(:)
    808814        q_ancien(1,:)=q(:,1)
    809         ql_ancien = 0.
    810         qs_ancien = 0.
    811         prlw_ancien = 0.
    812         prsw_ancien = 0.
    813         prw_ancien = 0.
    814815!jyg<
    815816!!        pbl_tke(:,:,:)=1.e-8
    … …  
    10631064         fcoriolis=0.0
    10641065         dt_cooling=0.0
    1065          d_th_adv=0.0
     1066         d_t_adv=0.0
    10661067         d_q_adv=0.0
    10671068       endif
    … …  
    10761077          dt_cooling=0.
    10771078       endif
    1078 
    1079 !CRio:Attention modif spécifique cas de Caroline
    1080       if (forcing_les) then
    1081          fcoriolis=0.
    1082 !Nudging
    1083        
    1084 !on calcule dt_cooling
    1085         do l=1,llm
    1086         if (play(l).ge.20000.) then
    1087             dt_cooling(l)=-1.5/86400.
    1088         elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then
    1089             dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)
    1090         else
    1091             dt_cooling(l)=-1.*(temp(l)-200.)/86400.
    1092         endif
    1093         enddo
    1094 
    1095       endif     
    1096 !RC
    10971079
    10981080      IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
    … …  
    11681150        if (prt_level.ge.3) then
    11691151          print *,                                                          &
    1170      &    'physiq-> temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1) ',         &
    1171      &              temp(1),dt_phys(1),d_th_adv(1),dt_cooling(1)
     1152     &    'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ',         &
     1153     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
    11721154           print* ,'dv_phys=',dv_phys
    11731155           print* ,'dv_age=',dv_age
    … …  
    11801162        temp(1:mxcalc)=temp(1:mxcalc)+timestep*(                            &
    11811163     &              dt_phys(1:mxcalc)                                       &
    1182      &             +d_th_adv(1:mxcalc)                                      &
     1164     &             +d_t_adv(1:mxcalc)                                      &
    11831165     &             +d_t_nudge(1:mxcalc)                                      &
    11841166     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
Note: See TracChangeset for help on using the changeset viewer.