Changeset 1529 for trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F

Timestamp:

Apr 5, 2016, 10:51:51 AM (9 years ago)

Author:

emillour

Message:

Generic GCM: Towards a cleaner separation between physics and dynamics

• Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
• Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
• Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

EM

File:

• trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F (modified) (17 diffs)

trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F

@@ -40 +40 @@
 !=================================================================
       use surfdat_h, only: phisfi
+      use comgeomphy, only: airephy
       use time_phylmdz_mod, only: ecritphy, day_step, iphysiq, day_ini
       USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
      &                               is_master, gather
-      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat, nbp_lev
       implicit none
 
 ! Commons
-      include "dimensions.h"
-      include "paramet.h"
-      include "comgeom.h"
       include "netcdf.inc"
 
@@ -56 +55 @@
       character (len=*),intent(in) :: nom,titre,unite
       integer,intent(in) :: dim
-      real,intent(in) :: px(ngrid,llm)
+      real,intent(in) :: px(ngrid,nbp_lev)
 
 ! Local variables:
 
-      real*4 dx3(iip1,jjp1,llm) ! to store a 3D data set
-      real*4 dx2(iip1,jjp1)     ! to store a 2D (surface) data set
-      real*4 dx1(llm)           ! to store a 1D (column) data set
+      real*4 dx3(nbp_lon+1,nbp_lat,nbp_lev) ! to store a 3D data set
+      real*4 dx2(nbp_lon+1,nbp_lat)     ! to store a 2D (surface) data set
+      real*4 dx1(nbp_lev)           ! to store a 1D (column) data set
       real*4 dx0
 
@@ -68 +67 @@
 !$OMP THREADPRIVATE(date)
 
-      REAL phis(ip1jmp1)
+      REAL phis((nbp_lon+1),nbp_lat)
+      REAL area((nbp_lon+1),nbp_lat)
 
       integer irythme
       integer ierr,ierr2
-      integer iq
-      integer i,j,l,zmax , ig0
+      integer i,j,l, ig0
 
       integer,save :: zitau=0
@@ -102 +101 @@
 #ifdef CPP_PARA
 ! Added to work in parallel mode
-      real dx3_glop(klon_glo,llm)
-      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
+      real dx3_glop(klon_glo,nbp_lev)
+      real dx3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set
       real dx2_glop(klon_glo)
-      real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+      real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
       real px2(ngrid)
-!      real dx1_glo(llm)          ! to store a 1D (column) data set
+!      real dx1_glo(nbp_lev)          ! to store a 1D (column) data set
 !      real dx0_glo
       real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
 #else
       real phisfi_glo(ngrid) ! surface geopotential on global physics grid
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -117 +118 @@
 !Sortie des variables au rythme voulu
 
-      irythme = ecritphy ! sortie au rythme de ecritphy
+      irythme = int(ecritphy) ! sortie au rythme de ecritphy
 !     irythme = iconser  ! sortie au rythme des variables de controle
 !     irythme = iphysiq  ! sortie a tous les pas physique
@@ -188 +189 @@
           ! Gather phisfi() geopotential on physics grid
           call Gather(phisfi,phisfi_glo)
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
 #else
          phisfi_glo(:)=phisfi(:)
+         areafi_glo(:)=airephy(:)
 #endif
 
@@ -216 +220 @@
          ierr = NF_ENDDEF(nid)
 
+         ! Build phis() and area()
+         do i=1,nbp_lon+1 ! poles
+           phis(i,1)=phisfi_glo(1)
+           phis(i,nbp_lat)=phisfi_glo(klon_glo)
+           ! for area, divide at the poles by nbp_lon
+           area(i,1)=areafi_glo(1)/nbp_lon
+           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+         enddo
+         do j=2,nbp_lat-1
+           ig0= 1+(j-2)*nbp_lon
+           do i=1,nbp_lon
+              phis(i,j)=phisfi_glo(ig0+i)
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           phis(nbp_lon+1,j)=phis(1,j)
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+         
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
-         call iniwrite(nid,day_ini,phis)
+         call iniwrite(nid,day_ini,phis,area)
 
          endif ! of if (is_master)
@@ -234 +256 @@
 
       if (ngrid.eq.1) then
-        ! in testphys1d, for the 1d version of the GCM, iphysiq 
-        ! should be most likely 1 (because no dyn!)
+        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
+        ! are undefined; so set them to 1
         iphysiq=1 
+        irythme=1
         ! NB: 
       endif
@@ -296 +319 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -304 +327 @@
 !         Passage variable physique -->  variable dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-           DO l=1,llm
-             DO i=1,iip1
+           DO l=1,nbp_lev
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -344 +367 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
-           edges(3)=llm
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
+           edges(3)=nbp_lev
            edges(4)=1
 !#ifdef NC_DOUBLE
@@ -380 +403 @@
             call Grid1Dto2D_glo(dx2_glop,dx2_glo)
             ! copy dx2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=dx2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -390 +413 @@
 !         recast (copy) variable from physics grid to dynamics grid
 
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx2(i,1)=px(1,1)
-                dx2(i,jjp1)=px(ngrid,1)
+                dx2(i,nbp_lat)=px(ngrid,1)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx2(i,j)=px(ig0+i,1)
                 ENDDO
-                dx2(iip1,j)=dx2(1,j)
+                dx2(nbp_lon+1,j)=dx2(1,j)
              ENDDO
 #endif
@@ -426 +449 @@
            corner(2)=1
            corner(3)=ntime
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=1
 
@@ -457 +480 @@
 !         Passage variable physique -->  physique dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-          do l=1,llm
+          do l=1,nbp_lev
             dx1(l)=px(1,l)
           enddo
@@ -479 +502 @@
            corner(2)=ntime
            
-           edges(1)=llm
+           edges(1)=nbp_lev
            edges(2)=1
 !#ifdef NC_DOUBLE
@@ -543 +566 @@
 
 #endif
+! of #ifndef MESOSCALE
       end