Index: LMDZ4/branches/V3_test/libf/dyn3dpar/PVtheta.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/PVtheta.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/PVtheta.F (revision 708) @@ -0,0 +1,196 @@ + SUBROUTINE PVtheta(ilon,ilev,pucov,pvcov,pteta, + $ ztfi,zplay,zplev, + $ nbteta,theta,PVteta) + IMPLICIT none + +c======================================================================= +c +c Auteur: I. Musat +c ------- +c +c Objet: +c ------ +c +c ******************************************************************* +c Calcul de la vorticite potentielle PVteta sur des iso-theta selon +c la methodologie du NCEP/NCAR : +c 1) on calcule la stabilite statique N**2=g/T*(dT/dz+g/cp) sur les +c niveaux du modele => N2 +c 2) on interpole les vents, la temperature et le N**2 sur des isentropes +c (en fait sur des iso-theta) lineairement en log(theta) => +c ucovteta, vcovteta, N2teta +c 3) on calcule la vorticite absolue sur des iso-theta => vorateta +c 4) on calcule la densite rho sur des iso-theta => rhoteta +c +c rhoteta = (T/theta)**(cp/R)*p0/(R*T) +c +c 5) on calcule la vorticite potentielle sur des iso-theta => PVteta +c +c PVteta = (vorateta * N2 * theta)/(g * rhoteta) ! en PVU +c +c NB: 1PVU=10**(-6) K*m**2/(s * kg) +c +c PVteta = vorateta * N2/(g**2 * rhoteta) ! en 1/(Pa*s) +c +c +c ******************************************************************* +c +c +c Variables d'entree : ilon,ilev,pucov,pvcov,pteta,ztfi,zplay,zplev,nbteta,theta +c -> sur la grille dynamique +c Variable de sortie : PVteta +c -> sur la grille physique +c======================================================================= + +#include "dimensions.h" +#include "paramet.h" +c +c variables Input +c + INTEGER ilon, ilev + REAL pvcov(iip1,jjm,ilev) + REAL pucov(iip1,jjp1,ilev) + REAL pteta(iip1,jjp1,ilev) + REAL ztfi(ilon,ilev) + REAL zplay(ilon,ilev), zplev(ilon,ilev+1) + INTEGER nbteta + REAL theta(nbteta) +c +c variable Output +c + REAL PVteta(ilon,nbteta) +c +c variables locales +c + INTEGER i, j, l, ig0 + REAL SSUM + REAL teta(ilon, ilev) + REAL ptetau(ip1jmp1, ilev), ptetav(ip1jm, ilev) + REAL ucovteta(ip1jmp1,ilev), vcovteta(ip1jm,ilev) + REAL N2(ilon,ilev-1), N2teta(ilon,nbteta) + REAL ztfiteta(ilon,nbteta) + REAL rhoteta(ilon,nbteta) + REAL vorateta(iip1,jjm,nbteta) + REAL voratetafi(ilon,nbteta), vorpol(iim) +c +#include "comgeom2.h" +#include "comconst.h" +#include "comvert.h" +c +c projection teta sur la grille physique +c + DO l=1,llm + teta(1,l) = pteta(1,1,l) + ig0 = 2 + DO j = 2, jjm + DO i = 1, iim + teta(ig0,l) = pteta(i,j,l) + ig0 = ig0 + 1 + ENDDO + ENDDO + teta(ig0,l) = pteta(1,jjp1,l) + ENDDO +c +c calcul pteta sur les grilles U et V +c + DO l=1, llm + DO j=1, jjp1 + DO i=1, iip1 + ig0=i+(j-1)*iip1 + ptetau(ig0,l)=pteta(i,j,l) + ENDDO !i + ENDDO !j + DO j=1, jjm + DO i=1, iip1 + ig0=i+(j-1)*iip1 + ptetav(ig0,l)=0.5*(pteta(i,j,l)+pteta(i,j+1,l)) + ENDDO !i + ENDDO !j + ENDDO !l +c +c projection pucov, pvcov sur une surface de theta constante +c + DO l=1, nbteta +cIM 1rout CALL tetaleveli1j1(ip1jmp1,llm,.true.,ptetau,theta(l), + CALL tetalevel(ip1jmp1,llm,.true.,ptetau,theta(l), + . pucov,ucovteta(:,l)) +cIM 1rout CALL tetaleveli1j(ip1jm,llm,.true.,ptetav,theta(l), + CALL tetalevel(ip1jm,llm,.true.,ptetav,theta(l), + . pvcov,vcovteta(:,l)) + ENDDO !l +c +c calcul vorticite absolue sur une iso-theta : vorateta +c + CALL tourabs(nbteta,vcovteta,ucovteta,vorateta) +c +c projection vorateta sur la grille physique => voratetafi +c + DO l=1,nbteta + DO j=2,jjm + ig0=1+(j-2)*iim + DO i=1,iim + voratetafi(ig0+i+1,l) = vorateta( i ,j-1,l) * alpha4(i+1,j) + + $ vorateta(i+1,j-1,l) * alpha1(i+1,j) + + $ vorateta(i ,j ,l) * alpha3(i+1,j) + + $ vorateta(i+1,j ,l) * alpha2(i+1,j) + ENDDO + voratetafi(ig0 +1,l) = voratetafi(ig0 +1+ iim,l) + ENDDO + ENDDO +c + DO l=1,nbteta + DO i=1,iim + vorpol(i) = vorateta(i,1,l)*aire(i,1) + ENDDO + voratetafi(1,l)= SSUM(iim,vorpol,1)/apoln + ENDDO +c + DO l=1,nbteta + DO i=1,iim + vorpol(i) = vorateta(i,jjm,l)* aire(i,jjm +1) + ENDDO + voratetafi(ilon,l)= SSUM(iim,vorpol,1)/apols + ENDDO +c +c calcul N**2 sur la grille physique => N2 +c + DO l=1, llm-1 + DO i=1, ilon + N2(i,l) = (g**2 * zplay(i,l) * + $ (ztfi(i,l+1)-ztfi(i,l)) )/ + $ (R*ztfi(i,l)*ztfi(i,l)* + $ (zplev(i,l)-zplev(i,l+1)) )+ + $ (g**2)/(ztfi(i,l)*CPP) + ENDDO !i + ENDDO !l +c +c calcul N2 sur une iso-theta => N2teta +c + DO l=1, nbteta + CALL tetalevel(ilon,llm-1,.true.,teta,theta(l), + $ N2,N2teta(:,l)) + CALL tetalevel(ilon,llm,.true.,teta,theta(l), + $ ztfi,ztfiteta(:,l)) + ENDDO !l=1, nbteta +c +c calcul rho et PV sur une iso-theta : rhoteta, PVteta +c + DO l=1, nbteta + DO i=1, ilon + rhoteta(i,l)=(ztfiteta(i,l)/theta(l))**(CPP/R)* + $ (preff/(R*ztfiteta(i,l))) +c +c PVteta en PVU +c + PVteta(i,l)=(theta(l)*g*voratetafi(i,l)*N2teta(i,l))/ + $ (g**2*rhoteta(i,l)) +c +c PVteta en 1/(Pa*s) +c + PVteta(i,l)=(voratetafi(i,l)*N2teta(i,l))/ + $ (g**2*rhoteta(i,l)) + ENDDO !i + ENDDO !l +c + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/advect_new_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/advect_new_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/advect_new_p.F (revision 708) @@ -0,0 +1,281 @@ + SUBROUTINE advect_new_p(ucov,vcov,teta,w,massebx,masseby, + & du,dv,dteta) + USE parallel + USE write_field_p + IMPLICIT NONE +c======================================================================= +c +c Auteurs: P. Le Van , Fr. Hourdin . +c ------- +c +c Objet: +c ------ +c +c ************************************************************* +c .... calcul des termes d'advection vertic.pour u,v,teta,q ... +c ************************************************************* +c ces termes sont ajoutes a du,dv,dteta et dq . +c Modif F.Forget 03/94 : on retire q de advect +c +c======================================================================= +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comvert.h" +#include "comgeom.h" +#include "logic.h" +#include "ener.h" + +c Arguments: +c ---------- + + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm) + REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),w(ip1jmp1,llm) + REAL dv(ip1jm,llm),du(ip1jmp1,llm),dteta(ip1jmp1,llm) + REAL,SAVE :: dv1(ip1jm,llm),du1(ip1jmp1,llm),dteta1(ip1jmp1,llm) + REAL,SAVE :: dv2(ip1jm,llm),du2(ip1jmp1,llm),dteta2(ip1jmp1,llm) +c Local: +c ------ + + REAL,SAVE :: uav(ip1jmp1,llm),vav(ip1jm,llm) + REAL wsur2(ip1jmp1) + REAL unsaire2(ip1jmp1), ge(ip1jmp1) + REAL deuxjour, ww, gt, uu, vv + + INTEGER ij,l,ijb,ije + + EXTERNAL SSUM + REAL SSUM + +c----------------------------------------------------------------------- +c 2. Calculs preliminaires: +c ------------------------- + + IF (conser) THEN + deuxjour = 2. * daysec + + DO 1 ij = 1, ip1jmp1 + unsaire2(ij) = unsaire(ij) * unsaire(ij) + 1 CONTINUE + END IF + + +c------------------ -yy ---------------------------------------------- +c . Calcul de u + +c$OMP MASTER + ijb=ij_begin + ije=ij_end + if (pole_nord) ijb=ijb+iip1 + if (pole_sud) ije=ije-iip1 + + DO ij=ijb,ije + du2(ij,1)=0. + ENDDO + + ijb=ij_begin + ije=ij_end + if (pole_sud) ije=ij_end-iip1 + + DO ij=ijb,ije + dv2(ij,1)=0. + ENDDO + + ijb=ij_begin + ije=ij_end + + DO ij=ijb,ije + dteta2(ij,1)=0. + ENDDO +c$OMP END MASTER + + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + + ijb=ij_begin + ije=ij_end + if (pole_nord) ijb=ijb+iip1 + if (pole_sud) ije=ije-iip1 + +c DO ij = iip2, ip1jmp1 +c uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) ) +c ENDDO + +c DO ij = iip2, ip1jm +c uav(ij,l) = uav(ij,l) + uav(ij+iip1,l) +c ENDDO + + DO ij = ijb, ije + + uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l)) + . +0.25*(ucov(ij+iip1,l)+ucov(ij,l)) + ENDDO + + if (pole_nord) then + DO ij = 1, iip1 + uav(ij ,l) = 0. + ENDDO + endif + + if (pole_sud) then + DO ij = 1, iip1 + uav(ip1jm+ij,l) = 0. + ENDDO + endif + + ENDDO +c$OMP END DO +c call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/))) + +c------------------ -xx ---------------------------------------------- +c . Calcul de v + + ijb=ij_begin + ije=ij_end + if (pole_sud) ije=ij_end-iip1 + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + + DO ij = ijb+1, ije + vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) ) + ENDDO + + DO ij = ijb,ije,iip1 + vav(ij,l) = vav(ij+iim,l) + ENDDO + + + DO ij = ijb, ije-1 + vav(ij,l) = vav(ij,l) + vav(ij+1,l) + ENDDO + + DO ij = ijb, ije, iip1 + vav(ij+iim,l) = vav(ij,l) + ENDDO + + ENDDO +c$OMP END DO +c call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/))) + +c----------------------------------------------------------------------- +c$OMP BARRIER + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO 20 l = 1, llmm1 + + +c ...... calcul de - w/2. au niveau l+1 ....... + ijb=ij_begin + ije=ij_end+iip1 + if (pole_sud) ije=ij_end + + DO 5 ij = ijb, ije + wsur2( ij ) = - 0.5 * w( ij,l+1 ) + 5 CONTINUE + + +c ..................... calcul pour du .................. + + ijb=ij_begin + ije=ij_end + if (pole_nord) ijb=ijb+iip1 + if (pole_sud) ije=ije-iip1 + + DO 6 ij = ijb ,ije-1 + ww = wsur2 ( ij ) + wsur2( ij+1 ) + uu = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) ) + du1(ij,l) = ww * ( uu - uav(ij, l ) )/massebx(ij, l ) + du2(ij,l+1)= ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1) + 6 CONTINUE + +c ................. calcul pour dv ..................... + ijb=ij_begin + ije=ij_end + if (pole_sud) ije=ij_end-iip1 + + DO 8 ij = ijb, ije + ww = wsur2( ij+iip1 ) + wsur2( ij ) + vv = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) ) + dv1(ij,l) = ww * (vv - vav(ij, l ) )/masseby(ij, l ) + dv2(ij,l+1)= ww * (vv - vav(ij,l+1) )/masseby(ij,l+1) + 8 CONTINUE + +c + +c ............................................................ +c ............... calcul pour dh ................... +c ............................................................ + +c ---z +c calcul de - d( teta * w ) qu'on ajoute a dh +c ............... + ijb=ij_begin + ije=ij_end + + DO 15 ij = ijb, ije + ww = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) ) + dteta1(ij, l ) = ww + dteta2(ij,l+1) = ww + 15 CONTINUE + +c ym ---> conser a voir plus tard + +c IF( conser) THEN +c +c DO 17 ij = 1,ip1jmp1 +c ge(ij) = wsur2(ij) * wsur2(ij) * unsaire2(ij) +c 17 CONTINUE +c gt = SSUM( ip1jmp1,ge,1 ) +c gtot(l) = deuxjour * SQRT( gt/ip1jmp1 ) +c END IF + + 20 CONTINUE +c$OMP END DO + + ijb=ij_begin + ije=ij_end + if (pole_nord) ijb=ijb+iip1 + if (pole_sud) ije=ije-iip1 + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije-1 + du(ij,l)=du(ij,l)+du2(ij,l)-du1(ij,l) + ENDDO + + DO ij = ijb+iip1-1, ije, iip1 + du( ij, l ) = du( ij -iim, l ) + ENDDO + ENDDO +c$OMP END DO NOWAIT + + ijb=ij_begin + ije=ij_end + if (pole_sud) ije=ij_end-iip1 + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije + dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l) + ENDDO + ENDDO +c$OMP END DO NOWAIT + ijb=ij_begin + ije=ij_end + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije + dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l) + ENDDO + ENDDO +c$OMP END DO NOWAIT + + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/control.inc =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/control.inc (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/control.inc (revision 708) @@ -0,0 +1,18 @@ +! +! $Header$ +! +!----------------------------------------------------------------------- +! INCLUDE 'control.h' + + COMMON/control/nday,day_step, & + & iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq , & + & periodav,ecritphy,iecrimoy,dayref,anneeref, & + & raz_date,offline,ip_ebil_dyn + + INTEGER nday,day_step,iperiod,iapp_tracvl,iconser,iecri, & + & idissip,iphysiq,iecrimoy,dayref,anneeref, raz_date, & + & ip_ebil_dyn + REAL periodav, ecritphy + logical offline + +!----------------------------------------------------------------------- Index: LMDZ4/branches/V3_test/libf/dyn3dpar/convmas1_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/convmas1_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/convmas1_p.F (revision 708) @@ -0,0 +1,62 @@ + SUBROUTINE convmas1_p (pbaru, pbarv, convm ) +c + USE parallel + IMPLICIT NONE + +c======================================================================= +c +c Auteurs: P. Le Van , F. Hourdin . +c ------- +c +c Objet: +c ------ +c +c ******************************************************************** +c .... calcul de la convergence du flux de masse aux niveaux p ... +c ******************************************************************** +c +c +c pbaru et pbarv sont des arguments d'entree pour le s-pg .... +c ..... convm est un argument de sortie pour le s-pg .... +c +c le calcul se fait de haut en bas, +c la convergence de masse au niveau p(llm+1) est egale a 0. et +c n'est pas stockee dans le tableau convm . +c +c +c======================================================================= +c +c Declarations: +c ------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "logic.h" + + REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) + REAL, target :: convm( ip1jmp1,llm ) + INTEGER l,ij + + INTEGER ijb,ije,jjb,jje + + +c----------------------------------------------------------------------- +c .... calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ...... + + CALL convflu_p( pbaru, pbarv, llm, convm ) + +c----------------------------------------------------------------------- +c filtrage: +c --------- + + jjb=jj_begin + jje=jj_end+1 + if (pole_sud) jje=jj_end + + CALL filtreg_p( convm, jjb, jje, jjp1, llm, 2, 2, .true., 1 ) + +c integration de la convergence de masse de haut en bas ...... +c + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/convmas2_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/convmas2_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/convmas2_p.F (revision 708) @@ -0,0 +1,56 @@ + SUBROUTINE convmas2_p ( convm ) +c + USE parallel + IMPLICIT NONE + +c======================================================================= +c +c Auteurs: P. Le Van , F. Hourdin . +c ------- +c +c Objet: +c ------ +c +c ******************************************************************** +c .... calcul de la convergence du flux de masse aux niveaux p ... +c ******************************************************************** +c +c +c pbaru et pbarv sont des arguments d'entree pour le s-pg .... +c ..... convm est un argument de sortie pour le s-pg .... +c +c le calcul se fait de haut en bas, +c la convergence de masse au niveau p(llm+1) est egale a 0. et +c n'est pas stockee dans le tableau convm . +c +c +c======================================================================= +c +c Declarations: +c ------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "logic.h" + + REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) + REAL :: convm( ip1jmp1,llm ) + INTEGER l,ij + INTEGER ijb,ije,jjb,jje + +c$OMP MASTER +c integration de la convergence de masse de haut en bas ...... + ijb=ij_begin + ije=ij_end+iip1 + if (pole_sud) ije=ij_end + + DO l = llmm1, 1, -1 + DO ij = ijb, ije + convm(ij,l) = convm(ij,l) + convm(ij,l+1) + ENDDO + ENDDO +c +c$OMP END MASTER + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/filtreg_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/filtreg_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/filtreg_p.F (revision 708) @@ -0,0 +1,303 @@ + SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv, + . ifiltre, iaire, griscal ,iter) + USE Parallel, only : OMP_CHUNK + IMPLICIT NONE + +c======================================================================= +c +c Auteur: P. Le Van 07/10/97 +c ------ +c +c Objet: filtre matriciel longitudinal ,avec les matrices precalculees +c pour l'operateur Filtre . +c ------ +c +c Arguments: +c ---------- +c +c +c ibeg..iend lattitude a filtrer +c nlat nombre de latitudes du champ +c nbniv nombre de niveaux verticaux a filtrer +c champ(iip1,nblat,nbniv) en entree : champ a filtrer +c en sortie : champ filtre +c ifiltre +1 Transformee directe +c -1 Transformee inverse +c +2 Filtre directe +c -2 Filtre inverse +c +c iaire 1 si champ intensif +c 2 si champ extensif (pondere par les aires) +c +c iter 1 filtre simple +c +c======================================================================= +c +c +c Variable Intensive +c ifiltre = 1 filtre directe +c ifiltre =-1 filtre inverse +c +c Variable Extensive +c ifiltre = 2 filtre directe +c ifiltre =-2 filtre inverse +c +c +#include "dimensions.h" +#include "paramet.h" +#include "parafilt.h" +#include "coefils.h" +c + INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter + INTEGER i,j,l,k + INTEGER iim2,immjm + INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil + + REAL champ( iip1,nlat,nbniv) + REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs + COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus) + , , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs) + , , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus) + REAL eignq(iim), sdd1(iim),sdd2(iim) + LOGICAL griscal + INTEGER hemisph, iaire +c + + IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) + * STOP'Pas de transformee simple dans cette version' + + IF( iter.EQ. 2 ) THEN + PRINT *,' Pas d iteration du filtre dans cette version !' + * , ' Utiliser old_filtreg et repasser !' + STOP + ENDIF + + IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN + PRINT *,' Cette routine ne calcule le filtre inverse que ', + * ' sur la grille des scalaires !' + STOP + ENDIF + + IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN + PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' + *,' corriger et repasser !' + STOP + ENDIF +c + + iim2 = iim * iim + immjm = iim * jjm +c +c + IF( griscal ) THEN + IF( nlat. NE. jjp1 ) THEN + PRINT 1111 + STOP + ELSE +c + IF( iaire.EQ.1 ) THEN + CALL SCOPY( iim, sddv, 1, sdd1, 1 ) + CALL SCOPY( iim, unsddv, 1, sdd2, 1 ) + ELSE + CALL SCOPY( iim, unsddv, 1, sdd1, 1 ) + CALL SCOPY( iim, sddv, 1, sdd2, 1 ) + END IF +c + jdfil1 = 2 + jffil1 = jfiltnu + jdfil2 = jfiltsu + jffil2 = jjm + END IF + ELSE + IF( nlat.NE.jjm ) THEN + PRINT 2222 + STOP + ELSE +c + IF( iaire.EQ.1 ) THEN + CALL SCOPY( iim, sddu, 1, sdd1, 1 ) + CALL SCOPY( iim, unsddu, 1, sdd2, 1 ) + ELSE + CALL SCOPY( iim, unsddu, 1, sdd1, 1 ) + CALL SCOPY( iim, sddu, 1, sdd2, 1 ) + END IF +c + jdfil1 = 1 + jffil1 = jfiltnv + jdfil2 = jfiltsv + jffil2 = jjm + END IF + END IF +c +c + DO 100 hemisph = 1, 2 +c + IF ( hemisph.EQ.1 ) THEN +c ym + jdfil = max(jdfil1,ibeg) + jffil = min(jffil1,iend) + ELSE +c ym + jdfil = max(jdfil2,ibeg) + jffil = min(jffil2,iend) + END IF + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO 50 l = 1, nbniv + DO 30 j = jdfil,jffil + + + DO 5 i = 1, iim + champ(i,j,l) = champ(i,j,l) * sdd1(i) + 5 CONTINUE +c + + IF( hemisph. EQ. 1 ) THEN + + IF( ifiltre. EQ. -2 ) THEN +#ifdef CRAY + CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , + * 1, iim, iim ) +#else +#ifdef BLAS + CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim, + . champ(1,j,l), 1, 0.0, eignq, 1) +#else + DO k = 1, iim + eignq(k) = 0.0 + ENDDO + DO k = 1, iim + DO i = 1, iim + eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l) + ENDDO + ENDDO +#endif +#endif + ELSE IF ( griscal ) THEN +#ifdef CRAY + CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq , + * 1, iim, iim ) +#else +#ifdef BLAS + CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim, + . champ(1,j,l), 1, 0.0, eignq, 1) +#else + DO k = 1, iim + eignq(k) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim + eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l) + ENDDO + ENDDO +#endif +#endif + ELSE +#ifdef CRAY + CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , + * 1, iim, iim ) +#else +#ifdef BLAS + CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim, + . champ(1,j,l), 1, 0.0, eignq, 1) +#else + DO k = 1, iim + eignq(k) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim + eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l) + ENDDO + ENDDO +#endif +#endif + ENDIF + + ELSE + + IF( ifiltre. EQ. -2 ) THEN +#ifdef CRAY + CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , + * eignq, 1, iim, iim ) +#else +#ifdef BLAS + CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim, + . champ(1,j,l), 1, 0.0, eignq, 1) +#else + DO k = 1, iim + eignq(k) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim + eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) + ENDDO + ENDDO +#endif +#endif + ELSE IF ( griscal ) THEN +#ifdef CRAY + CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , + * eignq, 1, iim, iim ) +#else +#ifdef BLAS + CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim, + . champ(1,j,l), 1, 0.0, eignq, 1) +#else + DO k = 1, iim + eignq(k) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim + eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) + ENDDO + ENDDO +#endif +#endif + ELSE +#ifdef CRAY + CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 , + * eignq, 1, iim, iim ) +#else +#ifdef BLAS + CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim, + . champ(1,j,l), 1, 0.0, eignq, 1) +#else + DO k = 1, iim + eignq(k) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim + eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) + ENDDO + ENDDO +#endif +#endif + ENDIF + + ENDIF +c + IF( ifiltre.EQ. 2 ) THEN + DO 15 i = 1, iim + champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) + 15 CONTINUE + ELSE + DO 16 i=1,iim + champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) +16 CONTINUE + ENDIF +c + champ( iip1,j,l ) = champ( 1,j,l ) +c + 30 CONTINUE +c + 50 CONTINUE +c$OMP END DO NOWAIT +c + 100 CONTINUE +c +1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a + *filtrer, sur la grille des scalaires'/) +2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi + *ltrer, sur la grille de V ou de Z'/) + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/initdynav_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/initdynav_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/initdynav_p.F (revision 708) @@ -0,0 +1,176 @@ +! +! $Header$ +! +c +c + subroutine initdynav_p(infile,day0,anne0,tstep,t_ops,t_wrt + . ,nq,fileid) + + USE IOIPSL + use parallel + use Write_field + use misc_mod + implicit none + +C +C Routine d'initialisation des ecritures des fichiers histoires LMDZ +C au format IOIPSL. Initialisation du fichier histoire moyenne. +C +C Appels succesifs des routines: histbeg +C histhori +C histver +C histdef +C histend +C +C Entree: +C +C infile: nom du fichier histoire a creer +C day0,anne0: date de reference +C tstep : frequence d'ecriture +C t_ops: frequence de l'operation pour IOIPSL +C t_wrt: frequence d'ecriture sur le fichier +C nq: nombre de traceurs +C +C Sortie: +C fileid: ID du fichier netcdf cree +C +C L. Fairhead, LMD, 03/99 +C +C ===================================================================== +C +C Declarations +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comvert.h" +#include "comgeom.h" +#include "temps.h" +#include "ener.h" +#include "logic.h" +#include "description.h" +#include "serre.h" +#include "advtrac.h" + +C Arguments +C + character*(*) infile + integer*4 day0, anne0 + real tstep, t_ops, t_wrt + integer fileid + integer nq + integer thoriid, zvertiid + +C Variables locales +C + integer tau0 + real zjulian + integer iq + real rlong(iip1,jjp1), rlat(iip1,jjp1) + integer ii,jj + integer zan, dayref + integer :: jjb,jje,jjn + + if (adjust) return +C +C Initialisations +C + pi = 4. * atan (1.) +C +C Appel a histbeg: creation du fichier netcdf et initialisations diverses +C + + zan = anne0 + dayref = day0 + CALL ymds2ju(zan, 1, dayref, 0.0, zjulian) + tau0 = itau_dyn + + do jj = 1, jjp1 + do ii = 1, iip1 + rlong(ii,jj) = rlonv(ii) * 180. / pi + rlat(ii,jj) = rlatu(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + + call histbeg(trim(infile)//'_'//trim(int2str(mpi_rank))//'.nc', + . iip1, rlong(:,1), jjn, rlat(1,jjb:jje), + . 1, iip1, 1, jjn, + . tau0, zjulian, tstep, thoriid, fileid) + +C +C Appel a histvert pour la grille verticale +C + call histvert(fileid, 'sigss', 'Niveaux sigma','Pa', + . llm, nivsigs, zvertiid) +C +C Appels a histdef pour la definition des variables a sauvegarder +C +C Vents U +C + write(6,*)'inithistave',tstep + call histdef(fileid, 'u', 'vents u scalaires moyennes', + . 'm/s', iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'ave(X)', t_ops, t_wrt) + +C +C Vents V +C + call histdef(fileid, 'v', 'vents v scalaires moyennes', + . 'm/s', iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'ave(X)', t_ops, t_wrt) + +C +C Temperature +C + call histdef(fileid, 'temp', 'temperature moyennee', 'K', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'ave(X)', t_ops, t_wrt) +C +C Temperature potentielle +C + call histdef(fileid, 'theta', 'temperature potentielle', 'K', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'ave(X)', t_ops, t_wrt) + + +C +C Geopotentiel +C + call histdef(fileid, 'phi', 'geopotentiel moyenne', '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'ave(X)', t_ops, t_wrt) +C +C Traceurs +C + DO iq=1,nq + call histdef(fileid, ttext(iq), ttext(iq), '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'ave(X)', t_ops, t_wrt) + enddo +C +C Masse +C + call histdef(fileid, 'masse', 'masse', 'kg', + . iip1, jjn, thoriid, 1, 1, 1, -99, + . 32, 'ave(X)', t_ops, t_wrt) +C +C Pression au sol +C + call histdef(fileid, 'ps', 'pression naturelle au sol', 'Pa', + . iip1, jjn, thoriid, 1, 1, 1, -99, + . 32, 'ave(X)', t_ops, t_wrt) +C +C Pression au sol +C + call histdef(fileid, 'phis', 'geopotentiel au sol', '-', + . iip1, jjn, thoriid, 1, 1, 1, -99, + . 32, 'ave(X)', t_ops, t_wrt) +C +C Fin +C + call histend(fileid) + return + end Index: LMDZ4/branches/V3_test/libf/dyn3dpar/initfluxsto_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/initfluxsto_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/initfluxsto_p.F (revision 708) @@ -0,0 +1,255 @@ +! +! $Header$ +! + subroutine initfluxsto_p + . (infile,tstep,t_ops,t_wrt,nq, + . fileid,filevid,filedid) + + USE IOIPSL + use parallel + use Write_field + use misc_mod + + implicit none + +C +C Routine d'initialisation des ecritures des fichiers histoires LMDZ +C au format IOIPSL +C +C Appels succesifs des routines: histbeg +C histhori +C histver +C histdef +C histend +C +C Entree: +C +C infile: nom du fichier histoire a creer +C day0,anne0: date de reference +C tstep: duree du pas de temps en seconde +C t_ops: frequence de l'operation pour IOIPSL +C t_wrt: frequence d'ecriture sur le fichier +C nq: nombre de traceurs +C +C Sortie: +C fileid: ID du fichier netcdf cree +C filevid:ID du fichier netcdf pour la grille v +C +C L. Fairhead, LMD, 03/99 +C +C ===================================================================== +C +C Declarations +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comvert.h" +#include "comgeom.h" +#include "temps.h" +#include "ener.h" +#include "logic.h" +#include "description.h" +#include "serre.h" + +C Arguments +C + character*(*) infile + integer*4 itau + real tstep, t_ops, t_wrt + integer fileid, filevid,filedid + integer nq,ndex(1) + real nivd(1) + +C Variables locales +C + integer tau0 + real zjulian + character*3 str + character*10 ctrac + integer iq + real rlong(iip1,jjp1), rlat(iip1,jjp1),rl(1,1) + integer uhoriid, vhoriid, thoriid, zvertiid,dhoriid,dvertiid + integer ii,jj + integer zan, idayref + logical ok_sync + integer :: jjb,jje,jjn +C +C Initialisations +C + pi = 4. * atan (1.) + str='q ' + ctrac = 'traceur ' + ok_sync = .true. +C +C Appel a histbeg: creation du fichier netcdf et initialisations diverses +C + + zan = annee_ref + idayref = day_ref + CALL ymds2ju(zan, 1, idayref, 0.0, zjulian) + tau0 = itau_dyn + + do jj = 1, jjp1 + do ii = 1, iip1 + rlong(ii,jj) = rlonu(ii) * 180. / pi + rlat(ii,jj) = rlatu(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + + call histbeg(trim(infile)//'_'//trim(int2str(mpi_rank))//'.nc', + . iip1, rlong(:,1), jjp1, rlat(1,jjb:jje), + . 1, iip1, 1, jjn, + . tau0, zjulian, tstep, uhoriid, fileid) +C +C Creation du fichier histoire pour la grille en V (oblige pour l'instant, +C IOIPSL ne permet pas de grilles avec des nombres de point differents dans +C un meme fichier) + + + do jj = 1, jjm + do ii = 1, iip1 + rlong(ii,jj) = rlonv(ii) * 180. / pi + rlat(ii,jj) = rlatv(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + if (pole_sud) jje=jj_end-1 + if (pole_sud) jjn=jj_nb-1 + + call histbeg('fluxstokev_'//trim(int2str(mpi_rank))//'.nc', + . iip1, rlong(:,1), jjm, rlat(1,jjb:jje), + . 1, iip1, 1, jjn, + . tau0, zjulian, tstep, vhoriid, filevid) + + rl(1,1) = 1. + + if (mpi_rank==0) then + + call histbeg('defstoke.nc', 1, rl, 1, rl, + . 1, 1, 1, 1, + . tau0, zjulian, tstep, dhoriid, filedid) + + endif +C +C Appel a histhori pour rajouter les autres grilles horizontales +C + do jj = 1, jjp1 + do ii = 1, iip1 + rlong(ii,jj) = rlonv(ii) * 180. / pi + rlat(ii,jj) = rlatu(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + + call histhori(fileid, iip1, rlong(:,jjb:jje),jjn,rlat(:,jjb:jje), + . 'scalar','Grille points scalaires', thoriid) + +C +C Appel a histvert pour la grille verticale +C + call histvert(fileid, 'sig_s', 'Niveaux sigma', + . 'sigma_level', + . llm, nivsigs, zvertiid) +C Pour le fichier V + call histvert(filevid, 'sig_s', 'Niveaux sigma', + . 'sigma_level', + . llm, nivsigs, zvertiid) +c pour le fichier def + nivd(1) = 1 + call histvert(filedid, 'sig_s', 'Niveaux sigma', + . 'sigma_level', + . 1, nivd, dvertiid) + +C +C Appels a histdef pour la definition des variables a sauvegarder + + CALL histdef(fileid, "phis", "Surface geop. height", "-", + . iip1,jjn,thoriid, 1,1,1, -99, 32, + . "once", t_ops, t_wrt) + + CALL histdef(fileid, "aire", "Grid area", "-", + . iip1,jjn,thoriid, 1,1,1, -99, 32, + . "once", t_ops, t_wrt) + + if (mpi_rank==0) then + + CALL histdef(filedid, "dtvr", "tps dyn", "s", + . 1,1,dhoriid, 1,1,1, -99, 32, + . "once", t_ops, t_wrt) + + CALL histdef(filedid, "istdyn", "tps stock", "s", + . 1,1,dhoriid, 1,1,1, -99, 32, + . "once", t_ops, t_wrt) + + CALL histdef(filedid, "istphy", "tps stock phy", "s", + . 1,1,dhoriid, 1,1,1, -99, 32, + . "once", t_ops, t_wrt) + + endif +C +C Masse +C + call histdef(fileid, 'masse', 'Masse', 'kg', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Pbaru +C + call histdef(fileid, 'pbaru', 'flx de masse zonal', 'kg m/s', + . iip1, jjn, uhoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) + +C +C Pbarv +C + if (pole_sud) jjn=jj_nb-1 + + call histdef(filevid, 'pbarv', 'flx de masse mer', 'kg m/s', + . iip1, jjn, vhoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C +C w +C + if (pole_sud) jjn=jj_nb + call histdef(fileid, 'w', 'flx de masse vert', 'kg m/s', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) + +C +C Temperature potentielle +C + call histdef(fileid, 'teta', 'temperature potentielle', '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C + +C +C Geopotentiel +C + call histdef(fileid, 'phi', 'geopotentiel instantane', '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Fin +C + call histend(fileid) + call histend(filevid) + call histend(filedid) + if (ok_sync) then + call histsync(fileid) + call histsync(filevid) + call histsync(filedid) + endif + + return + end Index: LMDZ4/branches/V3_test/libf/dyn3dpar/inithist_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/inithist_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/inithist_p.F (revision 708) @@ -0,0 +1,214 @@ +! +! $Header$ +! + subroutine inithist_p(infile,day0,anne0,tstep,t_ops,t_wrt,nq, + . fileid,filevid) + + USE IOIPSL + use parallel + use Write_field + use misc_mod + + implicit none + +C +C Routine d'initialisation des ecritures des fichiers histoires LMDZ +C au format IOIPSL +C +C Appels succesifs des routines: histbeg +C histhori +C histver +C histdef +C histend +C +C Entree: +C +C infile: nom du fichier histoire a creer +C day0,anne0: date de reference +C tstep: duree du pas de temps en seconde +C t_ops: frequence de l'operation pour IOIPSL +C t_wrt: frequence d'ecriture sur le fichier +C nq: nombre de traceurs +C +C Sortie: +C fileid: ID du fichier netcdf cree +C filevid:ID du fichier netcdf pour la grille v +C +C L. Fairhead, LMD, 03/99 +C +C ===================================================================== +C +C Declarations +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comvert.h" +#include "comgeom.h" +#include "temps.h" +#include "ener.h" +#include "logic.h" +#include "description.h" +#include "serre.h" +#include "advtrac.h" + +C Arguments +C + character*(*) infile + integer*4 day0, anne0 + real tstep, t_ops, t_wrt + integer fileid, filevid + integer nq + +C Variables locales +C + integer tau0 + real zjulian + integer iq + real rlong(iip1,jjp1), rlat(iip1,jjp1) + integer uhoriid, vhoriid, thoriid, zvertiid + integer ii,jj + integer zan, dayref + integer :: jjb,jje,jjn +C +C Initialisations +C + if (adjust) return + + pi = 4. * atan (1.) +C +C Appel a histbeg: creation du fichier netcdf et initialisations diverses +C + + zan = anne0 + dayref = day0 + CALL ymds2ju(zan, 1, dayref, 0.0, zjulian) + tau0 = itau_dyn + + do jj = 1, jjp1 + do ii = 1, iip1 + rlong(ii,jj) = rlonu(ii) * 180. / pi + rlat(ii,jj) = rlatu(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + + call histbeg(trim(infile)//'_'//trim(int2str(mpi_rank))//'.nc', + . iip1, rlong(:,1), jjn, rlat(1,jjb:jje), + . 1, iip1, 1, jjn, + . tau0, zjulian, tstep, uhoriid, fileid) +C +C Creation du fichier histoire pour la grille en V (oblige pour l'instant, +C IOIPSL ne permet pas de grilles avec des nombres de point differents dans +C un meme fichier) + + do jj = 1, jjm + do ii = 1, iip1 + rlong(ii,jj) = rlonv(ii) * 180. / pi + rlat(ii,jj) = rlatv(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + if (pole_sud) jje=jj_end-1 + if (pole_sud) jjn=jj_nb-1 + + call histbeg('dyn_histv_'//trim(int2str(mpi_rank))//'.nc', + . iip1, rlong(:,1), jjn, rlat(1,jjb:jje), + . 1, iip1, 1, jjn, + . tau0, zjulian, tstep, vhoriid, filevid) +C +C Appel a histhori pour rajouter les autres grilles horizontales +C + + do jj = 1, jjp1 + do ii = 1, iip1 + rlong(ii,jj) = rlonv(ii) * 180. / pi + rlat(ii,jj) = rlatu(jj) * 180. / pi + enddo + enddo + + jjb=jj_begin + jje=jj_end + jjn=jj_nb + + call histhori(fileid, iip1, rlong(:,jjb:jje),jjn,rlat(:,jjb:jje), + . 'scalar','Grille points scalaires', thoriid) +C +C Appel a histvert pour la grille verticale +C + call histvert(fileid, 'sig_s', 'Niveaux sigma','-', + . llm, nivsigs, zvertiid) +C Pour le fichier V + call histvert(filevid, 'sig_s', 'Niveaux sigma','-', + . llm, nivsigs, zvertiid) +C +C Appels a histdef pour la definition des variables a sauvegarder +C +C Vents U +C + jjn=jj_nb + + call histdef(fileid, 'ucov', 'vents u covariants', 'm/s', + . iip1, jjn, uhoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Vents V +C + if (pole_sud) jjn=jj_nb-1 + + call histdef(filevid, 'vcov', 'vents v covariants', 'm/s', + . iip1, jjn, vhoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) + +C +C Temperature potentielle +C + jjn=jj_nb + + call histdef(fileid, 'teta', 'temperature potentielle', '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Geopotentiel +C + call histdef(fileid, 'phi', 'geopotentiel instantane', '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Traceurs +C + DO iq=1,nq + call histdef(fileid, ttext(iq), ttext(iq), '-', + . iip1, jjn, thoriid, llm, 1, llm, zvertiid, + . 32, 'inst(X)', t_ops, t_wrt) + enddo +C +C Masse +C + call histdef(fileid, 'masse', 'masse', 'kg', + . iip1, jjn, thoriid, 1, 1, 1, -99, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Pression au sol +C + call histdef(fileid, 'ps', 'pression naturelle au sol', 'Pa', + . iip1, jjn, thoriid, 1, 1, 1, -99, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Pression au sol +C + call histdef(fileid, 'phis', 'geopotentiel au sol', '-', + . iip1, jjn, thoriid, 1, 1, 1, -99, + . 32, 'inst(X)', t_ops, t_wrt) +C +C Fin +C + call histend(fileid) + call histend(filevid) + return + end Index: LMDZ4/branches/V3_test/libf/dyn3dpar/omp_chunk.h =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/omp_chunk.h (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/omp_chunk.h (revision 708) @@ -0,0 +1,1 @@ +#define OMP_CHUNK 5 Index: LMDZ4/branches/V3_test/libf/dyn3dpar/paramet90.h =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/paramet90.h (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/paramet90.h (revision 708) @@ -0,0 +1,17 @@ +!----------------------------------------------------------------------- +! INCLUDE 'paramet.h' (sauce fortran 90) + + INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1 + INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm + INTEGER ijmllm,mvar + INTEGER jcfil,jcfllm + + PARAMETER( iip1= iim+1-1/iim,iip2=iim+2,iip3=iim+3, & + jjp1=jjm+1-1/jjm) + PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 ) + PARAMETER( kftd = iim/2 -ndm ) + PARAMETER( ip1jm = iip1*jjm, ip1jmp1= iip1*jjp1 ) + PARAMETER( ip1jmi1= ip1jm - iip1 ) + PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm ) + PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm ) + PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm ) Index: LMDZ4/branches/V3_test/libf/dyn3dpar/tourabs.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/tourabs.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/tourabs.F (revision 708) @@ -0,0 +1,98 @@ + SUBROUTINE tourabs ( ntetaSTD,vcov, ucov, vorabs ) + IMPLICIT NONE + +c======================================================================= +c +c Modif: I. Musat (28/10/04) +c ------- +c adaptation du code tourpot.F pour le calcul de la vorticite absolue +c cf. P. Le Van +c +c Objet: +c ------ +c +c ******************************************************************* +c ............. calcul de la vorticite absolue ................. +c ******************************************************************* +c +c ntetaSTD, vcov,ucov sont des argum. d'entree pour le s-pg . +c vorabs est un argum.de sortie pour le s-pg . +c +c======================================================================= + +#include "dimensions.h" +#include "paramet.h" +#include "comgeom.h" +#include "logic.h" +#include "comconst.h" +c + INTEGER ntetaSTD + REAL vcov( ip1jm,ntetaSTD ), ucov( ip1jmp1,ntetaSTD ) + REAL vorabs( ip1jm,ntetaSTD ) +c +c variables locales + INTEGER l, ij, i, j + REAL rot( ip1jm,ntetaSTD ) + + + +c ... vorabs = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) .. + + + +c ........ Calcul du rotationnel du vent V puis filtrage ........ + + DO 5 l = 1,ntetaSTD + + DO 2 i = 1, iip1 + DO 2 j = 1, jjm +c + ij=i+(j-1)*iip1 + IF(ij.LE.ip1jm - 1) THEN +c + IF(cv(ij).EQ.0..OR.cv(ij+1).EQ.0..OR. + $ cu(ij).EQ.0..OR.cu(ij+iip1).EQ.0.) THEN + rot( ij,l ) = 0. + continue + ELSE + rot( ij,l ) = (vcov(ij+1,l)/cv(ij+1)-vcov(ij,l)/cv(ij))/ + $ (2.*pi*RAD*cos(rlatv(j)))*float(iim) + $ +(ucov(ij+iip1,l)/cu(ij+iip1)-ucov(ij,l)/cu(ij))/ + $ (pi*RAD)*(float(jjm)-1.) +c + ENDIF + ENDIF !(ij.LE.ip1jm - 1) THEN + 2 CONTINUE + +c .... correction pour rot( iip1,j,l ) ..... +c .... rot(iip1,j,l) = rot(1,j,l) ..... + +CDIR$ IVDEP + + DO 3 ij = iip1, ip1jm, iip1 + rot( ij,l ) = rot( ij -iim, l ) + 3 CONTINUE + + 5 CONTINUE + + + CALL filtreg( rot, jjm, ntetaSTD, 2, 1, .FALSE., 1 ) + + + DO 10 l = 1, ntetaSTD + + DO 6 ij = 1, ip1jm - 1 + vorabs( ij,l ) = ( rot(ij,l) + fext(ij)*unsairez(ij) ) + 6 CONTINUE + +c ..... correction pour vorabs( iip1,j,l) ..... +c .... vorabs(iip1,j,l)= vorabs(1,j,l) .... +CDIR$ IVDEP + DO 8 ij = iip1, ip1jm, iip1 + vorabs( ij,l ) = vorabs( ij -iim,l ) + 8 CONTINUE + + 10 CONTINUE + + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/vlspltgen_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/vlspltgen_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/vlspltgen_p.F (revision 708) @@ -0,0 +1,426 @@ +! +! $Header$ +! + SUBROUTINE vlspltgen_p( q,iadv,pente_max,masse,w,pbaru,pbarv,pdt, + , p,pk,teta ) +c +c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron +c +c ******************************************************************** +c Shema d'advection " pseudo amont " . +c + test sur humidite specifique: Q advecte< Qsat aval +c (F. Codron, 10/99) +c ******************************************************************** +c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... +c +c pente_max facteur de limitation des pentes: 2 en general +c 0 pour un schema amont +c pbaru,pbarv,w flux de masse en u ,v ,w +c pdt pas de temps +c +c teta temperature potentielle, p pression aux interfaces, +c pk exner au milieu des couches necessaire pour calculer Qsat +c -------------------------------------------------------------------- + USE parallel + USE mod_hallo + USE Write_Field_p + USE VAMPIR + IMPLICIT NONE + +c +#include "dimensions.h" +#include "paramet.h" +#include "logic.h" +#include "comvert.h" +#include "comconst.h" + +c +c Arguments: +c ---------- + INTEGER iadv(nqmx) + REAL masse(ip1jmp1,llm),pente_max + REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm) + REAL q(ip1jmp1,llm,nqmx) + REAL w(ip1jmp1,llm),pdt + REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm) +c +c Local +c --------- +c + INTEGER ij,l +c + REAL,SAVE :: qsat(ip1jmp1,llm) + REAL,SAVE :: zm(ip1jmp1,llm,nqmx) + REAL,SAVE :: mu(ip1jmp1,llm) + REAL,SAVE :: mv(ip1jm,llm) + REAL,SAVE :: mw(ip1jmp1,llm+1) + REAL,SAVE :: zq(ip1jmp1,llm,nqmx) + REAL zzpbar, zzw + + REAL qmin,qmax + DATA qmin,qmax/0.,1.e33/ + +c--pour rapport de melange saturant-- + + REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play + REAL ptarg,pdelarg,foeew,zdelta + REAL tempe(ip1jmp1) + INTEGER ijb,ije,iq + type(request) :: MyRequest1 + type(request) :: MyRequest2 + +c fonction psat(T) + + FOEEW ( PTARG,PDELARG ) = EXP ( + * (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) + * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) ) + + r2es = 380.11733 + r3les = 17.269 + r3ies = 21.875 + r4les = 35.86 + r4ies = 7.66 + retv = 0.6077667 + rtt = 273.16 + +c-- Calcul de Qsat en chaque point +c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2 +c pour eviter une exponentielle. +c$OMP MASTER + call SetTag(MyRequest1,100) + call SetTag(MyRequest2,101) +c$OMP END MASTER + + ijb=ij_begin-iip1 + ije=ij_end+iip1 + if (pole_nord) ijb=ij_begin + if (pole_sud) ije=ij_end + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l = 1, llm + DO ij = ijb, ije + tempe(ij) = teta(ij,l) * pk(ij,l) /cpp + ENDDO + DO ij = ijb, ije + zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) ) + play = 0.5*(p(ij,l)+p(ij,l+1)) + qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play ) + qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) ) + ENDDO + ENDDO +c$OMP END DO NOWAIT +c PRINT*,'Debut vlsplt version debug sans vlyqs' + + zzpbar = 0.5 * pdt + zzw = pdt + + ijb=ij_begin + ije=ij_end + if (pole_nord) ijb=ijb+iip1 + if (pole_sud) ije=ije-iip1 + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij = ijb,ije + mu(ij,l)=pbaru(ij,l) * zzpbar + ENDDO + ENDDO +c$OMP END DO NOWAIT + + ijb=ij_begin-iip1 + ije=ij_end + if (pole_nord) ijb=ij_begin + if (pole_sud) ije=ij_end-iip1 + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije + mv(ij,l)=pbarv(ij,l) * zzpbar + ENDDO + ENDDO +c$OMP END DO NOWAIT + + ijb=ij_begin + ije=ij_end + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije + mw(ij,l)=w(ij,l) * zzw + ENDDO + ENDDO +c$OMP END DO NOWAIT + +c$OMP MASTER + DO ij=ijb,ije + mw(ij,llm+1)=0. + ENDDO +c$OMP END MASTER + +c CALL SCOPY(ijp1llm,q,1,zq,1) +c CALL SCOPY(ijp1llm,masse,1,zm,1) + + ijb=ij_begin + ije=ij_end + + DO iq=1,nqmx +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + zq(ijb:ije,l,iq)=q(ijb:ije,l,iq) + zm(ijb:ije,l,iq)=masse(ijb:ije,l) + ENDDO +c$OMP END DO NOWAIT + ENDDO + + +c$OMP BARRIER + DO iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10) then + + call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu, + & ij_begin,ij_begin+2*iip1-1) + call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu, + & ij_end-2*iip1+1,ij_end) + +c$OMP MASTER + call VTb(VTHallo) + call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1) + call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1) + call VTe(VTHallo) +c$OMP END MASTER + else if (iadv(iq)==14) then + + call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat, + & ij_begin,ij_begin+2*iip1-1) + call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat, + & ij_end-2*iip1+1,ij_end) +c$OMP MASTER + call VTb(VTHallo) + call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1) + call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1) + call VTe(VTHallo) +c$OMP END MASTER + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo + + +c$OMP BARRIER +c$OMP MASTER + call VTb(VTHallo) + call SendRequest(MyRequest1) + call VTe(VTHallo) +c$OMP END MASTER +c$OMP BARRIER + do iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10) then + + call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu, + & ij_begin+2*iip1,ij_end-2*iip1) + + else if (iadv(iq)==14) then + + call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat, + & ij_begin+2*iip1,ij_end-2*iip1) + + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo +c$OMP BARRIER +c$OMP MASTER + call VTb(VTHallo) + call WaitRecvRequest(MyRequest1) + call WaitSendRequest(MyRequest1) + call VTe(VTHallo) +c$OMP END MASTER +c$OMP BARRIER + + do iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10) then + + call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv) + + else if (iadv(iq)==14) then + + call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat) + + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo + + do iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10 .or. iadv(iq)==14 ) then + +c$OMP BARRIER + call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw, + & ij_begin,ij_begin+2*iip1-1) + call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw, + & ij_end-2*iip1+1,ij_end) +c$OMP BARRIER + +c$OMP MASTER + call VTb(VTHallo) + call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2) + call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2) + call VTe(VTHallo) +c$OMP END MASTER +c$OMP BARRIER + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo +c$OMP BARRIER +c$OMP MASTER + call VTb(VTHallo) + call SendRequest(MyRequest2) + call VTe(VTHallo) +c$OMP END MASTER +c$OMP BARRIER + do iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10 .or. iadv(iq)==14 ) then +c$OMP BARRIER + call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw, + & ij_begin+2*iip1,ij_end-2*iip1) +c$OMP BARRIER + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo + +c$OMP BARRIER +c$OMP MASTER + call VTb(VTHallo) + call WaitRecvRequest(MyRequest2) + call WaitSendRequest(MyRequest2) + call VTe(VTHallo) +c$OMP END MASTER +c$OMP BARRIER + + + do iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10) then + + call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv) + + else if (iadv(iq)==14) then + + call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat) + + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo + + do iq=1,nqmx + + if(iadv(iq) == 0) then + + cycle + + else if (iadv(iq)==10) then + + call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu, + & ij_begin,ij_end) + + else if (iadv(iq)==14) then + + call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat, + & ij_begin,ij_end) + + else + + stop 'vlspltgen_p : schema non parallelise' + + endif + + enddo + + + ijb=ij_begin + ije=ij_end +c$OMP BARRIER + + DO iq=1,nqmx + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije +c print *,'zq-->',ij,l,iq,zq(ij,l,iq) +c print *,'q-->',ij,l,iq,q(ij,l,iq) + q(ij,l,iq)=zq(ij,l,iq) + ENDDO + ENDDO +c$OMP END DO NOWAIT + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO l=1,llm + DO ij=ijb,ije-iip1+1,iip1 + q(ij+iim,l,iq)=q(ij,l,iq) + ENDDO + ENDDO +c$OMP END DO NOWAIT + + ENDDO + + +c$OMP BARRIER + +cc$OMP MASTER +c call WaitSendRequest(MyRequest1) +c call WaitSendRequest(MyRequest2) +cc$OMP END MASTER +cc$OMP BARRIER + + RETURN + END Index: LMDZ4/branches/V3_test/libf/dyn3dpar/write_field_p.F90 =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/write_field_p.F90 (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/write_field_p.F90 (revision 708) @@ -0,0 +1,73 @@ +module write_field_p +implicit none + + interface WriteField_p + module procedure Write_field3d_p,Write_Field2d_p,Write_Field1d_p + end interface WriteField_p + + contains + + subroutine write_field1D_p(name,Field) + USE parallel + USE write_field + implicit none + + integer, parameter :: MaxDim=1 + character(len=*) :: name + real, dimension(:) :: Field + real, dimension(:),allocatable :: New_Field + integer, dimension(MaxDim) :: Dim + + + Dim=shape(Field) + allocate(New_Field(Dim(1))) + New_Field(:)=Field(:) + call Gather_Field(New_Field,dim(1),1,0) + + if (MPI_Rank==0) call WriteField(name,New_Field) + + end subroutine write_field1D_p + + subroutine write_field2D_p(name,Field) + USE parallel + USE write_field + implicit none + + integer, parameter :: MaxDim=2 + character(len=*) :: name + real, dimension(:,:) :: Field + real, dimension(:,:),allocatable :: New_Field + integer, dimension(MaxDim) :: Dim + + Dim=shape(Field) + allocate(New_Field(Dim(1),Dim(2))) + New_Field(:,:)=Field(:,:) + call Gather_Field(New_Field(1,1),dim(1)*dim(2),1,0) + + if (MPI_Rank==0) call WriteField(name,New_Field) + + + end subroutine write_field2D_p + + subroutine write_field3D_p(name,Field) + USE parallel + USE write_field + implicit none + + integer, parameter :: MaxDim=3 + character(len=*) :: name + real, dimension(:,:,:) :: Field + real, dimension(:,:,:),allocatable :: New_Field + integer, dimension(MaxDim) :: Dim + + Dim=shape(Field) + allocate(New_Field(Dim(1),Dim(2),Dim(3))) + New_Field(:,:,:)=Field(:,:,:) + call Gather_Field(New_Field(1,1,1),dim(1)*dim(2),dim(3),0) + + if (MPI_Rank==0) call WriteField(name,New_Field) + + end subroutine write_field3D_p + +end module write_field_p + Index: LMDZ4/branches/V3_test/libf/dyn3dpar/writedynav_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/writedynav_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/writedynav_p.F (revision 708) @@ -0,0 +1,165 @@ +! +! $Header$ +! + subroutine writedynav_p( histid, nq, time, vcov, + , ucov,teta,ppk,phi,q,masse,ps,phis) + + USE ioipsl + USE parallel + USE misc_mod + implicit none + +C +C Ecriture du fichier histoire au format IOIPSL +C +C Appels succesifs des routines: histwrite +C +C Entree: +C histid: ID du fichier histoire +C nqmx: nombre maxi de traceurs +C time: temps de l'ecriture +C vcov: vents v covariants +C ucov: vents u covariants +C teta: temperature potentielle +C phi : geopotentiel instantane +C q : traceurs +C masse: masse +C ps :pression au sol +C phis : geopotentiel au sol +C +C +C Sortie: +C fileid: ID du fichier netcdf cree +C +C L. Fairhead, LMD, 03/99 +C +C ===================================================================== +C +C Declarations +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comvert.h" +#include "comgeom.h" +#include "temps.h" +#include "ener.h" +#include "logic.h" +#include "description.h" +#include "serre.h" +#include "advtrac.h" + +C +C Arguments +C + + INTEGER histid, nq + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) + REAL teta(ip1jmp1,llm),phi(ip1jmp1,llm),ppk(ip1jmp1,llm) + REAL ps(ip1jmp1),masse(ip1jmp1,llm) + REAL phis(ip1jmp1) + REAL q(ip1jmp1,llm,nq) + integer time + + +C Variables locales +C + integer ndex2d(iip1*jjp1),ndex3d(iip1*jjp1*llm),iq, ii, ll + real us(ip1jmp1,llm), vs(ip1jmp1,llm) + real tm(ip1jmp1,llm) + REAL vnat(ip1jm,llm),unat(ip1jmp1,llm) + logical ok_sync + integer itau_w + integer :: ijb,ije,jjn +C +C Initialisations +C + if (adjust) return + + ndex3d = 0 + ndex2d = 0 + ok_sync = .TRUE. + us = 999.999 + vs = 999.999 + tm = 999.999 + vnat = 999.999 + unat = 999.999 + itau_w = itau_dyn + time + +C Passage aux composantes naturelles du vent + call covnat_p(llm, ucov, vcov, unat, vnat) + +C +C Appels a histwrite pour l'ecriture des variables a sauvegarder +C +C Vents U scalaire +C + call gr_u_scal_p(llm, unat, us) + + ijb=ij_begin + ije=ij_end + jjn=jj_nb + + call histwrite(histid, 'u', itau_w, us(ijb:ije,:), + . iip1*jjn*llm, ndex3d) +C +C Vents V scalaire +C + if (pole_sud) ije=jj_end-iip1 + if (pole_sud) jjn=jj_nb-1 + + call gr_v_scal_p(llm, vnat, vs) + call histwrite(histid, 'v', itau_w, vs(ijb::ije,:), + . iip1*jjn*llm, ndex3d) +C +C Temperature potentielle moyennee +C + ijb=ij_begin + ije=ij_end + jjn=jj_nb + + call histwrite(histid, 'theta', itau_w, teta(ijb::ije,:), + . iip1*jjn*llm, ndex3d) +C +C Temperature moyennee +C + do ll=1,llm + do ii = ijb, ije + tm(ii,ll) = teta(ii,ll) * ppk(ii,ll)/cpp + enddo + enddo + + call histwrite(histid, 'temp', itau_w, tm(ijb:ije,:), + . iip1*jjn*llm, ndex3d) +C +C Geopotentiel +C + call histwrite(histid, 'phi', itau_w, phi(ijb:ije,:), + . iip1*jjn*llm, ndex3d) +C +C Traceurs +C + DO iq=1,nq + call histwrite(histid, ttext(iq), itau_w, q(ijb:ije,:,iq), + . iip1*jjn*llm, ndex3d) + enddo +C +C Masse +C + call histwrite(histid, 'masse', itau_w, masse(ijb:ije,1), + . iip1*jjn, ndex2d) +C +C Pression au sol +C + call histwrite(histid, 'ps', itau_w, ps(ijb:ije), + . iip1*jjn, ndex2d) +C +C Geopotentiel au sol +C + call histwrite(histid, 'phis', itau_w, phis(ijb:ije), + . iip1*jjn, ndex2d) +C +C Fin +C + if (ok_sync) call histsync(histid) + return + end Index: LMDZ4/branches/V3_test/libf/dyn3dpar/writehist_p.F =================================================================== --- LMDZ4/branches/V3_test/libf/dyn3dpar/writehist_p.F (revision 708) +++ LMDZ4/branches/V3_test/libf/dyn3dpar/writehist_p.F (revision 708) @@ -0,0 +1,148 @@ +! +! $Header$ +! + subroutine writehist_p( histid, histvid, nq, time, vcov, + , ucov,teta,phi,q,masse,ps,phis) + + USE ioipsl + USE parallel + USE misc_mod + implicit none + +C +C Ecriture du fichier histoire au format IOIPSL +C +C Appels succesifs des routines: histwrite +C +C Entree: +C histid: ID du fichier histoire +C histvid:ID du fichier histoire pour les vents V (appele a disparaitre) +C nqmx: nombre maxi de traceurs +C time: temps de l'ecriture +C vcov: vents v covariants +C ucov: vents u covariants +C teta: temperature potentielle +C phi : geopotentiel instantane +C q : traceurs +C masse: masse +C ps :pression au sol +C phis : geopotentiel au sol +C +C +C Sortie: +C fileid: ID du fichier netcdf cree +C +C L. Fairhead, LMD, 03/99 +C +C ===================================================================== +C +C Declarations +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comvert.h" +#include "comgeom.h" +#include "temps.h" +#include "ener.h" +#include "logic.h" +#include "description.h" +#include "serre.h" +#include "advtrac.h" + +C +C Arguments +C + + INTEGER histid, nq, histvid + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) + REAL teta(ip1jmp1,llm),phi(ip1jmp1,llm) + REAL ps(ip1jmp1),masse(ip1jmp1,llm) + REAL phis(ip1jmp1) + REAL q(ip1jmp1,llm,nq) + integer time + + +C Variables locales +C + integer iq, ii, ll + integer ndexu(ip1jmp1*llm),ndexv(ip1jm*llm),ndex2d(ip1jmp1) + logical ok_sync + integer itau_w + integer :: ijb,ije,jjn +C +C Initialisations +C + if (adjust) return + + + ndexu = 0 + ndexv = 0 + ndex2d = 0 + ok_sync =.TRUE. + itau_w = itau_dyn + time +C +C Appels a histwrite pour l'ecriture des variables a sauvegarder +C +C Vents U +C + ijb=ij_begin + ije=ij_end + jjn=jj_nb + + call histwrite(histid, 'ucov', itau_w, ucov, + . iip1*jjn*llm, ndexu) + +C +C Vents V +C + if (pole_sud) ije=jj_end-iip1 + if (pole_sud) jjn=jj_nb-1 + + call histwrite(histvid, 'vcov', itau_w, vcov(ijb:ije,:), + . iip1*jjn*llm, ndexv) + +C +C Temperature potentielle +C + ijb=ij_begin + ije=ij_end + jjn=jj_nb + + call histwrite(histid, 'teta', itau_w, teta(ijb:ije,:), + . iip1*jjn*llm, ndexu) +C +C Geopotentiel +C + call histwrite(histid, 'phi', itau_w, phi(ijb:ije,:), + . iip1*jjn*llm, ndexu) +C +C Traceurs +C + DO iq=1,nq + call histwrite(histid, ttext(iq), itau_w, q(ijb:ije,:,iq), + . iip1*jjn*llm, ndexu) + enddo +C +C Masse +C + call histwrite(histid, 'masse', itau_w, masse(ijb:ije,1), + . iip1*jjn, ndex2d) +C +C Pression au sol +C + call histwrite(histid, 'ps', itau_w, ps(ijb:ije), + . iip1*jjn, ndex2d) +C +C Geopotentiel au sol +C + call histwrite(histid, 'phis', itau_w, phis(ijb:ije), + . iip1*jjn, ndex2d) +C +C Fin +C + if (ok_sync) then + call histsync(histid) + call histsync(histvid) + endif + return + end