Changeset 3275 for trunk/LMDZ.PLUTO.old

Timestamp:

Mar 20, 2024, 3:05:14 PM (8 months ago)

Author:

afalco

Message:

Pluto PCM:
Changed _vap to _gas;
Included surfprop.F90;
callcorrk includes methane
AF

Location:

trunk/LMDZ.PLUTO.old

Files:

• . (modified) (1 prop)
• compile1d (modified) (1 diff)
• deftank (modified) (1 prop)
• deftank/gcm (modified) (1 prop)
• deftank/kbo_def (modified) (1 prop)
• deftank/nogcm (modified) (1 prop)
• deftank/nogcm_simple (modified) (1 prop)
• deftank/testphys1d (modified) (1 prop)
• libf (modified) (1 prop)
• libf/bibio (modified) (1 prop)
• libf/dyn3d (modified) (1 prop)
• libf/dyn3d/poubelle (modified) (1 prop)
• libf/dyn3d/startsHD (modified) (1 prop)
• libf/dyn3d/stock (modified) (1 prop)
• libf/filtrez (modified) (1 prop)
• libf/grid (modified) (1 prop)
• libf/grid/dimension (modified) (1 prop)
• libf/phypluto (modified) (1 prop)
• libf/phypluto/callcorrk.F (modified) (39 diffs)
• libf/phypluto/cooling_stock (modified) (1 prop)
• libf/phypluto/physiq.F (modified) (2 diffs)
• libo (modified) (1 prop)

trunk/LMDZ.PLUTO.old/compile1d

@@ -1 +1 @@
 #!/bin/bash
-
-./makegcm_cygwin -debug -fdefault-real-8 -d 25 -b 17x23 -t 7 -s 1 -p pluto testphys1d
+export PATH=$PATH:.
+./makegcm_spirit_gfortran -debug -fdefault-real-8 -d 25 -t 7 -s 2 -b 17x23 -p pluto testphys1d
+# ./makegcm_cygwin -debug -fdefault-real-8 -d 25 -b 17x23 -t 7 -s 1 -p pluto testphys1d

trunk/LMDZ.PLUTO.old/libf/phypluto/callcorrk.F

@@ -5 +5 @@
      &     fluxsurf_sw,fluxtop_lw,fluxtop_sw,fluxtop_dn,
      &     reffrad,tau_col,ptime,pday,firstcall,lastcall,zzlay)
-         
+
       use radinc_h
       use radcommon_h
-      use ioipsl_getincom 
+      use ioipsl_getincom
       use radii_mod
       use aerosol_mod
@@ -24 +24 @@
 !
 !     Authors
-!     ------- 
+!     -------
 !     Emmanuel 01/2001, Forget 09/2001
 !     Robin Wordsworth (2009)
@@ -37 +37 @@
 !-----------------------------------------------------------------------
 !     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
-!     Layer #1 is the layer near the ground. 
+!     Layer #1 is the layer near the ground.
 !     Layer #nlayermx is the layer at the top.
 
@@ -44 +44 @@
       INTEGER ngrid,nlayer
       INTEGER igout
-      REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau 
+      REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau
       REAL albedo(ngrid)                    ! SW albedo
       REAL emis(ngrid)                      ! LW emissivity
@@ -110 +110 @@
       real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
       REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
-      real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) ! 
-      real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) ! 
+      real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) !
+      real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) !
       REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
       REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
@@ -137 +137 @@
       save qxvaer, qsvaer, gvaer
       save qxiaer, qsiaer, giaer
-      save QREFvis3d, QREFir3d 
+      save QREFvis3d, QREFir3d
 
       REAL tau_col(ngrid) ! diagnostic from aeropacity
@@ -220 +220 @@
          call setspv            ! basic visible properties
 
-         ! Radiative Hazes 
+         ! Radiative Hazes
          if (aerohaze) then
 
@@ -231 +231 @@
            !--------------------------------------------------
            do iaer=1,naerkind
-              if ((iaer.eq.iaero_haze)) then 
-               call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer), 
+              if ((iaer.eq.iaero_haze)) then
+               call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer),
      &             nueffrad(1,1,iaer))
               endif
@@ -238 +238 @@
            if (haze_radproffix) then
               print*, 'haze_radproffix=T : fixed profile for haze rad'
-           else 
+           else
               print*,'reffrad haze:',reffrad(1,1,iaero_haze)
               print*,'nueff haze',nueffrad(1,1,iaero_haze)
@@ -272 +272 @@
 !-----------------------------------------------------------------------
 !     Get 3D aerosol optical properties.
-      ! ici on selectionne les proprietes opt correspondant a reffrad 
+      ! ici on selectionne les proprietes opt correspondant a reffrad
       if (aerohaze) then
         !--------------------------------------------------
@@ -283 +283 @@
         endif
 
-        call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         
-     &       QVISsQREF3d,omegaVIS3d,gVIS3d,                         
-     &      QIRsQREF3d,omegaIR3d,gIR3d,                             
+        call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
+     &       QVISsQREF3d,omegaVIS3d,gVIS3d,
+     &      QIRsQREF3d,omegaIR3d,gIR3d,
      &      QREFvis3d,QREFir3d)
 
         ! Get aerosol optical depths.
