Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 2573) +++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 2574) @@ -408,5 +408,6 @@ allocate(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev)) allocate(u_seri(klon,klev),v_seri(klon,klev)) - allocate(l_mixmin(klon,klev,nbsrf), l_mix(klon,klev,nbtr)) + allocate(l_mixmin(klon,klev,nbsrf), l_mix(klon,klev,nbsrf)) + l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ! doit etre initialse car pas toujours remplis allocate(tr_seri(klon,klev,nbtr)) Index: LMDZ5/trunk/libf/phylmd/yamada4.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/yamada4.F90 (revision 2573) +++ LMDZ5/trunk/libf/phylmd/yamada4.F90 (revision 2574) @@ -212,7 +212,7 @@ ! initialize arrays: - m2(:, :) = 0.0 - sm(:, :) = 0.0 - rif(:, :) = 0.0 + m2(1:ngrid, :) = 0.0 + sm(1:ngrid, :) = 0.0 + rif(1:ngrid, :) = 0.0 !------------------------------------------------------------ @@ -264,5 +264,5 @@ DO k = 2, klev - q2(:, k) = l(:, k)**2*zz(:, k) + q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k) END DO @@ -349,9 +349,9 @@ DO k = 2, klev - 1 - km(:, k) = l(:, k)*sqrt(q2(:,k))*sm(:, k) - q2(:, k) = q2(:, k) + dt*km(:, k)*m2(:, k)*(1.-rif(:,k)) - q2(:, k) = min(max(q2(:,k),1.E-10), 1.E4) - q2(:, k) = 1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1)) - q2(:, k) = q2(:, k)*q2(:, k) + km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) + q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) + q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) + q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) + q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) END DO @@ -392,5 +392,5 @@ IF (iflag_pbl>=12) THEN - q2(:, 1) = q2(:, 2) + q2(1:ngrid, 1) = q2(1:ngrid, 2) CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2) END IF @@ -803,7 +803,5 @@ famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1. - - - + IF (ngrid==0) RETURN IF (firstcall) THEN @@ -811,5 +809,4 @@ firstcall = .FALSE. END IF - @@ -853,5 +850,5 @@ !...................................................................... - l0(:) = 150. + l0(1:ngrid) = 150. DO k = 2, klev DO ig = 1, ngrid @@ -868,7 +865,7 @@ !================================================================================= - l2(:,:)=0.0 - l_mixmin(:,:,nsrf)=0. - l_mix(:,:,nsrf)=0. + l2(1:ngrid,:)=0.0 + l_mixmin(1:ngrid,:,nsrf)=0. + l_mix(1:ngrid,:,nsrf)=0. IF (nsrf .EQ. 1) THEN @@ -877,23 +874,19 @@ !-------------- - extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. - lmax=150. ! Longueur de m??lange max dans l'absolu + extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. + lmax=150. ! Longueur de m??lange max dans l'absolu ! calculs !--------- - DO ig=1,ngrid + DO ig=1,ngrid ! On calcule la hauteur du relief !................................. - ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon) ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief ! (en gros, une maille de taille xcell avec une pente constante zstdslope) - - jg=ni(ig) - ! IF (zsig(jg) .EQ. 0.) THEN ! zstdslope(ig)=0. @@ -907,49 +900,30 @@ h_oro(ig)=zstd(jg) hlim(ig)=extent*h_oro(ig) - - - - - END DO - - l2limit(1:ngrid)=0. - - DO k=2,klev - DO ig=1,ngrid - - winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) - - IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! Si on est sous l'orographie - - - 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) - - 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 - l2neutre=MIN(l2neutre,lmax) ! On majore par lmax - l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) - - l2(ig,k)=l2limit(ig) + ENDDO + + l2limit(1:ngrid)=0. + + DO k=2,klev + DO ig=1,ngrid + winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) + IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie + 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) + 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 + l2neutre=MIN(l2neutre,lmax) ! On majore par lmax + l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) + l2(ig,k)=l2limit(ig) - ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles + ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z) ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim - - l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) - - - ELSE ! Au dessus de extent*h, on prend l2=l0 - - l2(ig,k)=0. - - - END IF - - END DO - END DO - - - END IF ! pbl_lmixmin_alpha + l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) + ELSE ! Au dessus de extent*h, on prend l2=l0 + l2(ig,k)=0. + END IF + ENDDO + ENDDO + ENDIF ! pbl_lmixmin_alpha !================================================================================== @@ -961,21 +935,21 @@ DO k=2,klev DO ig=1,ngrid - - lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin) - - END DO - END DO + lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin) + ENDDO + ENDDO + +! Diagnostics DO k=2,klev - DO ig=1,ngrid - jg=ni(ig) - l_mix(jg,k,nsrf)=lmix(ig,k) - l_mixmin(jg,k,nsrf)=l2(ig,k) - ENDDO + DO ig=1,ngrid + jg=ni(ig) + l_mix(jg,k,nsrf)=lmix(ig,k) + l_mixmin(jg,k,nsrf)=l2(ig,k) + ENDDO ENDDO - DO ig=1,ngrid - jg=ni(ig) - l_mix(jg,1,nsrf)=hlim(ig) - ENDDO + DO ig=1,ngrid + jg=ni(ig) + l_mix(jg,1,nsrf)=hlim(ig) + ENDDO