Changeset 2935 for trunk/LMDZ.MARS/libf

Timestamp:

Apr 12, 2023, 1:44:55 PM (20 months ago)

Author:

romain.vande

Message:

Mars newstart : Adapt routine lect_start_archive.F to the subslope parametrisation

File:

• trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F (modified) (30 diffs)

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F

@@ -23 +23 @@
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff
       USE comconst_mod, ONLY: kappa,g,pi
+      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
+     &             subslope_dist,end_comslope_h,ini_comslope_h
+      use netcdf
       implicit none
 
@@ -44 +47 @@
       REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1)
       REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot)
-      REAL,INTENT(OUT) :: tsurf(ngrid) ! surface temperature
-      REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx) ! soil temperature
-      REAL,INTENT(OUT) :: albedo(ngrid,2) ! surface albedo
-      REAL,INTENT(OUT) :: emis(ngrid) ! ground emissivity
-      REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
+      REAL,INTENT(OUT) :: tsurf(ngrid,nslope) ! surface temperature
+      REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx,nslope) ! soil temperature
+      REAL,INTENT(OUT) :: albedo(ngrid,2,nslope) ! surface albedo
+      REAL,INTENT(OUT) :: emis(ngrid,nslope) ! ground emissivity
+      REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot,nslope)
       REAL,INTENT(OUT) :: tauscaling(ngrid) ! dust conversion factor
       REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! sub grid cloud fraction
-      REAL,INTENT(OUT) :: watercap(ngrid) ! infinite polar cap
+      REAL,INTENT(OUT) :: watercap(ngrid,nslope) ! infinite polar cap
       REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1)
       REAL,INTENT(OUT) :: t(iip1,jjp1,llm)
@@ -128 +131 @@
 c------------------------------------------------------
       real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
-      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
+      real tsurfS(iip1,jjp1,nslope),tsoilS(iip1,jjp1,nsoilmx,nslope)
       real inertiedatS(iip1,jjp1,nsoilmx)
-      real co2iceS(iip1,jjp1),emisS(iip1,jjp1)
-      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
+      real co2iceS(iip1,jjp1,nslope),emisS(iip1,jjp1,nslope)
+      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot,nslope)
       real tauscalingS(iip1,jjp1)
       real totcloudfracS(iip1,jjp1)
-      real watercapS(iip1,jjp1)
-      real albedoS(iip1,jjp1)
+      real watercapS(iip1,jjp1,nslope)
+      real albedoS(iip1,jjp1,nslope)
 
       real ptotal, co2icetotal
@@ -153 +156 @@
       real, dimension(:), allocatable :: mlayerold
       real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
-      real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
-      real, dimension(:,:,:),allocatable :: tsoiloldnew
+      real, dimension(:,:,:,:), allocatable :: tsoilold,qsurfold
+      real, dimension(:,:,:),allocatable :: tsoiloldnew_noslope
+      real, dimension(:,:,:), allocatable :: tsoilold_noslope
+      real, dimension(:,:,:), allocatable ::qsurfold_noslope
+      real, dimension(:,:,:,:),allocatable :: tsoiloldnew
 ! tsoiloldnew: old soil values, but along new subterranean grid
       real, dimension(:,:,:), allocatable :: inertiedatold
@@ -160 +166 @@
       real, dimension(:,:,:), allocatable :: inertiedatoldnew
       real, dimension(:,:), allocatable :: psold,phisold
-      real, dimension(:,:), allocatable :: co2iceold,tsurfold
-      real, dimension(:,:), allocatable :: emisold
+      real, dimension(:,:,:), allocatable :: co2iceold,tsurfold
+      real, dimension(:,:), allocatable :: co2iceold_noslope
+      real, dimension(:,:), allocatable :: tsurfold_noslope
+      real, dimension(:,:,:), allocatable :: emisold
+      real, dimension(:,:), allocatable :: emisold_noslope
       real, dimension(:,:,:,:), allocatable :: qold
       real, dimension(:,:), allocatable :: tauscalingold
       real, dimension(:,:), allocatable :: totcloudfracold
-      real, dimension(:,:), allocatable :: watercapold
-      real, dimension(:,:), allocatable :: albedoold
+      real, dimension(:,:,:), allocatable :: watercapold
+      real, dimension(:,:), allocatable :: watercapold_noslope
+      real, dimension(:,:,:), allocatable :: albedoold
+      real, dimension(:,:), allocatable :: albedoold_noslope
+
 
       real tab_cntrl(100)
@@ -200 +212 @@
       ! start_archive) and then .false.
       logical :: old_co2ice=.false. 