-        call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      
-     &      reffrad,QREFvis3d,QREFir3d,                             
+        call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,
+     &      reffrad,QREFvis3d,QREFir3d,
      &      tau_col)
       endif
@@ -298 +298 @@
       IF (methane) then
         vmrch4(:,:)=0.
-        
+
         if (ch4fix) then
            if (vmrch4_proffix) then
             !! Interpolate on the model vertical grid
              do ig=1,ngridmx
-               CALL interp_line(levdat,vmrdat,Nfine,
-     &                  zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
+    !            CALL interp_line(levdat,vmrdat,Nfine,
+    !  &                  zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
              enddo
            else
@@ -317 +317 @@
 !     Prepare NON LTE correction in Pluto atmosphere
       IF (nlte) then
-        CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
-     &             eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
+    !     CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,
+    !  &             eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw)
       ENDIF
 c     Net atmospheric radiative cooling rate from C2H2 (K.s-1):
@@ -343 +343 @@
 !     shortwave
             do iaer=1,naerkind
-               DO nw=1,L_NSPECTV 
+               DO nw=1,L_NSPECTV
                   do l=1,nlayermx
 
-                     temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) 
+                     temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer)
      $                    *QREFvis3d(ig,nlayermx+1-l,iaer)
 
-                     temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) 
+                     temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
      $                    *QREFvis3d(ig,max(nlayermx-l,1),iaer)
                      qxvaer(2*l,nw,iaer)  = temp1
@@ -378 +378 @@
 
 !     longwave
-               DO nw=1,L_NSPECTI 
+               DO nw=1,L_NSPECTI
                   do l=1,nlayermx
 
-                     temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) 
+                     temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer)
      $                    *QREFir3d(ig,nlayermx+1-l,iaer)
 
-                     temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) 
+                     temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer)
      $                    *QREFir3d(ig,max(nlayermx-l,1),iaer)
 
@@ -421 +421 @@
                   do nw=1,L_NSPECTV
                      if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
-                        print*,'Serious problems with qsvaer values' 
+                        print*,'Serious problems with qsvaer values'
                         print*,'in callcorrk'
                         call abort
@@ -430 +430 @@
                   end do
 
-                  do nw=1,L_NSPECTI 
+                  do nw=1,L_NSPECTI
                      if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
                         print*,'Serious problems with qsiaer values'
@@ -448 +448 @@
 !-----------------------------------------------------------------------
 !     Aerosol optical depths
-         IF (aerohaze) THEN   
-          do iaer=1,naerkind     ! heritage generic        
+         IF (aerohaze) THEN
+          do iaer=1,naerkind     ! heritage generic
             do k=0,nlayer-1
                pweight=
      $              (pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/
      $              (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
-               if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then 
+               if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then
                  temp=aerosol(ig,L_NLAYRAD-k,iaer)/
      $              QREFvis3d(ig,L_NLAYRAD-k,iaer)
@@ -479 +479 @@
 !     Albedo and emissivity
          albi=1-emis(ig)        ! longwave
-         albv=albedo(ig)        ! shortwave 
+         albv=albedo(ig)        ! shortwave
          acosz=mu0(ig)          ! cosine of sun incident angle
 
 !-----------------------------------------------------------------------
-!     Methane vapour 
+!     Methane vapour
 
 c     qvar = mixing ratio
@@ -490 +490 @@
 c     datagcm/composition.in for the k-coefficients.
          qvar(:)=0.
-         IF (methane) then 
+         IF (methane) then
 
            do l=1,nlayer
@@ -554 +554 @@
 !! following lines changed in 03/2015 to solve upper atmosphere bug
 !        plevrad(1) = 0.
-!        plevrad(2) = max(pgasmin,0.0001*plevrad(3))      
+!        plevrad(2) = max(pgasmin,0.0001*plevrad(3))
 !
 !        tlevrad(1) = tlevrad(2)
@@ -563 +563 @@
 !
 !        pmid(1) = plevrad(2)
-!        pmid(2) = plevrad(2) 
+!        pmid(2) = plevrad(2)
 
          DO l=1,L_NLAYRAD-1
@@ -574 +574 @@
          pmid(L_LEVELS) = plevrad(L_LEVELS)
          tmid(L_LEVELS) = tlevrad(L_LEVELS)
-      
+
       !TB
          if ((PMID(2).le.1.e-5).and.(ig.eq.1)) then
@@ -608 +608 @@
           endif
          enddo
-     
+
 !=======================================================================
 !     Calling the main radiative transfer subroutines
@@ -614 +614 @@
 !-----------------------------------------------------------------------
 !     Shortwave
-        
+
          IF(fract(ig) .GE. 1.0e-4) THEN ! only during daylight  IPM?! flux UV...
 
@@ -623 +623 @@
             END DO
 
-            !print*, 'starting optcv' 
+            !print*, 'starting optcv'
             call optcv(dtauv,tauv,taucumv,plevrad,
      $           qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,
@@ -666 +666 @@
 !        IR spectral output from top of the atmosphere
          if(specOLR)then
-            do nw=1,L_NSPECTI 
+            do nw=1,L_NSPECTI
                OLR_nu(ig,nw)=nfluxtopi_nu(nw)
             end do
@@ -673 +673 @@
 ! **********************************************************
 !     Finally, the heating rates
-!     g/cp*DF/DP 
+!     g/cp*DF/DP
 ! **********************************************************
 
@@ -682 +682 @@
                !dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))*dpp   !averaged dtlw on each wavelength
                do nw=1,L_NSPECTV
-                 dtsw_nu(L_NLAYRAD+1-l,nw)= 
+                 dtsw_nu(L_NLAYRAD+1-l,nw)=
      &              (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp
                end do
@@ -689 +689 @@
                !dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))*dpp   !averaged dtlw on each wavelength
                do nw=1,L_NSPECTI
