Index: dynamico_lmdz/simple_physics/phyparam/param/multipl.F =================================================================== --- dynamico_lmdz/simple_physics/phyparam/param/multipl.F (revision 4187) +++ (revision ) @@ -1,17 +1,0 @@ - SUBROUTINE multipl(n,x1,x2,y) - IMPLICIT NONE -c======================================================================= -c -c multiplication de deux vecteurs -c -c======================================================================= -c - INTEGER n,i - REAL x1(n),x2(n),y(n) -c - DO 10 i=1,n - y(i)=x1(i)*x2(i) -10 CONTINUE -c - RETURN - END Index: dynamico_lmdz/simple_physics/phyparam/param/phyparam.F =================================================================== --- dynamico_lmdz/simple_physics/phyparam/param/phyparam.F (revision 4187) +++ dynamico_lmdz/simple_physics/phyparam/param/phyparam.F (revision 4188) @@ -12,4 +12,5 @@ USE phys_const USE planet + USE vdif_mod, ONLY : vdif c IMPLICIT NONE @@ -161,5 +162,5 @@ - EXTERNAL vdif,convadj + EXTERNAL convadj EXTERNAL orbite,mucorr EXTERNAL ismin,ismax Index: dynamico_lmdz/simple_physics/phyparam/param/vdif.F =================================================================== --- dynamico_lmdz/simple_physics/phyparam/param/vdif.F (revision 4187) +++ (revision ) @@ -1,299 +1,0 @@ - subroutiNE vdif(ngrid,nlay,ptime, - $ ptimestep,pcapcal,pz0, - $ pplay,pplev,pzlay,pzlev, - $ pu,pv,ph,ptsrf,pemis, - $ pdufi,pdvfi,pdhfi,pfluxsrf, - $ pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, - $ lwrite) - USE phys_const - USE vdif_mod - IMPLICIT NONE - -c======================================================================= -c -c Diffusion verticale -c Shema implicite -c On commence par rajouter au variables x la tendance physique -c et on resoult en fait: -c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) -c -c arguments: -c ---------- -c -c entree: -c ------- -c -c -c======================================================================= - -c----------------------------------------------------------------------- -c declarations: -c ------------- - -#include "description.h" -c -c arguments: -c ---------- - - INTEGER ngrid,nlay - REAL ptime,ptimestep - REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) - REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1) - REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay) - REAL ptsrf(ngrid),pemis(ngrid) - REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay) - REAL pfluxsrf(ngrid) - REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay) - REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid) - REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1) - LOGICAL lwrite -c -c local: -c ------ - - INTEGER ilev,ig,ilay,nlev - INTEGER unit,ierr,it1,it2 - INTEGER cluvdb,putdat,putvdim,setname,setvdim - REAL z4st,zdplanck(ngrid),zu2 - REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1) - REAL zcdv(ngrid),zcdh(ngrid) - REAL zu(ngrid,nlay),zv(ngrid,nlay) - REAL zh(ngrid,nlay) - REAL ztsrf2(ngrid) - REAL z1(ngrid),z2(ngrid) - REAL za(ngrid,nlay),zb(ngrid,nlay) - REAL zb0(ngrid,nlay) - REAL zc(ngrid,nlay),zd(ngrid,nlay) - REAL zcst1 - - EXTERNAL coefdifv - INTEGER, SAVE :: icount=0 -!$OMP THREADPRIVATE(icount) - -c -c----------------------------------------------------------------------- -c initialisations: -c ---------------- - - nlev=nlay+1 - -c computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp: -c with rho=p/RT=p/ (R Theta) (p/ps)**kappa -c --------------------------------- - - DO ilay=1,nlay - DO ig=1,ngrid - za(ig,ilay)= - s (pplev(ig,ilay)-pplev(ig,ilay+1))/g - ENDDO - ENDDO - - zcst1=4.*g*ptimestep/(r*r) - DO ilev=2,nlev-1 - DO ig=1,ngrid - zb0(ig,ilev)=pplev(ig,ilev)* - s (pplev(ig,1)/pplev(ig,ilev))**rcp / - s (ph(ig,ilev-1)+ph(ig,ilev)) - zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/ - s (pplay(ig,ilev-1)-pplay(ig,ilev)) - ENDDO - ENDDO - DO ig=1,ngrid - zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) - ENDDO - IF(lwrite) THEN - ig=ngrid/2+1 - PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz' - DO ilay=1,nlay - WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay), - s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay) - ENDDO - PRINT*,'Pression (mbar) ,altitude (km),zb' - DO ilev=1,nlay - WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev), - s zb0(ig,ilev) - ENDDO - ENDIF - -c----------------------------------------------------------------------- -c 2. ajout des tendances physiques: -c ------------------------------ - - DO ilev=1,nlay - DO ig=1,ngrid - zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep - zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep - zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c 3. calcul de cd : -c ---------------- -c - CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh) - CALL vdif_k(ngrid,nlay, - s ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh) - - DO ig=1,ngrid - zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) - zcdv(ig)=zcdv(ig)*sqrt(zu2) - zcdh(ig)=zcdh(ig)*sqrt(zu2) - ENDDO - - IF(lwrite) THEN - PRINT* - PRINT*,'Diagnostique diffusion verticale' - PRINT*,'coefficients Cd pour v et h' - PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1) - PRINT*,'coefficients K pour v et h' - DO ilev=1,nlay - PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev) - ENDDO - ENDIF - -c----------------------------------------------------------------------- -c integration verticale pour u: -c ----------------------------- -c - CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2)) - CALL multipl(ngrid,zcdv,zb0,zb) - DO ig=1,ngrid - z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) - zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig) - zd(ig,nlay)=zb(ig,nlay)*z1(ig) - ENDDO - - DO ilay=nlay-1,1,-1 - DO ig=1,ngrid - z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ - $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) - zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ - $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) - zd(ig,ilay)=zb(ig,ilay)*z1(ig) - ENDDO - ENDDO - - DO ig=1,ngrid - zu(ig,1)=zc(ig,1) - ENDDO - DO ilay=2,nlay - DO ig=1,ngrid - zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c integration verticale pour v: -c ----------------------------- -c - DO ig=1,ngrid - z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) - zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig) - zd(ig,nlay)=zb(ig,nlay)*z1(ig) - ENDDO - - DO ilay=nlay-1,1,-1 - DO ig=1,ngrid - z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ - $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) - zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ - $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) - zd(ig,ilay)=zb(ig,ilay)*z1(ig) - ENDDO - ENDDO - - DO ig=1,ngrid - zv(ig,1)=zc(ig,1) - ENDDO - DO ilay=2,nlay - DO ig=1,ngrid - zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c integration verticale pour h: -c ----------------------------- -c - CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) - CALL multipl(ngrid,zcdh,zb0,zb) - DO ig=1,ngrid - z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) - zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig) - zd(ig,nlay)=zb(ig,nlay)*z1(ig) - ENDDO - - DO ilay=nlay-1,1,-1 - DO ig=1,ngrid - z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ - $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) - zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ - $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) - zd(ig,ilay)=zb(ig,ilay)*z1(ig) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c rajout eventuel de planck dans le shema implicite: -c -------------------------------------------------- - - z4st=4.*5.67e-8*ptimestep -c z4st=0. - DO ig=1,ngrid - zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) - ENDDO - -c----------------------------------------------------------------------- -c calcul le l'evolution de la temperature du sol: -c ----------------------------------------------- - - DO ig=1,ngrid - z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) - s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep - z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) - ztsrf2(ig)=z1(ig)/z2(ig) - zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) - pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep - ENDDO - -c----------------------------------------------------------------------- -c integration verticale finale: -c ----------------------------- - - DO ilay=2,nlay - DO ig=1,ngrid - zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c calcul final des tendances de la diffusion verticale: -c ----------------------------------------------------- - - DO ilev = 1, nlay - DO ig=1,ngrid - pdudif(ig,ilev)=( zu(ig,ilev)- - $ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep - pdvdif(ig,ilev)=( zv(ig,ilev)- - $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep - pdhdif(ig,ilev)=( zh(ig,ilev)- - $ (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep - ENDDO - ENDDO - - IF(lwrite) THEN - PRINT* - PRINT*,'Diagnostique de la diffusion verticale' - PRINT*,'h avant et apres diffusion verticale' - PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1) - DO 3110 ilev=1,nlay - PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev) -3110 CONTINUE - ENDIF - - -c----------------------------------------------------------------------- - - RETURN - END