+      
+      ! flag to check if nslope is not present (old start_archive)
+      integer :: nslope_read
+      logical :: no_slope=.false. 
+      integer :: ndims
+      integer, dimension(:), allocatable :: dimids
 c=======================================================================
 
@@ -324 +342 @@
         depthinterpol=.true.
       endif
+      
+      
+! 1.2.1 find out the # of subgrid scale slope
+      ierr= NF_INQ_DIMID(nid,"nslope",dimid)
+      if (ierr.ne.NF_NOERR) then
+         write(*,*) "nslope not present, old startarchive"
+         write(*,*) "nslope=1"
+         nslope=1
+         no_slope=.true.
+      else
+         write(*,*) "nslope present"
+         ierr= NF_INQ_DIMLEN(nid,dimid,nslope_read)
+         if (ierr.ne.NF_NOERR) then
+          write(*,*) 'lect_start_archive error: cannot 
+     & find nslope length'
+          stop
+         else
+          write(*,*) "lect_start_archive: nslope_read=",nslope_read
+          if(nslope_read.ne.1) then
+            write(*,*) 'lect_start_archive only works with nslope=1'
+            stop
+          else
+            nslope=1
+            no_slope=.false.
+          endif
+         endif
+      endif
+      def_slope(1)=-90.
+      def_slope(2)=90.
+      def_slope_mean(1)=0.
+      subslope_dist(:,:)=1.
 
 ! 1.3 Report dimensions
@@ -331 +380 @@
       write(*,*) "latitude: ",jmold
       write(*,*) "altitude: ",lmold
+      write(*,*) "nslope: ",nslope
       if (nqold.gt.0) then
         write(*,*) "old tracers q*: ",nqold
@@ -358 +408 @@
       allocate(phisold(imold+1,jmold+1))
       allocate(qold(imold+1,jmold+1,lmold,nqtot))
-      allocate(co2iceold(imold+1,jmold+1))
-      allocate(tsurfold(imold+1,jmold+1))
-      allocate(emisold(imold+1,jmold+1))
+      allocate(co2iceold(imold+1,jmold+1,nslope))
+      allocate(tsurfold(imold+1,jmold+1,nslope))
+      allocate(emisold(imold+1,jmold+1,nslope))
       allocate(q2old(imold+1,jmold+1,lmold+1))
 !      allocate(tsoilold(imold+1,jmold+1,nsoilmx))
-      allocate(tsoilold(imold+1,jmold+1,nsoilold))
-      allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx))
+      allocate(tsoilold(imold+1,jmold+1,nsoilold,nslope))
+      allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx,nslope))
       allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia
       allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx))
@@ -370 +420 @@
       allocate(surfithold(imold+1,jmold+1))
       allocate(mlayerold(nsoilold))
-      allocate(qsurfold(imold+1,jmold+1,nqtot))
+      allocate(qsurfold(imold+1,jmold+1,nqtot,nslope))
       allocate(tauscalingold(imold+1,jmold+1))
       allocate(totcloudfracold(imold+1,jmold+1))
-      allocate(watercapold(imold+1,jmold+1))
-      allocate(albedoold(imold+1,jmold+1))
+      allocate(watercapold(imold+1,jmold+1,nslope))
+      allocate(albedoold(imold+1,jmold+1,nslope))
+      
+      if(no_slope) then
+      allocate(co2iceold_noslope(imold+1,jmold+1))
+      allocate(tsurfold_noslope(imold+1,jmold+1))
+      allocate(emisold_noslope(imold+1,jmold+1))
+      allocate(watercapold_noslope(imold+1,jmold+1))
+      allocate(albedoold_noslope(imold+1,jmold+1))
+      allocate(qsurfold_noslope(imold+1,jmold+1,nqtot))
+      allocate(tsoilold_noslope(imold+1,jmold+1,nsoilold))
+      allocate(tsoiloldnew_noslope(imold+1,jmold+1,nsoilmx))
+      endif
 
       allocate(var (imold+1,jmold+1,llm))
@@ -722 +783 @@
         ENDIF
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
-#endif
+      if(no_slope) then
+      ierr = nf90_get_var(nid, nvarid,co2iceold_noslope)
+      co2iceold(:,:,1)=co2iceold_noslope(:,:)
+      else
+      ierr = nf90_get_var(nid, nvarid,co2iceold)
+      endif
       IF (ierr .NE. NF_NOERR) THEN
         PRINT*, "lect_start_archive: Failed loading <co2ice>"