-                 dtlw_nu(L_NLAYRAD+1-l,nw)= 
-     &              (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp      
+                 dtlw_nu(L_NLAYRAD+1-l,nw)=
+     &              (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp
                end do
-         END DO      
-         
+         END DO
+
          ! values at top of atmosphere
          dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1)))
 
-         ! SW 
+         ! SW
          !dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)*dpp
          do nw=1,L_NSPECTV
@@ -704 +704 @@
          end do
 
-         ! LW 
-c        dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp 
+         ! LW
+c        dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp
          do nw=1,L_NSPECTI
              dtlw_nu(L_NLAYRAD,nw)=
@@ -717 +717 @@
 
          if (.not.nlte) then
-            eps_nlte_sw23(ig,:) =1. ! IF no NLTE  
-            eps_nlte_sw33(ig,:) =1. ! IF no NLTE  
-            eps_nlte_lw(ig,:) =1. ! IF no NLTE 
+            eps_nlte_sw23(ig,:) =1. ! IF no NLTE
+            eps_nlte_sw33(ig,:) =1. ! IF no NLTE
+            eps_nlte_lw(ig,:) =1. ! IF no NLTE
          endif
- 
+
          do l=1,nlayer
 
             !LW
-            dtlw(ig,l) =0. 
+            dtlw(ig,l) =0.
 !           dtlw_co(ig,l) =0.  ! only for diagnostic
             do nw=1,L_NSPECTI
               ! wewei : wavelength in micrometer
-              if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then 
+              if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then
                 dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l)
-              else 
+              else
                 !dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996)
                 dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996)
 !               dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic
               end if
-              dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength 
+              dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength
             end do
             ! adding c2h2 if cooling active
@@ -743 +743 @@
             !SW
             dtsw(ig,l) =0.
- 
+
             if (strobel) then
 
              do nw=1,L_NSPECTV
-              if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then 
+              if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then
                 dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l)
-              elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then 
+              elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then
                 dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l)
               else
                 dtsw_nu(l,nw)=dtsw_nu(l,nw)
               end if
-              dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)  
+              dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw)
              end do
 
@@ -764 +764 @@
              enddo
 
-            endif 
+            endif
 
 
@@ -771 +771 @@
 
 !     Diagnotics for last call for each grid point
-         !if (lastcall) then 
+         !if (lastcall) then
 
           !print*,'albedi vis=',albv
@@ -806 +806 @@
       endif
 
-      if(lastcall)then 
+      if(lastcall)then
 
         ! 1D Output
@@ -816 +816 @@
                open(116,file='surf_vals.out')
                write(116,*) tsurf(1),pplev(1,1),
-     &             fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1) 
+     &             fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1)
                do nw=1,L_NSPECTV
                  write(116,*) wavev(nw),fmnetv_nu(L_NLAYRAD,nw)
@@ -830 +830 @@
            if(diagrad_OLR)then
                open(117,file='OLRnu.out')
-               write(117,*) 'IR wavel - band width - OLR' 
+               write(117,*) 'IR wavel - band width - OLR'
                do nw=1,L_NSPECTI
                    write(117,*) wavei(nw),
      &        abs(1.e4/bwnv(nw)-1.e4/bwnv(nw+1)),OLR_nu(1,nw)
-               enddo   
+               enddo
                close(117)
            endif
@@ -846 +846 @@
                  write(118,*) plevrad(2*l)
                  do nw=1,L_NSPECTI
-                     write(119,*) fluxupi_nu(l,nw) 
+                     write(119,*) fluxupi_nu(l,nw)
                  enddo
-               enddo 
+               enddo
                close(118)
                close(119)

trunk/LMDZ.PLUTO.old/libf/phypluto/physiq.F

@@ -1849 +1849 @@
           call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
           call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
-         ! call WRITEDIAGFI(ngrid,"pressure","Pression","Pa",3,pplay)
+          call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)
           call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
      &                                        "W m-2",2,fluxrad)
@@ -2073 +2073 @@
           call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",0,ps)
           call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+          call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)
 
           call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",