Changeset 1395 for trunk/LMDZ.GENERIC

Timestamp:

Mar 12, 2015, 12:45:17 PM (10 years ago)

Author:

emillour

Message:

All GCMS:
Some cleanup and tidying on the dynamics/physics interface.
Essentially affects the "iniphysiq" routine in all physics packages.
EM

Location:

trunk/LMDZ.GENERIC/libf

Files:

• dyn3d/gcm.F (modified) (5 diffs)
• phystd/iniphysiq.F90 (modified) (6 diffs)
• phystd/newstart.F (modified) (1 diff)

trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F

@@ -5 +5 @@
       use control_mod, only: nday, day_step, iperiod, iphysiq,
      &                       iconser, ecritphy, idissip
-      use comgeomphy, only: initcomgeomphy
+!      use comgeomphy, only: initcomgeomphy
       IMPLICIT NONE
 
@@ -159 +159 @@
 c    variables pour l'initialisation de la physique :
 c    ------------------------------------------------
-      INTEGER ngridmx
-      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
-      REAL zcufi(ngridmx),zcvfi(ngridmx)
-      REAL latfi(ngridmx),lonfi(ngridmx)
-      REAL airefi(ngridmx)
-      SAVE latfi, lonfi, airefi
-      INTEGER i,j
+!      INTEGER ngridmx
+!      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
+!      REAL zcufi(ngridmx),zcvfi(ngridmx)
+!      REAL latfi(ngridmx),lonfi(ngridmx)
+!      REAL airefi(ngridmx)
+!      SAVE latfi, lonfi, airefi
+!      INTEGER i,j
 
 c-----------------------------------------------------------------------
@@ -183 +183 @@
 !#ifdef CPP_PHYS
       CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-      call initcomgeomphy
+!      call initcomgeomphy ! now done in iniphysiq
 !#endif
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -231 +231 @@
 
 !      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
-         latfi(1)=rlatu(1)
-         lonfi(1)=0.
-         zcufi(1) = cu(1)
-         zcvfi(1) = cv(1)
-         DO j=2,jjm
-            DO i=1,iim
-               latfi((j-2)*iim+1+i)= rlatu(j)
-               lonfi((j-2)*iim+1+i)= rlonv(i)
-               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
-               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
-            ENDDO
-         ENDDO
-         latfi(ngridmx)= rlatu(jjp1)
-         lonfi(ngridmx)= 0.
-         zcufi(ngridmx) = cu(ip1jm+1)
-         zcvfi(ngridmx) = cv(ip1jm-iim)
-
-         ! build airefi(), mesh area on physics grid
-         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
-         ! Poles are single points on physics grid
-         airefi(1)=airefi(1)*iim
-         airefi(ngridmx)=airefi(ngridmx)*iim
+!         latfi(1)=rlatu(1)
+!         lonfi(1)=0.
+!         zcufi(1) = cu(1)
+!         zcvfi(1) = cv(1)
+!         DO j=2,jjm
+!            DO i=1,iim
+!               latfi((j-2)*iim+1+i)= rlatu(j)
+!               lonfi((j-2)*iim+1+i)= rlonv(i)
+!               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
+!               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
+!            ENDDO
+!         ENDDO
+!         latfi(ngridmx)= rlatu(jjp1)
+!         lonfi(ngridmx)= 0.
+!         zcufi(ngridmx) = cu(ip1jm+1)
+!         zcvfi(ngridmx) = cv(ip1jm-iim)
+!
+!         ! build airefi(), mesh area on physics grid
+!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
+!         ! Poles are single points on physics grid
+!         airefi(1)=airefi(1)*iim
+!         airefi(ngridmx)=airefi(ngridmx)*iim
 
 ! Initialisation de la physique: pose probleme quand on tourne 
@@ -258 +258 @@
 ! Il faut une cle CPP_PHYS
 !#ifdef CPP_PHYS
-!         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
-         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys,
-     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
+!         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
+         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys,
+     &                rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,
      &                1)
 !     &                iflag_phys)

trunk/LMDZ.GENERIC/libf/phystd/iniphysiq.F90

@@ -1 @@
-subroutine iniphysiq(ngrid,nlayer, punjours, pdayref,ptimestep,          &
-                     plat,plon,parea,pcu,pcv,                            &
+subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
+                     rlatu,rlonv,aire,cu,cv,                             &
                      prad,pg,pr,pcpp,iflag_phys)
 
@@ -10 +10 @@
                                klon_omp_end, & ! end index of local omp subgrid
                                klon_mpi_begin ! start indes of columns (on local mpi grid)
-use comgeomphy, only : airephy, & ! physics grid area (m2)
+
+use comgeomphy, only : initcomgeomphy, &
+                       airephy, & ! physics grid area (m2)
                        cuphy, & ! cu coeff. (u_covariant = cu * u)
                        cvphy, & ! cv coeff. (v_covariant = cv * v)
@@ -21 +23 @@
 include "dimensions.h"
 include "comvert.h"
+include "iniprint.h"
 
 real,intent(in) :: prad ! radius of the planet (m)
@@ -27 +30 @@
 real,intent(in) :: pcpp ! specific heat Cp
 real,intent(in) :: punjours ! length (in s) of a standard day
-integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid)
+!integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid)
 integer,intent(in) :: nlayer ! number of atmospheric layers
