Changeset 2574
- Timestamp:
- Jun 15, 2016, 1:26:59 PM (8 years ago)
- Location:
- LMDZ5/trunk/libf/phylmd
- Files:
-
- 2 edited
-
phys_local_var_mod.F90 (modified) (1 diff)
-
yamada4.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
r2561 r2574 408 408 allocate(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev)) 409 409 allocate(u_seri(klon,klev),v_seri(klon,klev)) 410 allocate(l_mixmin(klon,klev,nbsrf), l_mix(klon,klev,nbtr)) 410 allocate(l_mixmin(klon,klev,nbsrf), l_mix(klon,klev,nbsrf)) 411 l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ! doit etre initialse car pas toujours remplis 411 412 412 413 allocate(tr_seri(klon,klev,nbtr)) -
LMDZ5/trunk/libf/phylmd/yamada4.F90
r2573 r2574 212 212 ! initialize arrays: 213 213 214 m2( :, :) = 0.0215 sm( :, :) = 0.0216 rif( :, :) = 0.0214 m2(1:ngrid, :) = 0.0 215 sm(1:ngrid, :) = 0.0 216 rif(1:ngrid, :) = 0.0 217 217 218 218 !------------------------------------------------------------ … … 264 264 265 265 DO k = 2, klev 266 q2( :, k) = l(:, k)**2*zz(:, k)266 q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k) 267 267 END DO 268 268 … … 349 349 350 350 DO k = 2, klev - 1 351 km( :, k) = l(:, k)*sqrt(q2(:,k))*sm(:, k)352 q2( :, k) = q2(:, k) + dt*km(:, k)*m2(:, k)*(1.-rif(:,k))353 q2( :, k) = min(max(q2(:,k),1.E-10), 1.E4)354 q2( :, k) = 1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1))355 q2( :, k) = q2(:, k)*q2(:, k)351 km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) 352 q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) 353 q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) 354 q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) 355 q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) 356 356 END DO 357 357 … … 392 392 393 393 IF (iflag_pbl>=12) THEN 394 q2( :, 1) = q2(:, 2)394 q2(1:ngrid, 1) = q2(1:ngrid, 2) 395 395 CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2) 396 396 END IF … … 803 803 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1. 804 804 805 806 807 805 IF (ngrid==0) RETURN 808 806 809 807 IF (firstcall) THEN … … 811 809 firstcall = .FALSE. 812 810 END IF 813 814 811 815 812 … … 853 850 !...................................................................... 854 851 855 l0( :) = 150.852 l0(1:ngrid) = 150. 856 853 DO k = 2, klev 857 854 DO ig = 1, ngrid … … 868 865 !================================================================================= 869 866 870 l2( :,:)=0.0871 l_mixmin( :,:,nsrf)=0.872 l_mix( :,:,nsrf)=0.867 l2(1:ngrid,:)=0.0 868 l_mixmin(1:ngrid,:,nsrf)=0. 869 l_mix(1:ngrid,:,nsrf)=0. 873 870 874 871 IF (nsrf .EQ. 1) THEN … … 877 874 !-------------- 878 875 879 extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1.880 lmax=150. ! Longueur de m??lange max dans l'absolu876 extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. 877 lmax=150. ! Longueur de m??lange max dans l'absolu 881 878 882 879 ! calculs 883 880 !--------- 884 881 885 DO ig=1,ngrid882 DO ig=1,ngrid 886 883 887 884 ! On calcule la hauteur du relief 888 885 !................................. 889 890 886 ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille 891 887 ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon) 892 888 ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief 893 889 ! (en gros, une maille de taille xcell avec une pente constante zstdslope) 894 895 896 890 jg=ni(ig) 897 898 891 ! IF (zsig(jg) .EQ. 0.) THEN 899 892 ! zstdslope(ig)=0. … … 907 900 h_oro(ig)=zstd(jg) 908 901 hlim(ig)=extent*h_oro(ig) 909 910 911 912 913 END DO 914 915 l2limit(1:ngrid)=0. 916 917 DO k=2,klev 918 DO ig=1,ngrid 919 920 winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) 921 922 IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! Si on est sous l'orographie 923 924 925 l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008) 926 927 l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar qui tend asymptotiquement vers h 928 l2neutre=MIN(l2neutre,lmax) ! On majore par lmax 929 l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) 930 931 l2(ig,k)=l2limit(ig) 902 ENDDO 903 904 l2limit(1:ngrid)=0. 905 906 DO k=2,klev 907 DO ig=1,ngrid 908 winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) 909 IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie 910 l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008) 911 l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h 912 l2neutre=MIN(l2neutre,lmax) ! On majore par lmax 913 l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) 914 l2(ig,k)=l2limit(ig) 932 915 933 ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles916 ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles 934 917 935 918 ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes 936 919 ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z) 937 920 ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim 938 939 l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) 940 941 942 ELSE ! Au dessus de extent*h, on prend l2=l0 943 944 l2(ig,k)=0. 945 946 947 END IF 948 949 END DO 950 END DO 951 952 953 END IF ! pbl_lmixmin_alpha 921 l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) 922 ELSE ! Au dessus de extent*h, on prend l2=l0 923 l2(ig,k)=0. 924 END IF 925 ENDDO 926 ENDDO 927 ENDIF ! pbl_lmixmin_alpha 954 928 955 929 !================================================================================== … … 961 935 DO k=2,klev 962 936 DO ig=1,ngrid 963 964 lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)965 966 END DO 967 END DO 937 lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin) 938 ENDDO 939 ENDDO 940 941 ! Diagnostics 968 942 969 943 DO k=2,klev 970 DO ig=1,ngrid971 jg=ni(ig)972 l_mix(jg,k,nsrf)=lmix(ig,k)973 l_mixmin(jg,k,nsrf)=l2(ig,k)974 ENDDO944 DO ig=1,ngrid 945 jg=ni(ig) 946 l_mix(jg,k,nsrf)=lmix(ig,k) 947 l_mixmin(jg,k,nsrf)=l2(ig,k) 948 ENDDO 975 949 ENDDO 976 977 jg=ni(ig)978 979 950 DO ig=1,ngrid 951 jg=ni(ig) 952 l_mix(jg,1,nsrf)=hlim(ig) 953 ENDDO 980 954 981 955
Note: See TracChangeset
for help on using the changeset viewer.
