Changeset 2325


Ignore:
Timestamp:
May 18, 2020, 1:36:01 PM (5 years ago)
Author:
tpierron
Message:

Mars GCM :

  • Update of the program xvik.F
  • Useless code removed
  • Add some comments about xvik outputs
  • Declaration of physics constants inside the program as variable and not as inputs like before
  • Temporal shift due to viking landing sol removed : xvik outputs are given in real sol
  • Add the possibility to choose the format of time dimension for xvik outputs : ls, sol or both
  • Path to input and output files removed from the code and put as input
  • List of input files removed from code and put as input

TP

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2324 r2325  
    30493049-simpleclouds.F,hdo_surfex_mod.F,physiq_mod.F: thresholds of 1.e-16 have been replaced by the variable qperemin defined in tracer_mod.F.
    30503050
     3051== 18/05/2020 == TP
     3052- Update of the program xvik.F
     3053- Useless code removed
     3054- Add some comments about xvik outputs
     3055- Declaration of physics constants inside the program as variable and not as inputs like before
     3056- Temporal shift due to viking landing sol removed : xvik outputs are given in real sol
     3057- Add the possibility to choose the format of time dimension for xvik outputs : ls, sol or both
     3058- Path to input and output files removed from the code and put as input
     3059- List of input files removed from code and put as input
     3060

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/xvik.F

    r1977 r2325  
    33      USE filtreg_mod, ONLY: inifilr
    44      USE comconst_mod, ONLY: dtvr,g,r,pi
    5 
     5     
     6     
    67      IMPLICIT NONE
     8     
     9     
    710c=======================================================================
    811c
    … …  
    1013c
    1114c=======================================================================
     15
     16
    1217c-----------------------------------------------------------------------
    1318c   declarations:
    14 c   -------------
     19c-----------------------------------------------------------------------
    1520
    1621
    … …  
    1924      include "comdissip.h"
    2025      include "comgeom2.h"
    21 !#include "control.h"
    2226      include "netcdf.inc"     
    2327
    2428
    25       INTEGER itau,nbpas,nbpasmx
     29      INTEGER itau,nbpas,nbpasmx 
    2630      PARAMETER(nbpasmx=1000000)
    2731      REAL temps(nbpasmx)
    … …  
    4145c   ------------------------------------
    4246      INTEGER ivik(2),jvik(2),ifile(2),iv
     47     
     48      REAL, PARAMETER ::  lonvik1 = -47.95
     49      REAL, PARAMETER ::  latvik1 =  22.27
     50      REAL, PARAMETER ::  lonvik2 =  134.29
     51      REAL, PARAMETER ::  latvik2 =  47.67
     52     
     53      REAL, PARAMETER :: phivik1 = -3637
     54      REAL, PARAMETER :: phivik2 = -4505
     55     
     56     
    4357      REAL lonvik(2),latvik(2),phivik(2),phisim(2)
    4458      REAL unanj
    … …  
    5367      REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta
    5468
    55       LOGICAL firstcal,lcal,latcal,lvent,day_ls
     69      LOGICAL firstcal
    5670      INTEGER*4 day0
    5771
    5872      REAL ziceco2(iip1,jjp1)
    59       REAL day,zt,sollong,sol,dayw
     73      REAL day,zt,sollong,sol,dayw,dayw_ls
    6074      REAL airtot1,gh
    6175
    … …  
    7084      CHARACTER file*80
    7185      CHARACTER pathchmp*80,pathsor*80,nomfich*80
     86     
     87      INTEGER Time_unit
     88     
    7289
    7390c   externe:
    … …  
    7693      EXTERNAL iniconst,inigeom,covcont,mywrite
    7794      EXTERNAL exner,pbar
    78       EXTERNAL solarlong,coordij,moy2
     95      EXTERNAL coordij,moy2
    7996      EXTERNAL SSUM
    8097      REAL SSUM
     98     
    8199
    82100cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    84102c-----------------------------------------------------------------------
    85103c   initialisations:
    86 c   ----------------
     104c-----------------------------------------------------------------------
    87105
    88106      chr2="0"
    … …  
    90108      print*,'WARNING!!!',unanj,'Jours/an'
    91109      nc=.true.
    92       lcal=.true.
    93       latcal=.true.
    94       lvent=.false.
    95       day_ls=.true.
    96 
    97 c lecture du fichier xvik.def
    98 
    99       phivik(1)=-3637
    100       phivik(2)=-4505
    101 
    102 
    103 
    104       OPEN(99,file='xvik.def',form='formatted')
    105 
    106       READ(99,*)
    107       READ(99,*,iostat=ierr) phivik
    108       IF(ierr.NE.0) GOTO 105
    109 
    110       READ(99,*,END=105)
    111       READ(99,'(a)',END=105) pathchmp
    112       READ(99,*,END=105)
    113       READ(99,'(a)',END=105) pathsor
    114       READ(99,*,END=105)
    115 c     READ(99,'(l1)',END=105) day_ls
    116       READ(99,'(l1)',END=105)
    117       READ(99,'(l1)',END=105) lcal
    118       READ(99,'(l1)',END=105)
    119       READ(99,'(l1)',END=105) lvent
    120       READ(99,'(l1)',END=105)
    121       READ(99,'(l1)',END=105) latcal
    122  
    123  105  CONTINUE
    124       CLOSE(99)
     110     
     111      phivik(1) = phivik1
     112      phivik(2) = phivik2
     113     
     114      print *, 'COORDVIKIIIN', latvik, lonvik
     115      print*, 'LES PHIVIK', phivik
     116     
     117     
     118
     119
     120
     121      WRITE(*,*) 'Chemin des fichiers histoires'
     122      READ (*,'(a)')  pathchmp
     123      WRITE(*,*) 'Chemin des fichiers sorties'
     124      READ (*,'(a)')  pathsor
     125     
     126      WRITE(*,*) 'Fichiers de sortie en sol (1)
     127     &,en ls (2) ,les deux (3)'
     128      READ (*,*)  Time_unit
     129     
     130     
    125131      write (*,*)'>>>>>>>>>>>>>>>>', phivik,g
    126132      DO iv=1,2
    … …  
    128134      END DO
    129135
    130       write(*,*) ' pathchmp:',trim(pathchmp)
    131       write(*,*) ' pathsor:',trim(pathsor)
    132      
    133 c-----------------------------------------------------------------------
    134136c-----------------------------------------------------------------------
    135137c   ouverture des fichiers xgraph:
    136 c   ------------------------------
    137 
     138c-----------------------------------------------------------------------
    138139      ifile(1)=12
    139140      ifile(2)=13
    140141      kyear=-1
    141 c      OPEN(77,file='xlongday',form='formatted')
    142 
    143142      unitlec=11
    144 c
    145       PRINT*,'entrer le nom du fichier NC'
    146       READ(5,'(a)') nomfich
    147 
    148       PRINT *,'nomfich : ',nomfich
    149 
     143     
     144     
     145      print*,'Entrer un fichier NC (sans le .nc)'
     146      READ(5,'(a)',err=9999) nomfich
     147     
    150148
    151149c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    157155      PRINT *,'>>>  nomfich : ',trim(nomfich)
    158156
    159 c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     157c----------------------------------------------------------------------
     158c   Ouverture des fichiers histoire:
     159c----------------------------------------------------------------------
    160160
    161161      file=pathchmp(1:len_trim(pathchmp))//'/'//
    … …  
    174174c----------------------------------------------------------------------
    175175c   initialisation de la physique:
    176 c   ------------------------------
     176c----------------------------------------------------------------------
    177177
    178178      CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR)
    … …  
    183183      CALL iniconst
    184184      CALL inigeom
    185 !      CALL inifilr
    186 
    187 
     185
     186c----------------------------------------------------------------------
    188187c   Lecture temps :
     188c----------------------------------------------------------------------
     189
    189190
    190191      ierr= NF_INQ_DIMID (nid,"Time",dimid)
    … …  
    208209
    209210        PRINT*,'temps(1:10)',(temps(itau),itau=1,10)
    210              
     211       
     212       
    211213c-----------------------------------------------------------------------
    212214c   coordonnees des point Viking:
    213 c   -----------------------------
    214 
    215       latvik(1)=22.27*pi/180.
    216       lonvik(1)=-47.95*pi/180.
    217       latvik(2)=47.67*pi/180.
    218       lonvik(2)=(360.-225.71)*pi/180.
    219 
     215c   --------------------------------------------------------------------
     216
     217      lonvik(1) = lonvik1 * pi/180.
     218      latvik(1) = latvik1 * pi/180.
     219      lonvik(2) = lonvik2 * pi/180.
     220      latvik(2) = latvik2 * pi/180.
     221     
     222                   
     223c----------------------------------------------------------------------   
    220224c   ponderations pour les 4 points autour de Viking
     225c----------------------------------------------------------------------
     226
     227
    221228      DO iv=1,2
    222229        ! locate index of GCM grid points near VL
    … …  
    247254      ENDDO
    248255
     256c----------------------------------------------------------------------
    249257c   altitude reelle et modele aux points Viking
     258c----------------------------------------------------------------------
     259
     260
    250261      DO iv=1,2
    251262         phisim(iv)=0.
    … …  
    259270      ENDDO
    260271      PRINT*,'relief aux points Viking pour les sorties:',phivik
     272           
    261273
    262274c----------------------------------------------------------------------
    263275c   lectures des etats:
    264 c   -------------------
     276c   -------------------------------------------------------------------
    265277
    266278       airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1))
    … …  
    269281c   debut de la boucle sur les etats dans un fichier histoire:
    270282c======================================================================
     283
     284
    271285       count_ps=(/iip1,jjp1,1/)
    272286       count_co2ice=(/iip1,jjp1,1/)
    … …  
    278292       start_co2ice=(/1,1,itau/)
    279293       start_temp=(/1,1,1,itau/)
     294       
     295c----------------------------------------------------------------------       
    280296c   lecture drs des champs:
    281 c   -----------------------
    282 c         varname='u'
    283 c         ierr=drsread (unitlec,varname,unat,itau)
    284 c         PRINT*,'unat',unat(iip1/2,jjp1/2,llm/2)
    285 c         varname='v'
    286 c         ierr=drsread (unitlec,varname,vnat,itau)
    287 c         PRINT*,'vnat',vnat(iip1/2,jjp1/2,llm/2)
     297c----------------------------------------------------------------------
     298
    288299
    289300ccccccccc  LECTURE Ps ccccccccccccccccccccccccccc
     301
     302
    290303          ierr = NF_INQ_VARID (nid, "ps", nvarid)
    291304#ifdef NC_DOUBLE
    … …  
    302315
    303316ccccccccc  LECTURE Temperature ccccccccccccccccccccccccccc
     317
     318
    304319          ierr = NF_INQ_VARID (nid, "temp", nvarid)
    305320#ifdef NC_DOUBLE
    … …  
    336351          ENDIF
    337352
    338 c          PRINT*,'t',t(iip1/2,jjp1/2,llm/2)
     353
    339354
    340355ccccccccc  LECTURE co2ice ccccccccccccccccccccccccccc
     356
     357
    341358          ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
    342359#ifdef NC_DOUBLE
    … …  
    352369          ENDIF
    353370
    354 
     371c----------------------------------------------------------------------
    355372c Gestion du temps
    356 c ----------------
     373c ---------------------------------------------------------------------
     374
    357375          day=temps(itau)
    358376          PRINT*,'day ',day
    359           CALL solarlong(day+day0,sollong)
    360           sol=day+day0+461. ! aded offset to sync with VL mission "sol 1"
     377          sol=day+day0
    361378          iyear=sol/unanj
    362379          WRITE (*,*) 'iyear',iyear
    363380          sol=sol-iyear*unanj
    364 c
     381
     382c----------------------------------------------------------------------
    365383c Ouverture / fermeture des fichiers
    366 c ----------------------------------
     384c ---------------------------------------------------------------------
     385
    367386          IF (iyear.NE.kyear) THEN
    368387             WRITE(chr2(1:1),'(i1)') iyear+1
    … …  
    386405             CLOSE(5+ifile(1))
    387406             OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted')
    388              OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted')
    389 c            OPEN(ifile(1)+8,file='xbpsol1'//chr2,form='formatted')
    390 c            OPEN(ifile(2)+8,file='xbpsol2'//chr2,form='formatted')
    391 c            OPEN(ifile(1)+2,file='xlps1'//chr2,form='formatted')
    392 c            OPEN(ifile(2)+2,file='xlps2'//chr2,form='formatted')
    393              IF(lcal) THEN
    394 c               OPEN(ifile(2)+4,file='xpressud'//chr2,form='formatted')
    395 c               OPEN(ifile(1)+4,file='xpresnord'//chr2,form='formatted')
    396 c               OPEN(ifile(1)+6,file='xpm2'//chr2,form='formatted')
    397              ENDIF
    398                          IF(latcal) THEN
    399 c               OPEN(ifile(2)+14,file='xlats'//chr2,form='formatted')
    400 c               OPEN(ifile(1)+14,file='xlatn'//chr2,form='formatted')
    401                          ENDIF
    402              IF(lvent) THEN
    403 c               OPEN(ifile(1)+16,file='xu1'//chr2,form='formatted')
    404 c               OPEN(ifile(2)+16,file='xu2'//chr2,form='formatted')
    405 c               OPEN(ifile(1)+12,file='xv1'//chr2,form='formatted')
    406 c               OPEN(ifile(2)+12,file='xv2'//chr2,form='formatted')
    407              ENDIF
     407             OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted')                                 
    408408             OPEN(97,file='xprestot'//chr2,form='formatted')
    409 c            OPEN(98,file='xlat37_'//chr2,form='formatted')
    410 c           WRITE(98,'(f5.1,16f7.1)') 0.,(rlonv(i)*180./pi,i=1,iim,4)
     409
    411410          ENDIF
    412411 
    413 
    414           sollong=sollong*180./pi
    415           IF(day_ls) THEN
    416              dayw=sol
    417              write(*,*) 'dayw', dayw
    418           ELSE
    419              dayw=sollong
    420           ENDIF
    421 
    422 c Calcul de la moyenne planetaire
    423 c -------------------------------
     412          dayw = sol
     413          call sol2ls(sol,sollong)
     414          dayw_ls = sollong
     415         
     416         
     417         
     418c----------------------------------------------------------------------
     419c Calcul de la moyenne de pression planetaire
     420c ---------------------------------------------------------------------
     421
     422
    424423          pstot=0.
    425424          captotS=0.
    … …  
    441440                 ENDDO
    442441              ENDDO
     442
     443
     444c --------------Ecriture fichier sortie xprestot-----------------------
     445c  Sol ou ls ou les deux
     446c  Ps_moy_planetaire (Pa)
     447c  Pequivalente de glace de CO2 au Nord (si entierement sublimee) (Pa)
     448c  Pequivalente de glace de CO2 au Sud (si entierement sublimee) (Pa)
     449
     450
     451         IF(Time_unit == 1) THEN
    443452              WRITE(97,'(4e16.6)') dayw,pstot*airtot1
    444      &       , captotN*g*airtot1, captotS*g*airtot1         
    445 
    446           IF(.NOT.firstcal) THEN
    447 c         WRITE(98,'(f5.1,16f7.3)')
    448 c     s    dayw,(ps(i,37),i=1,iim,4)
    449 
     453     &       , captotN*g*airtot1, captotS*g*airtot1       
     454 
     455         ELSEIF (Time_unit == 2) THEN   
     456              WRITE(97,'(4e16.6)') dayw_ls,pstot*airtot1
     457     &       , captotN*g*airtot1, captotS*g*airtot1
     458     
     459         ELSE
     460             WRITE(97,'(5e16.6)') dayw,dayw_ls,pstot*airtot1
     461     &       , captotN*g*airtot1,captotS*g*airtot1
     462     
     463                   
     464         ENDIF           
     465
     466c----------------------------------------------------------------------
    450467c boucle sur les sites vikings:
    451 c ----------------------------
    452 
    453           DO iv=1,2
    454 
     468c----------------------------------------------------------------------
     469
     470c----------------------------------------------------------------------
    455471c interpolation de la temperature dans la 7eme couche, de la pression
    456472c de surface et des vents aux points viking.
     473c----------------------------------------------------------------------
     474
     475         IF(.NOT.firstcal) THEN
     476         
     477          DO iv=1,2
    457478
    458479             zp1=0.
    … …  
    460481             zp2_sm=0.
    461482             zt=0.
    462 !             zu=0.
    463 !             zv=0.
     483
    464484             DO jj=0,1
     485             
    465486                j=jvik(iv)+jj
     487               
    466488                DO ii=0,1
     489               
    467490                   i=ivik(iv)+ii
    468 !                   zt=zt+zw(ii,jj,iv)*t(i,j,7)
    469491                   zt=zt+zw(ii,jj,iv)*t7(i,j)
    470 !                   zp1=zp1+zw(ii,jj,iv)*ps(i,j)
    471492                   zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P)
    472                     WRITE (*,*) 'ps autour iv',ps(i,j),iv
    473 !                   zu=zu+zw(ii,jj,iv)*unat(i,j,1)/cu(i,j)
    474 !                   zv=zv+zw(ii,jj,iv)*vnat(i,j,1)/cv(i,j)
     493                   WRITE (*,*) 'ps autour iv',ps(i,j),iv
     494
    475495                ENDDO
    476496             ENDDO
     497             
    477498             zp1=exp(zp1) ! because of the bilinear interpolation in log(P)
    478  
    479 c               pression au sol extrapolee a partir de la temp. 7eme couche
    480            WRITE (*,*) 'constR ',constR
    481            WRITE (*,*) 'zt ',zt
    482              gh=constR*zt
     499             WRITE (*,*) 'constR ',constR
     500             WRITE (*,*) 'zt ',zt
     501             gh=constR*zt           
     502             
     503c----------------------------------------------------------------------
     504c  pression au sol extrapolee a partir de la temp. 7eme couche
     505c----------------------------------------------------------------------
     506           
    483507             if (gh.eq.0) then ! if we don't have temperature values
    484508               ! assume a scale height of 10km
    … …  
    487511               zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh)
    488512             endif
    489            WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh'
    490            WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh
    491 !           WRITE(ifile(iv)+10,'(2e15.5)') dayw,zp1
    492            WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1
    493            
    494 c   sorties eventuelles de vent
    495 !             IF(lvent) THEN
    496 !                WRITE(ifile(iv)+16,'(2e15.5)')
    497 !     s          dayw,zu
    498 !                WRITE(ifile(iv)+12,'(2e15.5)')
    499 !     s          dayw,zv
    500 !             ENDIF
     513           
     514          WRITE (*,*) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh'
     515          WRITE (*,*) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh
     516             
     517
     518c ------Ecriture 2 fichiers (1 pour Vl1, 1 pour VL2) sortie xpsol ------
     519c  Sol ou ls ou les deux
     520c  Ps site VLi (i=1,2) a  l'altitude GCM (Pa)
     521c  Ps site VLi (i=1,2) a  l'altitude exacte  (interpolee) (Pa)
     522             
     523             IF(Time_unit == 1) THEN
     524                WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1
     525             ELSEIF (Time_unit == 2) THEN   
     526                WRITE(ifile(iv)+10,'(3e15.5)') dayw_ls,zp2,zp1
     527             ELSE   
     528                WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1
     529             ENDIF     
     530             
    501531          ENDDO
    502 c         IF (lcal) THEN
    503 c            WRITE(ifile(1)+4,'(2e15.6)') dayw,airtot1*g*.01*
    504 c    s       (SSUM(ip1jmp1/2,ziceco2,1)-SSUM(jjp1/2,ziceco2,iip1))
    505 c            WRITE(ifile(2)+4,'(2e15.6)') dayw,airtot1*g*.01*
    506 c    s       (SSUM(iip1*jjm/2,ziceco2(1,jjm/2+2),1)-
    507 c    s       SSUM(jjm/2,ziceco2(1,jjm/2+2),iip1))
    508 c         ENDIF
    509 c            IF(latcal) THEN
    510 c               CALL icelat(iim,jjm,ziceco2,rlatv,zicelat)
    511 c               WRITE(ifile(1)+14,'(2e15.6)') dayw,zicelat(1)*180./pi
    512 c               WRITE(ifile(2)+14,'(2e15.6)') dayw,zicelat(2)*180./pi
    513 c            ENDIF
     532
    514533         ENDIF
     534         
    515535         firstcal=.false.
     536
    516537
    517538c======================================================================
    518539c   Fin de la boucle sur les etats du fichier histoire:
    519540c======================================================================
     541
    520542      ENDDO
    521543
    … …  
    523545
    524546      PRINT*,'Fin du fichier',nomfich
    525       print*,'Entrer un nouveau fichier ou return pour finir'
     547      print*,'Entrer un nouveau fichier NC
     548     &(sans le .nc) ou return pour finir'
    526549      READ(5,'(a)',err=9999) nomfich
     550
    527551
    528552c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    552576
    5535777777  FORMAT ('latitude/longitude',4f7.1)
     578
     579
     580
    554581      END
     582
     583      subroutine sol2ls(sol,Ls)
     584!==============================================================================
     585! Purpose:
     586! Convert a date/time, given in sol (martian day),
     587! into solar longitude date/time, in Ls (in degrees),
     588! where sol=0 is (by definition) the northern hemisphere
     589!  spring equinox (where Ls=0).
     590!==============================================================================
     591! Notes:
     592! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
     593! "Ls" will be increased by N*360
     594! Won't work as expected if sol is negative (then again,
     595! why would that ever happen?)
     596!==============================================================================
     597
     598      implicit none
     599
     600!==============================================================================
     601! Arguments:
     602!==============================================================================
     603      real,intent(in) :: sol
     604      real,intent(out) :: Ls
     605
     606!==============================================================================
     607! Local variables:
     608!==============================================================================
     609      real year_day,peri_day,timeperi,e_elips,twopi,degrad
     610      data year_day /669./            ! # of sols in a martian year
     611      data peri_day /485.0/           
     612      data timeperi /1.9082314/
     613      data e_elips  /0.093358/
     614      data twopi       /6.2831853/    ! 2.*pi
     615      data degrad   /57.2957795/      ! pi/180
     616
     617      real zanom,xref,zx0,zdx,zteta,zz
     618
     619      integer count_years
     620      integer iter
     621
     622!==============================================================================
     623! 1. Compute Ls
     624!==============================================================================
     625
     626      zz=(sol-peri_day)/year_day
     627      zanom=twopi*(zz-nint(zz))
     628      xref=abs(zanom)
     629
     630!  The equation zx0 - e * sin (zx0) = xref, solved by Newton
     631      zx0=xref+e_elips*sin(xref)
     632      do iter=1,20 ! typically, 2 or 3 iterations are enough
     633         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
     634         zx0=zx0+zdx
     635         if(abs(zdx).le.(1.e-7)) then
     636!            write(*,*)'iter:',iter,'     |zdx|:',abs(zdx)
     637             exit
     638         endif
     639      enddo
     640
     641      if(zanom.lt.0.) zx0=-zx0
     642
     643      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
     644      Ls=zteta-timeperi
     645
     646      if(Ls.lt.0.) then
     647         Ls=Ls+twopi
     648      else
     649         if(Ls.gt.twopi) then
     650            Ls=Ls-twopi
     651         endif
     652      endif
     653
     654      Ls=degrad*Ls
     655! Ls is now in degrees
     656
     657!==============================================================================
     658! 1. Account for (eventual) years included in input date/time sol
     659!==============================================================================
     660
     661count_years=0 ! initialize
     662      zz=sol  ! use "zz" to store (and work on) the value of sol
     663      do while (zz.ge.year_day)
     664          count_years=count_years+1
     665          zz=zz-year_day
     666      enddo
     667
     668! Add 360 degrees to Ls for every year
     669      if (count_years.ne.0) then
     670         Ls=Ls+360.*count_years
     671      endif
     672
     673
     674      end subroutine sol2ls
Note: See TracChangeset for help on using the changeset viewer.