@@ -738 +800 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
-#endif
+      if(no_slope) then
+      ierr = nf90_get_var(nid, nvarid,emisold_noslope)
+      emisold(:,:,1)=emisold_noslope(:,:)
+      else
+       ierr = nf90_get_var(nid, nvarid,emisold)
+      endif
       IF (ierr .NE. NF_NOERR) THEN
          PRINT*, "lect_start_archive: Failed loading <emis>"
@@ -768 +831 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
-#endif
+      if(no_slope) then
+        ierr = nf90_get_var(nid, nvarid,tsurfold_noslope)
+        tsurfold(:,:,1)=tsurfold_noslope(:,:)
+      else
+        ierr = nf90_get_var(nid, nvarid,tsurfold)
+      endif
       IF (ierr .NE. NF_NOERR) THEN
          PRINT*, "lect_start_archive: Failed loading <tsurf>"
@@ -832 +896 @@
       IF (ierr .NE. NF_NOERR) THEN
          PRINT*, "lect_start_archive: <watercap> not in file"
-         watercapold(:,:) = 0.
+         watercapold(:,:,:) = 0.
       ELSE
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
-     &                            watercapold)
-#else
-        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
-     &                          watercapold)
-#endif
+      if(no_slope) then
+       ierr = nf90_get_var(nid, nvarid,watercapold_noslope)
+       watercapold(:,:,1)=watercapold_noslope(:,:)
+      else
+       ierr = nf90_get_var(nid, nvarid,watercapold)
+      endif
         IF (ierr .NE. NF_NOERR) THEN
            PRINT*, "lect_start_archive: Failed loading <watercap>"
@@ -851 +914 @@
       IF (ierr .NE. NF_NOERR) THEN
          PRINT*, "lect_start_archive: <albedo> not in file"
-         albedoold(:,:) = 0.
+         albedoold(:,:,:) = 0.
       ELSE
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,count,
-     &                            albedoold)
-#else
-        ierr = NF_GET_VARA_REAL(nid,nvarid,start,count,
-     &                          albedoold)
-#endif
+      if(no_slope) then
+       ierr = nf90_get_var(nid, nvarid,albedoold_noslope)
+      albedoold(:,:,1)=albedoold_noslope(:,:)
+      else
+       ierr = nf90_get_var(nid, nvarid,albedoold)
+      endif
         IF (ierr .NE. NF_NOERR) THEN
            PRINT*, "lect_start_archive: Failed loading <albedo>"
@@ -876 +938 @@
 c -------------------------------------------
 ! Surface tracers:      
-      qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0
+      qsurfold(1:imold+1,1:jmold+1,1:nqtot,1:nslope)=0
 
       DO iq=1,nqtot
@@ -906 +968 @@
           print*, "which (constant) value should it be initialized to?"
           read(*,*) tmpval
-          qsurfold(1:imold+1,1:jmold+1,iq)=tmpval
+          qsurfold(1:imold+1,1:jmold+1,iq,1:nslope)=tmpval
         ELSE ! tracer exists in file, load it
-#ifdef NC_DOUBLE
-          ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
-     &          qsurfold(1,1,iq))
-#else
-          ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
-     &          qsurfold(1,1,iq))
-#endif
+        if(no_slope) then
+          ierr = nf90_get_var(nid, nvarid,qsurfold_noslope(:,:,iq))
+          qsurfold(:,:,iq,1)=qsurfold_noslope(:,:,iq)
+         else
+           ierr = nf90_get_var(nid, nvarid,qsurfold(:,:,iq,:))
+         endif
           IF (ierr .NE. NF_NOERR) THEN
             PRINT*, "lect_start_archive: ",
@@ -946 +1007 @@
             CALL abort
          ENDIF
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
-     &          tsoilold(1,1,isoil))
-#else
-         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
-     &          tsoilold(1,1,isoil))
-#endif
+         if(no_slope) then
+          ierr = nf90_get_var(nid, nvarid,tsoilold_noslope(:,:,isoil))
+         tsoilold(:,:,:,1)=tsoilold_noslope(:,:,:)
+         else
+          ierr = nf90_get_var(nid, nvarid,tsoilold(:,:,isoil,:))
+         endif
          IF (ierr .NE. NF_NOERR) THEN
             PRINT*, "lect_start_archive: ",
@@ -972 +1032 @@
         call abort
        else
-#ifdef NC_DOUBLE
-      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
-#else
-      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
-#endif
+       if(no_slope) then
+        ierr = nf90_get_var(nid, nvarid,tsoilold_noslope)
+        tsoilold(:,:,:,1)=tsoilold_noslope(:,:,:)
+      else
+        ierr = nf90_get_var(nid, nvarid,tsoilold)
+       endif
        endif ! of if (ierr.ne.NF_NOERR)
        
