Changeset 2920
- Timestamp:
- Jun 29, 2017, 11:58:07 AM (7 years ago)
- Location:
- LMDZ5/trunk/libf/phylmd/dyn1d
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/dyn1d/1DUTILS.h
r2904 r2920 482 482 forc_omega =0 483 483 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) 484 498 485 499 !Config Key = forc_w … … 653 667 654 668 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 663 674 print*,'in dyn1deta0 ',fichnom,klon,klev,nqtot 664 675 CALL open_startphy(fichnom) … … 807 818 CALL open_restartphy(fichnom) 808 819 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" 816 824 817 825 modname = 'dyn1dredem' … … 1396 1404 cor(:) = rkappa*temp*(1.+q(:,1)*rv/rd)/(play*(1.+q(:,1))) 1397 1405 1398 1399 1406 do k=2,llm-1 1400 1407 … … 2789 2796 hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1)) 2790 2797 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))2792 2798 2793 2799 else !play>plev_prof_cas(1) … … 2816 2822 hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2) 2817 2823 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)2819 2824 2820 2825 endif ! play.le.plev_prof_cas(1) … … 2845 2850 hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg 2846 2851 vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg 2847 dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg2848 2852 2849 2853 endif ! play … … 5025 5029 REAL tau 5026 5030 !c DATA tau /3600./ 5027 ! DATA tau /5400./ 5028 DATA tau /43200./ 5031 DATA tau /5400./ 5029 5032 ! 5030 5033 INTEGER k,i … … 5038 5041 DO k = 1,klev 5039 5042 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 5042 5044 ! 5043 5045 d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k)) … … 5046 5048 print *,' k,u,d_u,v,d_v ', & 5047 5049 k,u(i,k),d_u(i,k),v(i,k),d_v(i,k) 5048 !ENDIF5050 ENDIF 5049 5051 ! 5050 5052 ENDDO … … 5173 5175 hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1)) 5174 5176 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))5176 5177 5177 5178 else !play>plev_prof_cas(1) … … 5210 5211 hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2) 5211 5212 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)5213 5213 5214 5214 endif ! play.le.plev_prof_cas(1) … … 5247 5247 hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg 5248 5248 vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg 5249 dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg5250 5249 5251 5250 endif ! play -
LMDZ5/trunk/libf/phylmd/dyn1d/1D_interp_cases.h
r2716 r2920 37 37 38 38 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 39 d_t h_adv(l) = ht_gcssold(l)39 d_t_adv(l) = ht_gcssold(l) 40 40 d_q_adv(l,1) = hq_gcssold(l) 41 41 dt_cooling(l) = 0.0 … … 84 84 85 85 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 86 d_t h_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)) 87 87 d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l)) 88 88 dt_cooling(l) = 0.0 … … 183 183 184 184 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 185 d_t h_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) 186 186 d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l) 187 187 d_u_adv(l) = hu_mod(l)-d_u_dyn_z(l) … … 224 224 ug(l)= ug_mod(l) 225 225 vg(l)= vg_mod(l) 226 d_t h_adv(l)=ht_mod(l)226 d_t_adv(l)=ht_mod(l) 227 227 d_q_adv(l,1)=hq_mod(l) 228 228 enddo … … 315 315 !calcul de l'advection totale 316 316 if (cptadvw) then 317 d_t h_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) 318 318 ! print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l) 319 319 d_q_adv(l,1) = hq_mod(l)-d_q_dyn_z(l) 320 320 ! print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l) 321 321 else 322 d_t h_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)) 323 323 d_q_adv(l,1) = (hq_mod(l)+vq_mod(l)) 324 324 endif … … 392 392 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 393 393 !calcul de l'advection totale 394 ! d_t h_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) 395 395 !attention: on impose dth 396 d_t h_adv(l) = alpha*omega(l)/rcpd+ &396 d_t_adv(l) = alpha*omega(l)/rcpd+ & 397 397 & ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l) 398 ! d_t h_adv(l) = 0.398 ! d_t_adv(l) = 0. 399 399 ! print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l) 400 400 d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l) … … 417 417 !--------------------------------------------------------------------- 418 418 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) 421 420 call lstendH(llm,nqtot,omega,dt_dyn,dq_dyn,q,temp,u,v,play) 422 421 423 422 do l=1,llm 424 d_t h_adv(l) = (dth_rico(l) + dt_dyn(l))423 d_t_adv(l) = (dth_rico(l) + dt_dyn(l)) 425 424 d_q_adv(l,1) = (dqh_rico(l) + dq_dyn(l,1)) 426 425 d_q_adv(l,2) = 0. … … 458 457 vg(l)= v_mod(l) 459 458 IF((phi(l)/RG).LT.1000) THEN 460 d_t h_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.459 d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600. 461 460 d_q_adv(l,1) = adv_qt_prof/1000./3600. 462 461 d_q_adv(l,2) = 0.0 463 462 ! print *,'INF1000: phi dth dq1 dq2', 464 ! : phi(l)/RG,d_t h_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) 465 464 ELSEIF ((phi(l)/RG).GE.1000.AND.(phi(l)/RG).lt.3000) THEN 466 465 fact=((phi(l)/RG)-1000.)/2000. 467 466 fact=1-fact 468 d_t h_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600.467 d_t_adv(l) = (adv_theta_prof + rad_theta_prof)*fact/3600. 469 468 d_q_adv(l,1) = adv_qt_prof*fact/1000./3600. 470 469 d_q_adv(l,2) = 0.0 471 470 ! print *,'SUP1000: phi fact dth dq1 dq2', 472 ! : phi(l)/RG,fact,d_t h_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) 473 472 ELSE 474 d_t h_adv(l) = 0.0473 d_t_adv(l) = 0.0 475 474 d_q_adv(l,1) = 0.0 476 475 d_q_adv(l,2) = 0.0 477 476 ! print *,'SUP3000: phi dth dq1 dq2', 478 ! : phi(l)/RG,d_t h_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) 479 478 ENDIF 480 479 dt_cooling(l) = 0.0 481 480 ! print *,'Interp armcu: phi dth dq1 dq2', 482 ! : l,phi(l),d_t h_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) 483 482 enddo 484 483 endif ! forcing_armcu … … 554 553 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 555 554 ! 556 ! d_t h_adv(l) = 0.0555 ! d_t_adv(l) = 0.0 557 556 ! d_q_adv(l,1) = 0.0 558 557 !CR:test advection=0 559 558 !calcul de l'advection verticale 560 d_t h_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) 561 560 ! print*,'temp adv',l,-d_t_dyn_z(l) 562 561 d_q_adv(l,1) = -d_q_dyn_z(l) … … 637 636 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 638 637 ! 639 ! d_t h_adv(l) = 0.0638 ! d_t_adv(l) = 0.0 640 639 ! d_q_adv(l,1) = 0.0 641 640 !CR:test advection=0 642 641 !calcul de l'advection verticale 643 d_t h_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) 644 643 ! print*,'temp adv',l,-d_t_dyn_z(l) 645 644 d_q_adv(l,1) = -d_q_dyn_z(l) … … 715 714 716 715 !Calcul de l advection verticale 716 717 717 d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:) 718 718 719 d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:) 719 720 d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:) … … 764 765 765 766 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) 767 768 omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 768 769 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) … … 782 783 783 784 if ((tend_t.eq.1).and.(tend_w.eq.0)) then 784 ! d_t h_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)785 d_t h_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) 786 787 else if ((tend_t.eq.1).and.(tend_w.eq.1)) then 787 ! d_t h_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)788 d_t h_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) 789 790 endif 790 791 … … 816 817 ENDIF 817 818 endif ! forcing_case 818 819 819 820 820 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 877 877 d_th_z(:)=0. 878 878 d_q_z(:)=0. 879 d_u_z(:)=0. 880 d_v_z(:)=0. 879 881 d_t_dyn_z(:)=0. 880 882 d_th_dyn_z(:)=0. 881 883 d_q_dyn_z(:)=0. 884 d_u_dyn_z(:)=0. 885 d_v_dyn_z(:)=0. 882 886 DO l=2,llm-1 883 887 d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1)) 884 888 d_th_z(l)=(teta(l+1)-teta(l-1))/(play(l+1)-play(l-1)) 885 889 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)) 886 892 ENDDO 887 893 d_t_z(1)=d_t_z(2) 888 894 d_th_z(1)=d_th_z(2) 889 895 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) 890 898 d_t_z(llm)=d_t_z(llm-1) 891 899 d_th_z(llm)=d_th_z(llm-1) 892 900 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) 893 903 894 904 !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 899 919 !wind nudging 900 920 if (nudging_u.gt.0.) then … … 902 922 u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u) 903 923 enddo 904 905 906 ug(l) = u_mod_cas(l)907 924 ! else 925 ! do l=1,llm 926 ! u(l) = u_mod_cas(l) 927 ! enddo 908 928 endif 909 929 … … 912 932 v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v) 913 933 enddo 914 915 916 vg(l) = v_mod_cas(l)917 934 ! else 935 ! do l=1,llm 936 ! v(l) = v_mod_cas(l) 937 ! enddo 918 938 endif 919 939 … … 922 942 w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w) 923 943 enddo 924 925 926 w(l) = w_mod_cas(l)927 944 ! else 945 ! do l=1,llm 946 ! w(l) = w_mod_cas(l) 947 ! enddo 928 948 endif 929 949 … … 941 961 942 962 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) 944 965 omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 945 966 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 946 967 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) 952 997 ! 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 979 1000 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 981 1003 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) 982 1005 endif 983 1006 -
LMDZ5/trunk/libf/phylmd/dyn1d/1D_read_forc_cases.h
r2716 r2920 76 76 dt_cooling(l)=thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1)) 77 77 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 78 79 frac = (height(k)-zlay(l))/(height(k)-height(k-1)) 79 80 if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k) … … 227 228 !? omega2(l)=-rho(l)*omega(l) 228 229 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 229 d_t h_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)) 230 231 d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l)) 231 232 d_q_adv(l,2) = 0.0 … … 280 281 !on applique le forcage total au premier pas de temps 281 282 !attention: signe different de toga 282 d_t h_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)) 283 284 d_q_adv(l,1) = (hq_mod(l)+vq_mod(l)) 284 285 d_q_adv(l,2) = 0.0 … … 341 342 !on applique le forcage total au premier pas de temps 342 343 !attention: signe different de toga 343 d_t h_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)344 d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 344 345 !forcage en th 345 ! d_t h_adv(l) = ht_mod(l)346 ! d_t_adv(l) = ht_mod(l) 346 347 d_q_adv(l,1) = hq_mod(l) 347 348 d_q_adv(l,2) = 0.0 … … 442 443 !on applique le forcage total au premier pas de temps 443 444 !attention: signe different de toga 444 d_t h_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)445 d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 445 446 !forcage en th 446 ! d_t h_adv(l) = ht_mod(l)447 ! d_t_adv(l) = ht_mod(l) 447 448 d_q_adv(l,1) = hq_mod(l) 448 449 d_q_adv(l,2) = 0.0 … … 557 558 !on applique le forcage total au premier pas de temps 558 559 !attention: signe different de toga 559 ! d_t h_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)560 ! d_t_adv(l) = alpha*omega(l)/rcpd+ht_mod(l) 560 561 !forcage en th 561 d_t h_adv(l) = ht_mod(l)562 d_t_adv(l) = ht_mod(l) 562 563 d_q_adv(l,1) = hq_mod(l) 563 564 d_q_adv(l,2) = 0.0 … … 657 658 ! Advective forcings are given in K or g/kg ... per HOUR 658 659 ! IF(height(l).LT.1000) THEN 659 ! d_t h_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.660 ! d_t_adv(l) = (adv_theta_prof + rad_theta_prof)/3600. 660 661 ! d_q_adv(l,1) = adv_qt_prof/1000./3600. 661 662 ! d_q_adv(l,2) = 0.0 662 663 ! ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN 663 ! d_t h_adv(l) = (adv_theta_prof + rad_theta_prof)*664 ! d_t_adv(l) = (adv_theta_prof + rad_theta_prof)* 664 665 ! : (1-(height(l)-1000.)/2000.) 665 ! d_t h_adv(l) = d_th_adv(l)/3600.666 ! d_t_adv(l) = d_t_adv(l)/3600. 666 667 ! d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.) 667 668 ! d_q_adv(l,1) = d_q_adv(l,1)/1000./3600. 668 669 ! d_q_adv(l,2) = 0.0 669 670 ! ELSE 670 ! d_t h_adv(l) = 0.0671 ! d_t_adv(l) = 0.0 671 672 ! d_q_adv(l,1) = 0.0 672 673 ! d_q_adv(l,2) = 0.0 … … 751 752 !? omega2(l)=-rho(l)*omega(l) 752 753 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 753 ! d_t h_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)754 ! d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) 754 755 ! d_q_adv(l,1) = vq_mod(l) 755 d_t h_adv(l) = alpha*omega(l)/rcpd756 d_t_adv(l) = alpha*omega(l)/rcpd 756 757 d_q_adv(l,1) = 0.0 757 758 d_q_adv(l,2) = 0.0 … … 827 828 ! omega2(l)=-rho(l)*omega(l) 828 829 alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l) 829 ! d_t h_adv(l) = alpha*omega(l)/rcpd+vt_mod(l)830 ! d_t_adv(l) = alpha*omega(l)/rcpd+vt_mod(l) 830 831 ! d_q_adv(l,1) = vq_mod(l) 831 d_t h_adv(l) = alpha*omega(l)/rcpd832 d_t_adv(l) = alpha*omega(l)/rcpd 832 833 d_q_adv(l,1) = 0.0 833 834 d_q_adv(l,2) = 0.0 … … 890 891 !on applique le forcage total au premier pas de temps 891 892 !attention: signe different de toga 892 d_t h_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)) 893 894 d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l)) 894 895 d_q_adv(l,2) = 0.0 895 896 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)) 897 899 enddo 898 900 … … 971 973 q(l,2) = ql_mod_cas(l) 972 974 u(l) = u_mod_cas(l) 973 ug(l)= u _mod_cas(l)975 ug(l)= ug_mod_cas(l) 974 976 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) 977 980 omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq 978 981 … … 980 983 !on applique le forcage total au premier pas de temps 981 984 !attention: signe different de toga 982 d_t h_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)) 983 986 d_t_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l)) 984 987 ! d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l)) … … 987 990 ! d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l)) 988 991 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) 991 995 enddo 992 996 -
LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90
r2904 r2920 20 20 wake_deltaq, wake_deltat, wake_s, wake_dens, & 21 21 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 24 23 25 24 USE dimphy … … 195 194 real :: du_phys(llm),dv_phys(llm),dt_phys(llm) 196 195 real :: dt_dyn(llm) 197 real :: dt_cooling(llm),d_t_adv(llm),d_t h_adv(llm),d_t_nudge(llm)196 real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm) 198 197 real :: d_u_nudge(llm),d_v_nudge(llm) 199 198 real :: du_adv(llm),dv_adv(llm) … … 207 206 REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv 208 207 REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge 209 ! REAL, ALLOCATABLE, DIMENSION(:):: d_t h_adv208 ! REAL, ALLOCATABLE, DIMENSION(:):: d_t_adv 210 209 211 210 !--------------------------------------------------------------------- … … 272 271 dt_dyn(:)=0. 273 272 dt_cooling(:)=0. 274 d_t h_adv(:)=0.273 d_t_adv(:)=0. 275 274 d_t_nudge(:)=0. 276 275 d_u_nudge(:)=0. … … 405 404 heure_ini_cas=0. 406 405 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 407 413 elseif (forcing_type .eq.40) THEN 408 414 forcing_GCSSold = .true. … … 574 580 allocate(d_q_adv(llm,nqtot)) 575 581 allocate(d_q_nudge(llm,nqtot)) 576 ! allocate(d_t h_adv(llm))582 ! allocate(d_t_adv(llm)) 577 583 578 584 q(:,:) = 0. … … 807 813 t_ancien(1,:)=temp(:) 808 814 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.814 815 !jyg< 815 816 !! pbl_tke(:,:,:)=1.e-8 … … 1063 1064 fcoriolis=0.0 1064 1065 dt_cooling=0.0 1065 d_t h_adv=0.01066 d_t_adv=0.0 1066 1067 d_q_adv=0.0 1067 1068 endif … … 1076 1077 dt_cooling=0. 1077 1078 endif 1078 1079 !CRio:Attention modif spécifique cas de Caroline1080 if (forcing_les) then1081 fcoriolis=0.1082 !Nudging1083 1084 !on calcule dt_cooling1085 do l=1,llm1086 if (play(l).ge.20000.) then1087 dt_cooling(l)=-1.5/86400.1088 elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then1089 dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)1090 else1091 dt_cooling(l)=-1.*(temp(l)-200.)/86400.1092 endif1093 enddo1094 1095 endif1096 !RC1097 1079 1098 1080 IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', & … … 1168 1150 if (prt_level.ge.3) then 1169 1151 print *, & 1170 & 'physiq-> temp(1),dt_phys(1),d_t h_adv(1),dt_cooling(1) ', &1171 & temp(1),dt_phys(1),d_t h_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) 1172 1154 print* ,'dv_phys=',dv_phys 1173 1155 print* ,'dv_age=',dv_age … … 1180 1162 temp(1:mxcalc)=temp(1:mxcalc)+timestep*( & 1181 1163 & dt_phys(1:mxcalc) & 1182 & +d_t h_adv(1:mxcalc) &1164 & +d_t_adv(1:mxcalc) & 1183 1165 & +d_t_nudge(1:mxcalc) & 1184 1166 & +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid.
Note: See TracChangeset
for help on using the changeset viewer.