source: trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/iniphysiq_mod.F90 @ 1540

MODULE iniphysiq_mod

CONTAINS

subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
               rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn,  &
                     prad,pg,pr,pcpp,iflag_phys)

use dimphy, only : klev ! number of atmospheric levels
use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns
                                       ! (on full grid)
use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid)
                               klon_omp_begin, & ! start index of local omp subgrid
                               klon_omp_end, & ! end index of local omp subgrid
                               klon_mpi_begin ! start indes of columns (on local mpi grid)
use control_mod, only: nday 
use comgeomphy, only : initcomgeomphy, &
                       airephy, & ! physics grid area (m2)
                       cuphy, & ! cu coeff. (u_covariant = cu * u)
                       cvphy, & ! cv coeff. (v_covariant = cv * v)
                       rlond, & ! longitudes
                       rlatd ! latitudes
use surf_heat_transp_mod, only: ini_surf_heat_transp
use infotrac, only : nqtot ! number of advected tracers
use planete_mod, only: ini_planete_mod
USE comvert_mod, ONLY: ap,bp,preff
use inifis_mod, only: inifis
use regular_lonlat_mod, only: init_regular_lonlat, &
                              east, west, north, south, &
                              north_east, north_west, &
                              south_west, south_east
use ioipsl_getin_p_mod, only: getin_p

implicit none
include "dimensions.h"
include "paramet.h"
include "comgeom.h"
include "iniprint.h"

real,intent(in) :: prad ! radius of the planet (m)
real,intent(in) :: pg ! gravitational acceleration (m/s2)
real,intent(in) :: pr ! ! reduced gas constant R/mu
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) :: nlayer ! number of atmospheric layers
integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
real,intent(in) :: rlatudyn(jj+1) ! latitudes of the physics grid
real,intent(in) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
real,intent(in) :: rlonvdyn(ii+1) ! longitudes of the physics grid
real,intent(in) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
real,intent(in) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
real,intent(in) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
real,intent(in) :: cvdyn((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)
integer,intent(in) :: iflag_phys ! type of physics to be called

integer :: ibegin,iend,offset
integer :: i,j
character(len=20) :: modname='iniphysiq'
character(len=80) :: abort_message
real :: total_area_phy, total_area_dyn
real :: pi
logical :: ok_slab_ocean

! boundaries, on global grid
real,allocatable :: boundslon_reg(:,:)
real,allocatable :: boundslat_reg(:,:)

! global array, on full physics grid:
real,allocatable :: latfi(:)
real,allocatable :: lonfi(:)
real,allocatable :: cufi(:)
real,allocatable :: cvfi(:)
real,allocatable :: airefi(:)

pi=2.*asin(1.0)

IF (nlayer.NE.klev) THEN
  write(*,*) 'STOP in ',trim(modname)
  write(*,*) 'Problem with dimensions :'
  write(*,*) 'nlayer     = ',nlayer
  write(*,*) 'klev   = ',klev
  abort_message = ''
  CALL abort_gcm (modname,abort_message,1)
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)
!ENDIF

!call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))

! init regular global longitude-latitude grid points and boundaries
ALLOCATE(boundslon_reg(ii,2))
ALLOCATE(boundslat_reg(jj+1,2))
  
DO i=1,ii
   boundslon_reg(i,east)=rlonudyn(i) 
   boundslon_reg(i,west)=rlonudyn(i+1) 
ENDDO

boundslat_reg(1,north)= PI/2 
boundslat_reg(1,south)= rlatvdyn(1)
DO j=2,jj
   boundslat_reg(j,north)=rlatvdyn(j-1) 
   boundslat_reg(j,south)=rlatvdyn(j) 
ENDDO
boundslat_reg(jj+1,north)= rlatvdyn(jj) 
boundslat_reg(jj+1,south)= -PI/2

! Write values in module regular_lonlat_mod
CALL init_regular_lonlat(ii,jj+1, rlonvdyn(1:ii), rlatudyn, &
                         boundslon_reg, boundslat_reg)

! Generate global arrays on full physics grid
allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
latfi(1)=rlatudyn(1)
lonfi(1)=0.
cufi(1) = cudyn(1)
cvfi(1) = cvdyn(1)
DO j=2,jj
  DO i=1,ii
    latfi((j-2)*ii+1+i)= rlatudyn(j)
    lonfi((j-2)*ii+1+i)= rlonvdyn(i)
    cufi((j-2)*ii+1+i) = cudyn((j-1)*(ii+1)+i)
    cvfi((j-2)*ii+1+i) = cvdyn((j-1)*(ii+1)+i)
  ENDDO
ENDDO
latfi(klon_glo)= rlatudyn(jj+1)
lonfi(klon_glo)= 0.
cufi(klon_glo) = cudyn((ii+1)*jj+1)
cvfi(klon_glo) = cvdyn((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,airedyn,airefi)
! Poles are single points on physics grid
airefi(1)=sum(airedyn(1:ii,1))
airefi(klon_glo)=sum(airedyn(1:ii,jj+1))

! Sanity check: do total planet area match between physics and dynamics?
total_area_dyn=sum(airedyn(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 
! 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)=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() 
call ini_planete_mod(nlayer,preff,ap,bp)

! for slab ocean, copy over some arrays
ok_slab_ocean=.false. ! default value
call getin_p("ok_slab_ocean",ok_slab_ocean)
if (ok_slab_ocean) then
  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
                            aireu,airev)
endif

! copy some fundamental parameters to physics 
! and do some initializations 
call inifis(klon_omp,nlayer,nqtot,pdayref,punjours,nday,ptimestep, &
            rlatd,rlond,airephy,prad,pg,pr,pcpp)

!$OMP END PARALLEL


end subroutine iniphysiq


END MODULE iniphysiq_mod