@@ -1130 +1191 @@
 
 c Glace CO2
-      call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call interp_horiz (co2iceold(:,:,1),co2ices(:,:,1),
+     & imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv)
 
 c Temperature de surface
-      call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs,tsurf)
+      call interp_horiz (tsurfold(:,:,1),tsurfs(:,:,1)
+     & ,imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs(:,:,1),tsurf(:,1))
 c     write(44,*) 'tsurf', tsurf
 
@@ -1147 +1208 @@
 
 c Emissivite de la surface
-      call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss,emis)
+      call interp_horiz (emisold(:,:,1),emiss(:,:,1)
+     & ,imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss(:,:,1),emis(:,1))
 
 c Dust conversion factor
@@ -1162 +1223 @@
 
 c Watercap
-      call interp_horiz (watercapold,watercaps,imold,jmold,iim,
-     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,watercaps,watercap)
+      call interp_horiz (watercapold(:,:,1),watercaps(:,:,1),
+     &       imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,watercaps(:,:,1),
+     &       watercap(:,1))
 
 c Surface albedo
-      call interp_horiz (albedoold,albedoS,imold,jmold,iim,
-     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,albedoS,albedo(:,1))
+      call interp_horiz (albedoold(:,:,1),albedoS(:,:,1),
+     &      imold,jmold,iim,jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,albedoS(:,:,1),albedo(:,1,1))
       
-      albedo(:,2)=albedo(:,1)
+      albedo(:,2,1)=albedo(:,1,1)
 
 c     write(46,*) 'emis',emis
@@ -1191 +1253 @@
          DO i=1,iim
             ptotal=ptotal+ps(i,j)*aire(i,j)/g
-            co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
+            co2icetotal = co2icetotal + co2iceS(i,j,1)*aire(i,j)
          END DO
       END DO
@@ -1217 +1279 @@
          DO j=1,jjp1
             DO i=1,iip1
-               co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
+               co2iceS(i,j,1)=co2iceS(i,j,1)*co2icetotalold/co2icetotal
             END DO
          END DO
@@ -1304 +1366 @@
          do j=1,jmold+1
            ! copy values
-           oldval(1)=tsurfold(i,j)
-           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
+           oldval(1)=tsurfold(i,j,1)
+           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold,1)
            ! build vertical coordinate
            oldgrid(1)=0. ! ground
@@ -1317 +1379 @@
      &                     mlayer,newval,nsoilmx)
          ! copy result in tsoilold
-         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
+         tsoiloldnew(i,j,1:nsoilmx,1)=newval(1:nsoilmx)
          enddo
         enddo
@@ -1336 +1398 @@
          do j=1,jmold+1
            ! copy values
-           oldval(1)=tsurfold(i,j)
-           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
+           oldval(1)=tsurfold(i,j,1)
+           oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold,1)
           ! interpolate
           call interp_line(oldgrid,oldval,nsoilold+1,
      &                     mlayer,newval,nsoilmx)
          ! copy result in tsoilold
-         tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
+         tsoiloldnew(i,j,1:nsoilmx,1)=newval(1:nsoilmx)
          enddo
         enddo
@@ -1352 +1414 @@
        
        else
-        tsoiloldnew(:,:,:)=tsoilold(:,:,:)
+        tsoiloldnew(:,:,:,1)=tsoilold(:,:,:,1)
        endif ! of if (depthinterpol)
       endif ! of if (olddepthdef)
 
       ! Do the horizontal interpolation
-       call interp_horiz(tsoiloldnew,tsoilS,
+       call interp_horiz(tsoiloldnew(:,:,:,1),tsoilS(:,:,:,1),
      &                  imold,jmold,iim,jjm,nsoilmx,
      &                   rlonuold,rlatvold,rlonu,rlatv)
 
       ! Reshape tsoilS to scalar grid as tsoil
-       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS,tsoil)
+       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS(:,:,:,1),
+     & tsoil(:,:,:))
 
 
@@ -1432 +1495 @@
 c traceurs surface
       do iq = 1, nqtot
-            call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
+            call interp_horiz(qsurfold(1,1,iq,1) ,qsurfs(1,1,iq,1),
      &                  imold,jmold,iim,jjm,1,
      &                  rlonuold,rlatvold,rlonu,rlatv)
       enddo
 
-      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs,qsurf)
+      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs(:,:,:,1),
+     &     qsurf(:,:,1))
 
 c traceurs 3D
@@ -1488 +1552 @@
       enddo
       
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,qsurf(:,igcm_co2))
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices(:,:,1),
+     & qsurf(:,igcm_co2,1))
 
 c-----------------------------------------------------------------------