-real,intent(in) :: plat(ngrid) ! latitudes of the physics grid
-real,intent(in) :: plon(ngrid) ! longitudes of the physics grid
-real,intent(in) :: parea(klon_glo) ! area (m2)
-real,intent(in) :: pcu(klon_glo) ! cu coeff. (u_covariant = cu * u)
-real,intent(in) :: pcv(klon_glo) ! cv coeff. (v_covariant = cv * v)
+integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
+integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
+real,intent(in) :: rlatu(jj+1) ! latitudes of the dynamics U grid
+real,intent(in) :: rlonv(ii+1) ! longitudes of the dynamics V grid
+real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
+real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
 integer,intent(in) :: pdayref ! reference day of for the simulation
 real,intent(in) :: ptimestep !physics time step (s)
@@ -39 +44 @@
 
 integer :: ibegin,iend,offset
+integer :: i,j
 character(len=20) :: modname='iniphysiq'
 character(len=80) :: abort_message
+real :: total_area_phy, total_area_dyn
+
+
+! global array, on full physics grid:
+real,allocatable :: latfi(:)
+real,allocatable :: lonfi(:)
+real,allocatable :: cufi(:)
+real,allocatable :: cvfi(:)
+real,allocatable :: airefi(:)
 
 IF (nlayer.NE.klev) THEN
@@ -51 +66 @@
 ENDIF
 
-IF (ngrid.NE.klon_glo) THEN
-  write(*,*) 'STOP in ',trim(modname)
-  write(*,*) 'Problem with dimensions :'
-  write(*,*) 'ngrid     = ',ngrid
-  write(*,*) 'klon   = ',klon_glo
-  abort_message = ''
-  CALL abort_gcm (modname,abort_message,1)
+!IF (ngrid.NE.klon_glo) THEN
+!  write(*,*) 'STOP in ',trim(modname)
+!  write(*,*) 'Problem with dimensions :'
+!  write(*,*) 'ngrid     = ',ngrid
+!  write(*,*) 'klon   = ',klon_glo
+!  abort_message = ''
+!  CALL abort_gcm (modname,abort_message,1)
+!ENDIF
+
+!call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
+
+! Generate global arrays on full physics grid
+allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
+latfi(1)=rlatu(1)
+lonfi(1)=0.
+cufi(1) = cu(1)
+cvfi(1) = cv(1)
+DO j=2,jj
+  DO i=1,ii
+    latfi((j-2)*ii+1+i)= rlatu(j)
+    lonfi((j-2)*ii+1+i)= rlonv(i)
+    cufi((j-2)*ii+1+i) = cu((j-1)*ii+1+i)
+    cvfi((j-2)*ii+1+i) = cv((j-1)*ii+1+i)
+  ENDDO
+ENDDO
+latfi(klon_glo)= rlatu(jj+1)
+lonfi(klon_glo)= 0.
+cufi(klon_glo) = cu((ii+1)*jj+1)
+cvfi(klon_glo) = cv((ii+1)*jj-ii)
+
+! build airefi(), mesh area on physics grid
+allocate(airefi(klon_glo))
+CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
+! Poles are single points on physics grid
+airefi(1)=sum(aire(1:ii,1))
+airefi(klon_glo)=sum(aire(1:ii,jj+1))
+
+! Sanity check: do total planet area match between physics and dynamics?
+total_area_dyn=sum(aire(1:ii,1:jj+1))
+total_area_phy=sum(airefi(1:klon_glo))
+IF (total_area_dyn/=total_area_phy) THEN
+  WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
+  WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
+  WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
+  IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
+    ! stop here if the relative difference is more than 0.001%
+    abort_message = 'planet total surface discrepancy'
+    CALL abort_gcm(modname, abort_message, 1)
+  ENDIF
 ENDIF
 
-!$OMP PARALLEL PRIVATE(ibegin,iend) & 
-        !$OMP SHARED(parea,pcu,pcv,plon,plat)
-      
+
+!$OMP PARALLEL 
+! Now generate local lon/lat/cu/cv/area arrays
+call initcomgeomphy
+
+!!!!$OMP PARALLEL PRIVATE(ibegin,iend) & 
+!!!     !$OMP SHARED(airefi,cufi,cvfi,lonfi,latfi)
+
 offset=klon_mpi_begin-1
-airephy(1:klon_omp)=parea(offset+klon_omp_begin:offset+klon_omp_end)
-cuphy(1:klon_omp)=pcu(offset+klon_omp_begin:offset+klon_omp_end)
-cvphy(1:klon_omp)=pcv(offset+klon_omp_begin:offset+klon_omp_end)
-rlond(1:klon_omp)=plon(offset+klon_omp_begin:offset+klon_omp_end)
-rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
+airephy(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end)
+cuphy(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end)
+cvphy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end)
+rlond(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end)
+rlatd(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end)
 
 ! copy over preff , ap() and bp() 

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

@@ -466 +466 @@
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
       ! Poles are single points on physics grid
-      airefi(1)=airefi(1)*iim
-      airefi(ngridmx)=airefi(ngridmx)*iim
+      airefi(1)=sum(aire(1:iim,1))
+      airefi(ngridmx)=sum(aire(1:iim,jjm+1))
 
 ! also initialize various physics flags/settings which might be needed