Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F (revision 2324) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F (revision 2325) @@ -3,6 +3,9 @@ USE filtreg_mod, ONLY: inifilr USE comconst_mod, ONLY: dtvr,g,r,pi - + + IMPLICIT NONE + + c======================================================================= c @@ -10,7 +13,9 @@ c c======================================================================= + + c----------------------------------------------------------------------- c declarations: -c ------------- +c----------------------------------------------------------------------- @@ -19,9 +24,8 @@ include "comdissip.h" include "comgeom2.h" -!#include "control.h" include "netcdf.inc" - INTEGER itau,nbpas,nbpasmx + INTEGER itau,nbpas,nbpasmx PARAMETER(nbpasmx=1000000) REAL temps(nbpasmx) @@ -41,4 +45,14 @@ c ------------------------------------ INTEGER ivik(2),jvik(2),ifile(2),iv + + REAL, PARAMETER :: lonvik1 = -47.95 + REAL, PARAMETER :: latvik1 = 22.27 + REAL, PARAMETER :: lonvik2 = 134.29 + REAL, PARAMETER :: latvik2 = 47.67 + + REAL, PARAMETER :: phivik1 = -3637 + REAL, PARAMETER :: phivik2 = -4505 + + REAL lonvik(2),latvik(2),phivik(2),phisim(2) REAL unanj @@ -53,9 +67,9 @@ REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta - LOGICAL firstcal,lcal,latcal,lvent,day_ls + LOGICAL firstcal INTEGER*4 day0 REAL ziceco2(iip1,jjp1) - REAL day,zt,sollong,sol,dayw + REAL day,zt,sollong,sol,dayw,dayw_ls REAL airtot1,gh @@ -70,4 +84,7 @@ CHARACTER file*80 CHARACTER pathchmp*80,pathsor*80,nomfich*80 + + INTEGER Time_unit + c externe: @@ -76,7 +93,8 @@ EXTERNAL iniconst,inigeom,covcont,mywrite EXTERNAL exner,pbar - EXTERNAL solarlong,coordij,moy2 + EXTERNAL coordij,moy2 EXTERNAL SSUM REAL SSUM + cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -84,5 +102,5 @@ c----------------------------------------------------------------------- c initialisations: -c ---------------- +c----------------------------------------------------------------------- chr2="0" @@ -90,37 +108,25 @@ print*,'WARNING!!!',unanj,'Jours/an' nc=.true. - lcal=.true. - latcal=.true. - lvent=.false. - day_ls=.true. - -c lecture du fichier xvik.def - - phivik(1)=-3637 - phivik(2)=-4505 - - - - OPEN(99,file='xvik.def',form='formatted') - - READ(99,*) - READ(99,*,iostat=ierr) phivik - IF(ierr.NE.0) GOTO 105 - - READ(99,*,END=105) - READ(99,'(a)',END=105) pathchmp - READ(99,*,END=105) - READ(99,'(a)',END=105) pathsor - READ(99,*,END=105) -c READ(99,'(l1)',END=105) day_ls - READ(99,'(l1)',END=105) - READ(99,'(l1)',END=105) lcal - READ(99,'(l1)',END=105) - READ(99,'(l1)',END=105) lvent - READ(99,'(l1)',END=105) - READ(99,'(l1)',END=105) latcal - - 105 CONTINUE - CLOSE(99) + + phivik(1) = phivik1 + phivik(2) = phivik2 + + print *, 'COORDVIKIIIN', latvik, lonvik + print*, 'LES PHIVIK', phivik + + + + + + WRITE(*,*) 'Chemin des fichiers histoires' + READ (*,'(a)') pathchmp + WRITE(*,*) 'Chemin des fichiers sorties' + READ (*,'(a)') pathsor + + WRITE(*,*) 'Fichiers de sortie en sol (1) + &,en ls (2) ,les deux (3)' + READ (*,*) Time_unit + + write (*,*)'>>>>>>>>>>>>>>>>', phivik,g DO iv=1,2 @@ -128,24 +134,16 @@ END DO - write(*,*) ' pathchmp:',trim(pathchmp) - write(*,*) ' pathsor:',trim(pathsor) - -c----------------------------------------------------------------------- c----------------------------------------------------------------------- c ouverture des fichiers xgraph: -c ------------------------------ - +c----------------------------------------------------------------------- ifile(1)=12 ifile(2)=13 kyear=-1 -c OPEN(77,file='xlongday',form='formatted') - unitlec=11 -c - PRINT*,'entrer le nom du fichier NC' - READ(5,'(a)') nomfich - - PRINT *,'nomfich : ',nomfich - + + + print*,'Entrer un fichier NC (sans le .nc)' + READ(5,'(a)',err=9999) nomfich + c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -157,5 +155,7 @@ PRINT *,'>>> nomfich : ',trim(nomfich) -c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +c---------------------------------------------------------------------- +c Ouverture des fichiers histoire: +c---------------------------------------------------------------------- file=pathchmp(1:len_trim(pathchmp))//'/'// @@ -174,5 +174,5 @@ c---------------------------------------------------------------------- c initialisation de la physique: -c ------------------------------ +c---------------------------------------------------------------------- CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR) @@ -183,8 +183,9 @@ CALL iniconst CALL inigeom -! CALL inifilr - - + +c---------------------------------------------------------------------- c Lecture temps : +c---------------------------------------------------------------------- + ierr= NF_INQ_DIMID (nid,"Time",dimid) @@ -208,15 +209,21 @@ PRINT*,'temps(1:10)',(temps(itau),itau=1,10) - + + c----------------------------------------------------------------------- c coordonnees des point Viking: -c ----------------------------- - - latvik(1)=22.27*pi/180. - lonvik(1)=-47.95*pi/180. - latvik(2)=47.67*pi/180. - lonvik(2)=(360.-225.71)*pi/180. - +c -------------------------------------------------------------------- + + lonvik(1) = lonvik1 * pi/180. + latvik(1) = latvik1 * pi/180. + lonvik(2) = lonvik2 * pi/180. + latvik(2) = latvik2 * pi/180. + + +c---------------------------------------------------------------------- c ponderations pour les 4 points autour de Viking +c---------------------------------------------------------------------- + + DO iv=1,2 ! locate index of GCM grid points near VL @@ -247,5 +254,9 @@ ENDDO +c---------------------------------------------------------------------- c altitude reelle et modele aux points Viking +c---------------------------------------------------------------------- + + DO iv=1,2 phisim(iv)=0. @@ -259,8 +270,9 @@ ENDDO PRINT*,'relief aux points Viking pour les sorties:',phivik + c---------------------------------------------------------------------- c lectures des etats: -c ------------------- +c ------------------------------------------------------------------- airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1)) @@ -269,4 +281,6 @@ c debut de la boucle sur les etats dans un fichier histoire: c====================================================================== + + count_ps=(/iip1,jjp1,1/) count_co2ice=(/iip1,jjp1,1/) @@ -278,14 +292,13 @@ start_co2ice=(/1,1,itau/) start_temp=(/1,1,1,itau/) + +c---------------------------------------------------------------------- c lecture drs des champs: -c ----------------------- -c varname='u' -c ierr=drsread (unitlec,varname,unat,itau) -c PRINT*,'unat',unat(iip1/2,jjp1/2,llm/2) -c varname='v' -c ierr=drsread (unitlec,varname,vnat,itau) -c PRINT*,'vnat',vnat(iip1/2,jjp1/2,llm/2) +c---------------------------------------------------------------------- + ccccccccc LECTURE Ps ccccccccccccccccccccccccccc + + ierr = NF_INQ_VARID (nid, "ps", nvarid) #ifdef NC_DOUBLE @@ -302,4 +315,6 @@ ccccccccc LECTURE Temperature ccccccccccccccccccccccccccc + + ierr = NF_INQ_VARID (nid, "temp", nvarid) #ifdef NC_DOUBLE @@ -336,7 +351,9 @@ ENDIF -c PRINT*,'t',t(iip1/2,jjp1/2,llm/2) + ccccccccc LECTURE co2ice ccccccccccccccccccccccccccc + + ierr = NF_INQ_VARID (nid, "co2ice", nvarid) #ifdef NC_DOUBLE @@ -352,17 +369,19 @@ ENDIF - +c---------------------------------------------------------------------- c Gestion du temps -c ---------------- +c --------------------------------------------------------------------- + day=temps(itau) PRINT*,'day ',day - CALL solarlong(day+day0,sollong) - sol=day+day0+461. ! aded offset to sync with VL mission "sol 1" + sol=day+day0 iyear=sol/unanj WRITE (*,*) 'iyear',iyear sol=sol-iyear*unanj -c + +c---------------------------------------------------------------------- c Ouverture / fermeture des fichiers -c ---------------------------------- +c --------------------------------------------------------------------- + IF (iyear.NE.kyear) THEN WRITE(chr2(1:1),'(i1)') iyear+1 @@ -386,40 +405,20 @@ CLOSE(5+ifile(1)) OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted') - OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted') -c OPEN(ifile(1)+8,file='xbpsol1'//chr2,form='formatted') -c OPEN(ifile(2)+8,file='xbpsol2'//chr2,form='formatted') -c OPEN(ifile(1)+2,file='xlps1'//chr2,form='formatted') -c OPEN(ifile(2)+2,file='xlps2'//chr2,form='formatted') - IF(lcal) THEN -c OPEN(ifile(2)+4,file='xpressud'//chr2,form='formatted') -c OPEN(ifile(1)+4,file='xpresnord'//chr2,form='formatted') -c OPEN(ifile(1)+6,file='xpm2'//chr2,form='formatted') - ENDIF - IF(latcal) THEN -c OPEN(ifile(2)+14,file='xlats'//chr2,form='formatted') -c OPEN(ifile(1)+14,file='xlatn'//chr2,form='formatted') - ENDIF - IF(lvent) THEN -c OPEN(ifile(1)+16,file='xu1'//chr2,form='formatted') -c OPEN(ifile(2)+16,file='xu2'//chr2,form='formatted') -c OPEN(ifile(1)+12,file='xv1'//chr2,form='formatted') -c OPEN(ifile(2)+12,file='xv2'//chr2,form='formatted') - ENDIF + OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted') OPEN(97,file='xprestot'//chr2,form='formatted') -c OPEN(98,file='xlat37_'//chr2,form='formatted') -c WRITE(98,'(f5.1,16f7.1)') 0.,(rlonv(i)*180./pi,i=1,iim,4) + ENDIF - - sollong=sollong*180./pi - IF(day_ls) THEN - dayw=sol - write(*,*) 'dayw', dayw - ELSE - dayw=sollong - ENDIF - -c Calcul de la moyenne planetaire -c ------------------------------- + dayw = sol + call sol2ls(sol,sollong) + dayw_ls = sollong + + + +c---------------------------------------------------------------------- +c Calcul de la moyenne de pression planetaire +c --------------------------------------------------------------------- + + pstot=0. captotS=0. @@ -441,18 +440,40 @@ ENDDO ENDDO + + +c --------------Ecriture fichier sortie xprestot----------------------- +c Sol ou ls ou les deux +c Ps_moy_planetaire (Pa) +c Pequivalente de glace de CO2 au Nord (si entierement sublimee) (Pa) +c Pequivalente de glace de CO2 au Sud (si entierement sublimee) (Pa) + + + IF(Time_unit == 1) THEN WRITE(97,'(4e16.6)') dayw,pstot*airtot1 - & , captotN*g*airtot1, captotS*g*airtot1 - - IF(.NOT.firstcal) THEN -c WRITE(98,'(f5.1,16f7.3)') -c s dayw,(ps(i,37),i=1,iim,4) - + & , captotN*g*airtot1, captotS*g*airtot1 + + ELSEIF (Time_unit == 2) THEN + WRITE(97,'(4e16.6)') dayw_ls,pstot*airtot1 + & , captotN*g*airtot1, captotS*g*airtot1 + + ELSE + WRITE(97,'(5e16.6)') dayw,dayw_ls,pstot*airtot1 + & , captotN*g*airtot1,captotS*g*airtot1 + + + ENDIF + +c---------------------------------------------------------------------- c boucle sur les sites vikings: -c ---------------------------- - - DO iv=1,2 - +c---------------------------------------------------------------------- + +c---------------------------------------------------------------------- c interpolation de la temperature dans la 7eme couche, de la pression c de surface et des vents aux points viking. +c---------------------------------------------------------------------- + + IF(.NOT.firstcal) THEN + + DO iv=1,2 zp1=0. @@ -460,25 +481,28 @@ zp2_sm=0. zt=0. -! zu=0. -! zv=0. + DO jj=0,1 + j=jvik(iv)+jj + DO ii=0,1 + i=ivik(iv)+ii -! zt=zt+zw(ii,jj,iv)*t(i,j,7) zt=zt+zw(ii,jj,iv)*t7(i,j) -! zp1=zp1+zw(ii,jj,iv)*ps(i,j) zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P) - WRITE (*,*) 'ps autour iv',ps(i,j),iv -! zu=zu+zw(ii,jj,iv)*unat(i,j,1)/cu(i,j) -! zv=zv+zw(ii,jj,iv)*vnat(i,j,1)/cv(i,j) + WRITE (*,*) 'ps autour iv',ps(i,j),iv + ENDDO ENDDO + zp1=exp(zp1) ! because of the bilinear interpolation in log(P) - -c pression au sol extrapolee a partir de la temp. 7eme couche - WRITE (*,*) 'constR ',constR - WRITE (*,*) 'zt ',zt - gh=constR*zt + WRITE (*,*) 'constR ',constR + WRITE (*,*) 'zt ',zt + gh=constR*zt + +c---------------------------------------------------------------------- +c pression au sol extrapolee a partir de la temp. 7eme couche +c---------------------------------------------------------------------- + if (gh.eq.0) then ! if we don't have temperature values ! assume a scale height of 10km @@ -487,35 +511,33 @@ zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh) endif - WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh' - WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh -! WRITE(ifile(iv)+10,'(2e15.5)') dayw,zp1 - WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1 - -c sorties eventuelles de vent -! IF(lvent) THEN -! WRITE(ifile(iv)+16,'(2e15.5)') -! s dayw,zu -! WRITE(ifile(iv)+12,'(2e15.5)') -! s dayw,zv -! ENDIF + + WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh' + WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh + + +c ------Ecriture 2 fichiers (1 pour Vl1, 1 pour VL2) sortie xpsol ------ +c Sol ou ls ou les deux +c Ps site VLi (i=1,2) a  l'altitude GCM (Pa) +c Ps site VLi (i=1,2) a  l'altitude exacte (interpolee) (Pa) + + IF(Time_unit == 1) THEN + WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1 + ELSEIF (Time_unit == 2) THEN + WRITE(ifile(iv)+10,'(3e15.5)') dayw_ls,zp2,zp1 + ELSE + WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1 + ENDIF + ENDDO -c IF (lcal) THEN -c WRITE(ifile(1)+4,'(2e15.6)') dayw,airtot1*g*.01* -c s (SSUM(ip1jmp1/2,ziceco2,1)-SSUM(jjp1/2,ziceco2,iip1)) -c WRITE(ifile(2)+4,'(2e15.6)') dayw,airtot1*g*.01* -c s (SSUM(iip1*jjm/2,ziceco2(1,jjm/2+2),1)- -c s SSUM(jjm/2,ziceco2(1,jjm/2+2),iip1)) -c ENDIF -c IF(latcal) THEN -c CALL icelat(iim,jjm,ziceco2,rlatv,zicelat) -c WRITE(ifile(1)+14,'(2e15.6)') dayw,zicelat(1)*180./pi -c WRITE(ifile(2)+14,'(2e15.6)') dayw,zicelat(2)*180./pi -c ENDIF + ENDIF + firstcal=.false. + c====================================================================== c Fin de la boucle sur les etats du fichier histoire: c====================================================================== + ENDDO @@ -523,6 +545,8 @@ PRINT*,'Fin du fichier',nomfich - print*,'Entrer un nouveau fichier ou return pour finir' + print*,'Entrer un nouveau fichier NC + &(sans le .nc) ou return pour finir' READ(5,'(a)',err=9999) nomfich + c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -552,3 +576,99 @@ 7777 FORMAT ('latitude/longitude',4f7.1) + + + END + + subroutine sol2ls(sol,Ls) +!============================================================================== +! Purpose: +! Convert a date/time, given in sol (martian day), +! into solar longitude date/time, in Ls (in degrees), +! where sol=0 is (by definition) the northern hemisphere +! spring equinox (where Ls=0). +!============================================================================== +! Notes: +! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year, +! "Ls" will be increased by N*360 +! Won't work as expected if sol is negative (then again, +! why would that ever happen?) +!============================================================================== + + implicit none + +!============================================================================== +! Arguments: +!============================================================================== + real,intent(in) :: sol + real,intent(out) :: Ls + +!============================================================================== +! Local variables: +!============================================================================== + real year_day,peri_day,timeperi,e_elips,twopi,degrad + data year_day /669./ ! # of sols in a martian year + data peri_day /485.0/ + data timeperi /1.9082314/ + data e_elips /0.093358/ + data twopi /6.2831853/ ! 2.*pi + data degrad /57.2957795/ ! pi/180 + + real zanom,xref,zx0,zdx,zteta,zz + + integer count_years + integer iter + +!============================================================================== +! 1. Compute Ls +!============================================================================== + + zz=(sol-peri_day)/year_day + zanom=twopi*(zz-nint(zz)) + xref=abs(zanom) + +! The equation zx0 - e * sin (zx0) = xref, solved by Newton + zx0=xref+e_elips*sin(xref) + do iter=1,20 ! typically, 2 or 3 iterations are enough + zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) + zx0=zx0+zdx + if(abs(zdx).le.(1.e-7)) then +! write(*,*)'iter:',iter,' |zdx|:',abs(zdx) + exit + endif + enddo + + if(zanom.lt.0.) zx0=-zx0 + + zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) + Ls=zteta-timeperi + + if(Ls.lt.0.) then + Ls=Ls+twopi + else + if(Ls.gt.twopi) then + Ls=Ls-twopi + endif + endif + + Ls=degrad*Ls +! Ls is now in degrees + +!============================================================================== +! 1. Account for (eventual) years included in input date/time sol +!============================================================================== + +count_years=0 ! initialize + zz=sol ! use "zz" to store (and work on) the value of sol + do while (zz.ge.year_day) + count_years=count_years+1 + zz=zz-year_day + enddo + +! Add 360 degrees to Ls for every year + if (count_years.ne.0) then + Ls=Ls+360.*count_years + endif + + + end subroutine sol2ls