Changeset 3438 for trunk


Ignore:
Timestamp:
Sep 25, 2024, 4:43:31 PM (3 months ago)
Author:
afalco
Message:

Pluto PCM: Minor fixes.
AF

Location:
trunk/LMDZ.PLUTO
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

    r3390 r3438  
    14541454     &              (lonfi(ig)*180./pi.ge.val4)  .and.
    14551455     &              (lonfi(ig)*180./pi.lt.val5)  .and.
    1456      &              (qsurf(ig,igcm_n2).lt.val7))  then     
     1456     &              (qsurf(ig,igcm_n2).gt.val7))  then     
    14571457                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
    14581458                    qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6
  • trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90

    r3421 r3438  
    171171
    172172     ! calculate global mean surface pressure for the fast mode
     173     IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
     174     DO ig=1,klon
     175        kp(ig) = exp(-phisfi(ig)/(r*38.))
     176     ENDDO
    173177     IF (fast) THEN
    174         IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
    175         DO ig=1,klon
    176            kp(ig) = exp(-phisfi(ig)/(r*38.))
    177         ENDDO
    178178        p00=glob_average2d(kp) ! mean pres at ref level
    179179     ENDIF
     
    382382        pdpsrf(ig)   = -pdicen2(ig)*g
    383383    !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
    384         IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
     384        IF (fast.and.(zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001)) then
    385385            pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
    386386            pdicen2(ig)=-pdpsrf(ig)/g
     
    855855         Mtot = masse(m)
    856856         MQtot = masse(m)*q(m)
    857          !if (m.lt.klev) then ! because some compilers will have problems
    858          !                      ! evaluating masse(klev+1)
     857         if (m.lt.klev) then ! because some compilers will have problems
     858                              ! evaluating masse(klev+1)
    859859         do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
    860860            m=m+1
    861861            Mtot = Mtot + masse(m)
    862862            MQtot = MQtot + masse(m)*q(m)
    863          !   if (m.eq.klev) exit
     863           if (m.eq.klev) exit
    864864         end do
    865          !endif
     865         endif
    866866         if (m.lt.klev) then
    867867            sigw=(w(l+1)-Mtot)/masse(m+1)
  • trunk/LMDZ.PLUTO/util/zrecast_plut.F90

    r3353 r3438  
    4040!              'phisinit.nc')
    4141! EM 02/2007 : Changed behavior for "altitude above surface" case
    42 !              (for MCD RMS computation). Number of levels is then set as 
     42!              (for MCD RMS computation). Number of levels is then set as
    4343!              number of levels in initial file,
    4444!              and the new set of above surface levels follow a more elaborate
     
    9999real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
    100100real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
    101 real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure 
     101real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
    102102real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
    103103real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
     
    195195    if (tmpndims.eq.4) then
    196196      nbvar=nbvar+1
    197       var(nbvar)=tmpvarname   
     197      var(nbvar)=tmpvarname
    198198    else
    199199      write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
     
    490490endif
    491491! look for 'phisinit' in current file
    492 ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid) 
     492ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
    493493if (ierr.ne.NF_NOERR) then
    494494  write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
     
    598598    if (nblev.eq.1) then ! in case only one level is asked for
    599599      zareoid(nblev)=zamin
    600     else 
     600    else
    601601      do i=1,nblev
    602602        zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
     
    625625    if (nblev.eq.1) then ! in case only one level is asked for
    626626      zradius(nblev)=zamin
    627     else 
     627    else
    628628      do i=1,nblev
    629629        zradius(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
     
    977977! 3.3.1 Define 1D variables
    978978
    979 ! longitude 
     979! longitude
    980980datashape(1)=lon_dimid
    981981ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
     
    11561156  stop "Error: Problem writing long_name for Time"
    11571157endif
    1158 text='days since 0000-01-1 00:00:00'
     1158!text='days since 0000-01-1 00:00:00'
     1159text='days since 00000'
    11591160ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
    11601161if (ierr.ne.NF_NOERR) then
     
    12591260    write(*,*) 'Error, could not define variable ',trim(var(i))
    12601261    write(*,*) NF_STRERROR(ierr)
    1261     stop 
     1262    stop
    12621263  endif
    12631264
     
    13201321    endif
    13211322  endif
    1322  
     1323
    13231324  ! look for a "units" attribute
    13241325  text='   '
     
    13341335    endif
    13351336  endif
    1336  
     1337
    13371338  ! look for a "missing_value" attribute
    13381339  ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
     
    13431344    miss_val=miss_val_def
    13441345  endif
    1345  
     1346
    13461347  ! write the missing_value attribute to output
    13471348  ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
     
    13561357    write(*,*)'nbattr:',nbattr,' j:',j
    13571358  endif
    1358    
     1359
    13591360enddo ! of do=1,nbvar
    13601361
     
    14991500    stop
    15001501  endif
    1501  
     1502
    15021503  ! interpolate dataset onto new grid
    15031504  call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
    15041505                      miss_val,ps,press,indata,plevel,outdata)
    1505  
     1506
    15061507  ! write new dataset to output file
    15071508  ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
     
    15111512    stop
    15121513  endif
    1513  
     1514
    15141515  enddo ! of do i=1,nbvar
    1515  
     1516
    15161517  ! interpolate zareoid onto new grid
    15171518  call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
     
    15401541    stop
    15411542  endif
    1542  
     1543
    15431544  ! interpolate dataset onto new grid
    15441545  ! check if variable is "rho" (to set flag for interpolation below)
     
    15481549    j=0
    15491550  endif
    1550  
     1551
    15511552  call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
    15521553                      miss_val,za_gcm,indata,j,zareoid,outdata)
     
    15611562
    15621563  enddo ! of do i=1,nbvar
    1563  
     1564
    15641565  ! interpolate pressure onto new grid
    15651566  write(*,*) ""
     
    15931594      stop
    15941595    endif
    1595  
     1596
    15961597    ! interpolate dataset onto new grid
    15971598    ! check if variable is "rho" (to set flag for interpolation below)
     
    16011602      j=0
    16021603    endif
    1603  
     1604
    16041605    call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
    16051606                        miss_val,zs_gcm,indata,phisinit,press,temp,rho,&
     
    16141615    endif
    16151616  enddo ! of do i=1,nbvar
    1616  
     1617
    16171618  ! interpolate pressure onto new grid
    16181619  write(*,*) ""
     
    16461647    stop
    16471648  endif
    1648  
     1649
    16491650  ! interpolate dataset onto new grid
    16501651  ! check if variable is "rho" (to set flag for interpolation below)
     
    16541655    j=0
    16551656  endif
    1656  
     1657
    16571658  call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
    16581659                      miss_val,ra_gcm,indata,j,zradius,outdata)
     
    16671668
    16681669  enddo ! of do i=1,nbvar
    1669  
     1670
    16701671  ! interpolate pressure onto new grid
    16711672  write(*,*) ""
     
    17341735
    17351736! Parameters needed to integrate hydrostatic equation:
    1736 real,parameter :: g0=0.62
     1737real,parameter :: g0=0.6198
    17371738!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
    1738 real,parameter :: a0=1187.E3
     1739real,parameter :: a0=1188.E3
    17391740!a0: 'mean' Mars radius=3396.km
    1740 real :: gz 
     1741real :: gz
    17411742! gz: gravitational acceleration at a given (zareoid) altitude
    17421743
     
    17491750
    17501751!==============================================================================
    1751 ! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
     1752! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
    17521753!   --later used to integrate the hydrostatic equation--
    17531754!==============================================================================
     
    17791780        ! compute sigma level of layer
    17801781        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
    1781        
     1782
    17821783        ! compute "mean" temperature of the layer
    17831784        if(temp(iloop,jloop,kloop,tloop).eq.  &
     
    17901791                     temp(iloop,jloop,kloop-1,tloop))
    17911792        endif
    1792        
     1793
    17931794        ! compute mean value of R of the layer
    17941795        Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
    1795        
     1796
    17961797        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
    17971798        ! NB: zareoid=zsurface+phis/g0
    17981799        gz=g0*(a0**2)/ &
    17991800           (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
    1800        
     1801
    18011802        ! compute zs_gcm(iloop,jloop,kloop,tloop)
    18021803        zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
    18031804                log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
    1804        
     1805
    18051806        ! compute za_gcm(kloop)
    18061807!        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
     
    18601861
    18611862! Parameters needed to integrate hydrostatic equation:
    1862 real,parameter :: g0=0.62
     1863real,parameter :: g0=0.6198
    18631864!g0: exact mean gravity
    1864 real,parameter :: a0=1187.E3
     1865real,parameter :: a0=1188.E3
    18651866!a0: 'mean' radius
    1866 real :: gz 
     1867real :: gz
    18671868! gz: gravitational acceleration at a given (zareoid) altitude
    18681869
     
    18751876
    18761877!==============================================================================
    1877 ! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
     1878! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
    18781879!   --later used to integrate the hydrostatic equation--
    18791880!==============================================================================
     
    19051906        ! compute sigma level of layer
    19061907        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
    1907        
     1908
    19081909        ! compute "mean" temperature of the layer
    19091910        if(temp(iloop,jloop,kloop,tloop).eq.  &
     
    19161917                     temp(iloop,jloop,kloop-1,tloop))
    19171918        endif
    1918        
     1919
    19191920        ! compute mean value of R of the layer
    19201921        Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
    1921        
     1922
    19221923        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
    19231924        gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
    1924        
     1925
    19251926        ! compute zlocal(kloop)
    19261927        zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
    19271928                log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
    1928        
     1929
    19291930        ! compute za_gcm(kloop)
    19301931        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
     
    19791980
    19801981real,parameter :: P_ref=1 ! reference surface pressure used to build zsurface
    1981 integer i 
     1982integer i
    19821983
    19831984! 1. Build scale heights
     
    20412042
    20422043! Parameters needed to integrate hydrostatic equation:
    2043 real,parameter :: g0=0.62
     2044real,parameter :: g0=0.6198
    20442045!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
    2045 real,parameter :: a0=1187.E3
     2046real,parameter :: a0=1188.E3
    20462047!a0: 'mean' Mars radius=3396.km
    2047 real :: gz 
     2048real :: gz
    20482049! gz: gravitational acceleration at a given (zareoid) altitude
    20492050
     
    20692070        ! compute sigma level of layer
    20702071        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
    2071        
     2072
    20722073        ! compute "mean" temperature of the layer
    20732074        if(temp(iloop,jloop,kloop,tloop).eq.  &
     
    20802081                     temp(iloop,jloop,kloop-1,tloop))
    20812082        endif
    2082                
     2083
    20832084        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
    20842085        ! NB: zareoid=zsurface+phis/g0
    20852086        gz=g0*(a0**2)/ &
    20862087           (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
    2087        
     2088
    20882089        ! compute zs_gcm(iloop,jloop,kloop,tloop)
    20892090        zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
    20902091                log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
    2091        
     2092
    20922093        ! compute za_gcm(kloop)
    20932094!        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
     
    21432144
    21442145! Parameters needed to integrate hydrostatic equation:
    2145 real,parameter :: g0=0.62
     2146real,parameter :: g0=0.6198
    21462147!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
    2147 real,parameter :: a0=1187.E3
     2148real,parameter :: a0=1188.E3
    21482149!a0: 'mean' Mars radius=3396.km
    2149 real :: gz 
     2150real :: gz
    21502151! gz: gravitational acceleration at a given (zareoid) altitude
    21512152
     
    21712172        ! compute sigma level of layer
    21722173        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
    2173        
     2174
    21742175        ! compute "mean" temperature of the layer
    21752176        if(temp(iloop,jloop,kloop,tloop).eq.  &
     
    21822183                     temp(iloop,jloop,kloop-1,tloop))
    21832184        endif
    2184                
     2185
    21852186        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
    21862187        gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
    2187        
     2188
    21882189        ! compute zlocal(kloop)
    21892190        zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
    21902191                log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
    2191        
     2192
    21922193        ! compute za_gcm(kloop)
    21932194        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
     
    22542255        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
    22552256      enddo
    2256      
    2257       ! Interpolation 
     2257
     2258      ! Interpolation
    22582259      do kloop=1,newaltlen
    22592260        ! compute, by interpolation, value at pressure level plevels(kloop)
     
    22682269        endif
    22692270      enddo
    2270      
     2271
    22712272    enddo !iloop
    22722273  enddo !jloop
     
    23392340        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
    23402341      enddo !kloop
    2341      
     2342
    23422343      ! Interpolation
    23432344      do kloop=1,newaltlen
     
    23512352      else ! z_new(kloop) is out of range
    23522353        newdata(iloop,jloop,kloop,tloop)=missing
    2353       endif 
     2354      endif
    23542355      enddo !kloop
    23552356    enddo !iloop
    23562357  enddo !jloop
    2357  enddo !tloop 
     2358 enddo !tloop
    23582359else ! case when flag=1 (i.e.: rho)
    23592360 do tloop=1,tlen
     
    23612362    do iloop=1,lonlen
    23622363      do kloop=1,altlen
    2363         ! extract the vertical coordinates 
     2364        ! extract the vertical coordinates
    23642365        za(kloop)=z_gcm(iloop,jloop,kloop,tloop)
    23652366        ! store log values along altitude
    23662367        logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
    23672368      enddo !kloop
    2368      
     2369
    23692370      ! Interpolation
    23702371      do kloop=1,newaltlen
     
    23822383    enddo !iloop
    23832384  enddo !jloop
    2384  enddo !tloop 
     2385 enddo !tloop
    23852386endif
    23862387
     
    24022403! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
    24032404! "Altitudes" must be above local surface altitudes.
    2404 ! Notes: 
     2405! Notes:
    24052406!       'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
    24062407!       If altitudes in 'znew(i)' fall below those in 'z_gcm(:,:,1,:)', then
     
    24082409!       flag=0, and extrapolated (exponentially) if flag=1.
    24092410!       If altitudes in 'znew(i)' are above those in 'z_gcm(:,:,altlen,:)', then
    2410 !       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)' 
     2411!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)'
    24112412!       if flag=0, and extrapolated (exponentially) if flag=1.
    24122413!==============================================================================
     
    24522453real :: x,y ! temporary variables
    24532454integer :: iloop,jloop,kloop,tloop
    2454 real,parameter :: g0=0.62
     2455real,parameter :: g0=0.6198
    24552456!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
    2456 real,parameter :: a0=1187.E3
     2457real,parameter :: a0=1188.E3
    24572458!a0: 'mean' Mars radius=3396.km
    24582459real,parameter :: Rmean=296.9 ! molecular gas constant
    2459 real :: gz 
     2460real :: gz
    24602461! gz: gravitational acceleration at a given (above areoid) altitude
    24612462
     
    24992500        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
    25002501      enddo !kloop
    2501      
     2502
    25022503      ! Interpolation
    25032504      do kloop=1,newaltlen
     
    25162517    enddo !iloop
    25172518  enddo !jloop
    2518  enddo !tloop 
     2519 enddo !tloop
    25192520else ! case when flag=1 (i.e.: rho or pressure)
    25202521 do tloop=1,tlen
     
    25232524      ! preliminary stuff
    25242525      do kloop=1,altlen
    2525         ! extract the vertical coordinates 
     2526        ! extract the vertical coordinates
    25262527        z(kloop)=z_gcm(iloop,jloop,kloop,tloop)
    25272528        ! store log values along altitude
    25282529        logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
    25292530      enddo !kloop
    2530      
     2531
    25312532      ! Interpolation
    25322533      do kloop=1,newaltlen
     
    25532554    enddo !iloop
    25542555  enddo !jloop
    2555  enddo !tloop 
     2556 enddo !tloop
    25562557endif
    25572558
     
    27142715        double precision x,xi,sum
    27152716
    2716           ae= 1187000.d0
     2717          ae= 1188000.d0
    27172718          gm =869.6 !stern et al 2015 gravitational parameter mu (km3/s-2) = G*M
    27182719          omega=1.138594423d-5    ! Pluto = 2*pi / P   rad/s
     
    92969297        enddo
    92979298        sum = sum+1.
    9298 ! centrifugal potential 
     9299! centrifugal potential
    92999300        v0=(gm*sum/r + 0.5*omega**2 * r**2)
    93009301!       write(out,*)'v0=', v0
     
    93179318        double precision dlon,dlat
    93189319        double precision rg
    9319        
     9320
    93209321        double precision pi,d2r
    93219322        parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
     
    93589359          enddo
    93599360          sum=sum+1.
    9360 ! centrifugal potential 
     9361! centrifugal potential
    93619362          rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
    93629363!          write(out,*) i,r,rg
Note: See TracChangeset for help on using the changeset viewer.