Changeset 3438

Timestamp:

Sep 25, 2024, 4:43:31 PM (2 months ago)

Author:

afalco

Message:

Pluto PCM: Minor fixes.
AF

Location:

trunk/LMDZ.PLUTO

Files:

• libf/dynphy_lonlat/phypluto/newstart.F (modified) (1 diff)
• libf/phypluto/condense_n2.F90 (modified) (3 diffs)
• util/zrecast_plut.F90 (modified) (56 diffs)

trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

@@ -1454 +1454 @@
      &              (lonfi(ig)*180./pi.ge.val4)  .and.
      &              (lonfi(ig)*180./pi.lt.val5)  .and. 
-     &              (qsurf(ig,igcm_n2).lt.val7))  then     
+     &              (qsurf(ig,igcm_n2).gt.val7))  then     
                     qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)*val3
                     qsurf(ig,igcm_co_ice)=qsurf(ig,igcm_co_ice)+val6

trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90

@@ -171 +171 @@
 
      ! calculate global mean surface pressure for the fast mode
+     IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
+     DO ig=1,klon
+        kp(ig) = exp(-phisfi(ig)/(r*38.))
+     ENDDO
      IF (fast) THEN
-        IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
-        DO ig=1,klon
-           kp(ig) = exp(-phisfi(ig)/(r*38.))
-        ENDDO
         p00=glob_average2d(kp) ! mean pres at ref level
      ENDIF
@@ -382 +382 @@
         pdpsrf(ig)   = -pdicen2(ig)*g
     !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
-        IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
+        IF (fast.and.(zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001)) then
             pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
             pdicen2(ig)=-pdpsrf(ig)/g
@@ -855 +855 @@
          Mtot = masse(m)
          MQtot = masse(m)*q(m)
-         !if (m.lt.klev) then ! because some compilers will have problems
-         !                      ! evaluating masse(klev+1)
+         if (m.lt.klev) then ! because some compilers will have problems
+                              ! evaluating masse(klev+1)
          do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
             m=m+1
             Mtot = Mtot + masse(m)
             MQtot = MQtot + masse(m)*q(m)
-         !   if (m.eq.klev) exit
+           if (m.eq.klev) exit
          end do
-         !endif
+         endif
          if (m.lt.klev) then
             sigw=(w(l+1)-Mtot)/masse(m+1)

trunk/LMDZ.PLUTO/util/zrecast_plut.F90

@@ -40 +40 @@
 !              'phisinit.nc')
 ! EM 02/2007 : Changed behavior for "altitude above surface" case
-!              (for MCD RMS computation). Number of levels is then set as 
+!              (for MCD RMS computation). Number of levels is then set as
 !              number of levels in initial file,
 !              and the new set of above surface levels follow a more elaborate
@@ -99 +99 @@
 real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
 real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
-real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure 
+real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
 real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
 real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
@@ -195 +195 @@
     if (tmpndims.eq.4) then
       nbvar=nbvar+1
-      var(nbvar)=tmpvarname    
+      var(nbvar)=tmpvarname
     else
       write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
@@ -490 +490 @@
 endif
 ! look for 'phisinit' in current file
-ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid) 
+ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
 if (ierr.ne.NF_NOERR) then
   write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
@@ -598 +598 @@
     if (nblev.eq.1) then ! in case only one level is asked for
       zareoid(nblev)=zamin
-    else 
+    else
       do i=1,nblev
         zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
@@ -625 +625 @@
     if (nblev.eq.1) then ! in case only one level is asked for
       zradius(nblev)=zamin
-    else 
+    else
       do i=1,nblev
         zradius(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
@@ -977 +977 @@
 ! 3.3.1 Define 1D variables
 
-! longitude 
+! longitude
 datashape(1)=lon_dimid
 ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
@@ -1156 +1156 @@
   stop "Error: Problem writing long_name for Time"
 endif
-text='days since 0000-01-1 00:00:00'
+!text='days since 0000-01-1 00:00:00'
+text='days since 00000'
 ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
 if (ierr.ne.NF_NOERR) then
@@ -1259 +1260 @@
     write(*,*) 'Error, could not define variable ',trim(var(i))
     write(*,*) NF_STRERROR(ierr)
-    stop 
+    stop
   endif
 
@@ -1320 +1321 @@
     endif
   endif
-  
+
   ! look for a "units" attribute
   text='   '
@@ -1334 +1335 @@
     endif
   endif
-  
+
   ! look for a "missing_value" attribute
   ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
@@ -1343 +1344 @@
     miss_val=miss_val_def
   endif
-  
+
   ! write the missing_value attribute to output
   ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
@@ -1356 +1357 @@
     write(*,*)'nbattr:',nbattr,' j:',j
   endif
-    
+
 enddo ! of do=1,nbvar
 
@@ -1499 +1500 @@
     stop
   endif
-  
+
   ! interpolate dataset onto new grid
   call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                       miss_val,ps,press,indata,plevel,outdata)
-  
+
   ! write new dataset to output file
   ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
@@ -1511 +1512 @@
     stop
   endif
-  
+
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate zareoid onto new grid
   call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
@@ -1540 +1541 @@
     stop
   endif
-  
+
   ! interpolate dataset onto new grid
   ! check if variable is "rho" (to set flag for interpolation below)
@@ -1548 +1549 @@
     j=0
   endif
-  
+
   call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                       miss_val,za_gcm,indata,j,zareoid,outdata)
@@ -1561 +1562 @@
 
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate pressure onto new grid
   write(*,*) ""
