Changeset 1395 for trunk/LMDZ.MARS/libf/phymars/iniphysiq.F90

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

File:

• trunk/LMDZ.MARS/libf/phymars/iniphysiq.F90 (modified) (6 diffs)

trunk/LMDZ.MARS/libf/phymars/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)
@@ -20 +22 @@
 implicit none
 
+include "iniprint.h"
+
 real,intent(in) :: prad ! radius of the planet (m)
 real,intent(in) :: pg ! gravitational acceleration (m/s2)
@@ -25 +29 @@
 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)
@@ -37 +43 @@
 
 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
@@ -49 +65 @@
 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
+
+! 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(parea,pcu,pcv,plon,plat)
       
 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 some fundamental parameters to physics