@@ -1593 +1594 @@
       stop
     endif
-  
+
     ! interpolate dataset onto new grid
     ! check if variable is "rho" (to set flag for interpolation below)
@@ -1601 +1602 @@
       j=0
     endif
-  
+
     call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                         miss_val,zs_gcm,indata,phisinit,press,temp,rho,&
@@ -1614 +1615 @@
     endif
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate pressure onto new grid
   write(*,*) ""
@@ -1646 +1647 @@
     stop
   endif
-  
+
   ! interpolate dataset onto new grid
   ! check if variable is "rho" (to set flag for interpolation below)
@@ -1654 +1655 @@
     j=0
   endif
-  
+
   call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                       miss_val,ra_gcm,indata,j,zradius,outdata)
@@ -1667 +1668 @@
 
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate pressure onto new grid
   write(*,*) ""
@@ -1734 +1735 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -1749 +1750 @@
 
 !==============================================================================
-! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
+! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
 !   --later used to integrate the hydrostatic equation--
 !==============================================================================
@@ -1779 +1780 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -1790 +1791 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-        
+
         ! compute mean value of R of the layer
         Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
-        
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         ! NB: zareoid=zsurface+phis/g0
         gz=g0*(a0**2)/ &
            (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
-        
+
         ! compute zs_gcm(iloop,jloop,kloop,tloop)
         zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
 !        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -1860 +1861 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' radius
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -1875 +1876 @@
 
 !==============================================================================
-! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
+! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
 !   --later used to integrate the hydrostatic equation--
 !==============================================================================
@@ -1905 +1906 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -1916 +1917 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-        
+
         ! compute mean value of R of the layer
         Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
-        
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
-        
+
         ! compute zlocal(kloop)
         zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
         za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -1979 +1980 @@
 
 real,parameter :: P_ref=1 ! reference surface pressure used to build zsurface
-integer i 
+integer i
 
 ! 1. Build scale heights
@@ -2041 +2042 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -2069 +2070 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -2080 +2081 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-                
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         ! NB: zareoid=zsurface+phis/g0
         gz=g0*(a0**2)/ &
            (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
-        
+
         ! compute zs_gcm(iloop,jloop,kloop,tloop)
         zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
 !        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -2143 +2144 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -2171 +2172 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -2182 +2183 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-                
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
-        
+
         ! compute zlocal(kloop)
         zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
         za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -2254 +2255 @@
         q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
       enddo
-      
-      ! Interpolation  
+
+      ! Interpolation
       do kloop=1,newaltlen
         ! compute, by interpolation, value at pressure level plevels(kloop)
@@ -2268 +2269 @@
         endif
       enddo
-      
+
     enddo !iloop
   enddo !jloop
@@ -2339 +2340 @@
         q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2351 +2352 @@
       else ! z_new(kloop) is out of range
         newdata(iloop,jloop,kloop,tloop)=missing
-      endif 
+      endif
       enddo !kloop
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 else ! case when flag=1 (i.e.: rho)
  do tloop=1,tlen
@@ -2361 +2362 @@
     do iloop=1,lonlen
       do kloop=1,altlen
-        ! extract the vertical coordinates 
+        ! extract the vertical coordinates
         za(kloop)=z_gcm(iloop,jloop,kloop,tloop)
         ! store log values along altitude
         logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2382 +2383 @@
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 endif
 
@@ -2402 +2403 @@
 ! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
 ! "Altitudes" must be above local surface altitudes.
-! Notes: 
+! Notes:
 !       'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
 !       If altitudes in 'znew(i)' fall below those in 'z_gcm(:,:,1,:)', then
@@ -2408 +2409 @@
 !       flag=0, and extrapolated (exponentially) if flag=1.
 !       If altitudes in 'znew(i)' are above those in 'z_gcm(:,:,altlen,:)', then
-!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)' 
+!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)'
 !       if flag=0, and extrapolated (exponentially) if flag=1.
 !==============================================================================
@@ -2452 +2453 @@
 real :: x,y ! temporary variables
 integer :: iloop,jloop,kloop,tloop
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
 real,parameter :: Rmean=296.9 ! molecular gas constant
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (above areoid) altitude
 
@@ -2499 +2500 @@
         q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2516 +2517 @@
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 else ! case when flag=1 (i.e.: rho or pressure)
  do tloop=1,tlen
@@ -2523 +2524 @@
       ! preliminary stuff
       do kloop=1,altlen
-        ! extract the vertical coordinates 
+        ! extract the vertical coordinates
         z(kloop)=z_gcm(iloop,jloop,kloop,tloop)
         ! store log values along altitude
         logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2553 +2554 @@
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 endif
 
@@ -2714 +2715 @@
         double precision x,xi,sum
 
-          ae= 1187000.d0
+          ae= 1188000.d0
           gm =869.6 !stern et al 2015 gravitational parameter mu (km3/s-2) = G*M
           omega=1.138594423d-5    ! Pluto = 2*pi / P   rad/s
@@ -9296 +9297 @@
         enddo
         sum = sum+1.
-! centrifugal potential 
+! centrifugal potential
         v0=(gm*sum/r + 0.5*omega**2 * r**2)
 !       write(out,*)'v0=', v0
@@ -9317 +9318 @@
         double precision dlon,dlat
         double precision rg
-        
+
         double precision pi,d2r
         parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
@@ -9358 +9359 @@
           enddo
           sum=sum+1.
-! centrifugal potential 
+! centrifugal potential
           rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
 !          write(out,*) i,r,rg