! ! $Id: startvar.F90 1952 2014-01-28 13:05:47Z musat $ ! !******************************************************************************* ! MODULE startvar ! !******************************************************************************* ! !------------------------------------------------------------------------------- ! Purpose: Access data from the database of atmospheric to initialize the model. !------------------------------------------------------------------------------- ! Comments: ! ! * This module is designed to work for Earth (and with ioipsl) ! ! * There are three ways to acces data, depending on the type of field ! which needs to be extracted. In any case the call should come after a restget ! and should be of the type : CALL startget(...) ! ! - A 1D variable on the physical grid : ! CALL startget_phys1d((varname, iml, jml, lon_in, lat_in, nbindex, & ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) ! ! - A 2D variable on the dynamical grid : ! CALL startget_phys2d(varname, iml, jml, lon_in, lat_in, & ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) ! ! - A 3D variable on the dynamical grid : ! CALL startget_dyn((varname, iml, jml, lon_in, lat_in, lml, pls, workvar, & ! champ, val_exp, jml2, lon_in2, lat_in2, ibar ) ! ! * Data needs to be in NetCDF format ! ! * Variables should have the following names in the files: ! 'RELIEF' : High resolution orography ! 'ST' : Surface temperature ! 'CDSW' : Soil moisture ! 'Z' : Surface geopotential ! 'SP' : Surface pressure ! 'U' : East ward wind ! 'V' : Northward wind ! 'TEMP' : Temperature ! 'R' : Relative humidity ! ! * There is a big mess with the longitude size. Should it be iml or iml+1 ? ! I have chosen to use the iml+1 as an argument to this routine and we declare ! internaly smaller fields when needed. This needs to be cleared once and for ! all in LMDZ. A convention is required. !------------------------------------------------------------------------------- #ifdef CPP_EARTH USE ioipsl IMPLICIT NONE PRIVATE PUBLIC startget_phys2d, startget_phys1d, startget_dyn ! INTERFACE startget ! MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn ! END INTERFACE REAL, SAVE :: deg2rad, pi INTEGER, SAVE :: iml_rel, jml_rel INTEGER, SAVE :: fid_phys, iml_phys, jml_phys INTEGER, SAVE :: fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel REAL, DIMENSION(:), ALLOCATABLE, TARGET, SAVE :: levdyn_ini REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: psol_dyn REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d CONTAINS !------------------------------------------------------------------------------- ! SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, & val_exp ,jml2, lon_in2, lat_in2, ibar) ! !------------------------------------------------------------------------------- ! Comment: ! This routine only works if the variable does not exist or is constant. !------------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: varname INTEGER, INTENT(IN) :: iml, jml REAL, DIMENSION(iml), INTENT(IN) :: lon_in REAL, DIMENSION(jml), INTENT(IN) :: lat_in INTEGER, INTENT(IN) :: nbindex REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ REAL, INTENT(IN) :: val_exp INTEGER, INTENT(IN) :: jml2 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 LOGICAL, INTENT(IN) :: ibar !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" REAL, DIMENSION(:,:), POINTER :: v2d !------------------------------------------------------------------------------- v2d=>NULL() IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN !--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES SELECT CASE(varname) CASE('tsol') IF(.NOT.ALLOCATED(tsol)) & CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) CASE('qsol') IF(.NOT.ALLOCATED(qsol)) & CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) CASE('psol') IF(.NOT.ALLOCATED(psol_dyn)) & CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval') IF(.NOT.ALLOCATED(relief)) & CALL start_init_orog(iml,jml,lon_in,lat_in) CASE('rads','snow','tslab','seaice','rugmer','agsno') CASE DEFAULT WRITE(lunout,*)'startget_phys1d' WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& //' from any data set'; STOP END SELECT !--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE SELECT CASE(varname) CASE('rads','snow','tslab','seaice'); champ=0.0 CASE('rugmer'); champ(:)=0.001 CASE('agsno'); champ(:)=50.0 CASE DEFAULT SELECT CASE(varname) CASE('tsol'); v2d=>tsol CASE('qsol'); v2d=>qsol CASE('psol'); v2d=>psol_dyn CASE('zmea'); v2d=>relief CASE('zstd'); v2d=>zstd CASE('zsig'); v2d=>zsig CASE('zgam'); v2d=>zgam CASE('zthe'); v2d=>zthe CASE('zpic'); v2d=>zpic CASE('zval'); v2d=>zval END SELECT IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size' STOP END IF CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ) END SELECT ELSE !--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION SELECT CASE(varname) CASE('tsol') IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml)) CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol) END SELECT END IF END SUBROUTINE startget_phys1d ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp, & jml2, lon_in2, lat_in2 , ibar, msk) ! !------------------------------------------------------------------------------- ! Comment: ! This routine only works if the variable does not exist or is constant. !------------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: varname INTEGER, INTENT(IN) :: iml, jml REAL, DIMENSION(iml), INTENT(IN) :: lon_in REAL, DIMENSION(jml), INTENT(IN) :: lat_in REAL, DIMENSION(iml,jml), INTENT(INOUT) :: champ REAL, INTENT(IN) :: val_exp INTEGER, INTENT(IN) :: jml2 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 LOGICAL, INTENT(IN) :: ibar REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" REAL, DIMENSION(:,:), POINTER :: v2d=>NULL() LOGICAL :: lrelief1, lrelief2 !------------------------------------------------------------------------------- v2d=>NULL() lrelief1=(.NOT.ALLOCATED(relief).AND. PRESENT(msk)) lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk)) IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN !--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES SELECT CASE(varname) CASE('psol') IF(.NOT.ALLOCATED(psol_dyn)) & CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) CASE('relief') IF(lrelief1) CALL start_init_orog(iml,jml,lon_in,lat_in,msk) IF(lrelief2) CALL start_init_orog(iml,jml,lon_in,lat_in) CASE('surfgeo') IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in) CASE('rugosite','masque') IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in) CASE DEFAULT WRITE(lunout,*)'startget_phys2d' WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& //' from any data set'; STOP END SELECT !--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE SELECT CASE(varname) CASE('psol'); v2d=>psol_dyn CASE('relief'); v2d=>relief CASE('rugosite'); v2d=>rugo CASE('masque'); v2d=>masque CASE('surfgeo'); v2d=>phis END SELECT IF(SIZE(champ)/=SIZE(v2d)) THEN WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size' STOP END IF champ(:,:)=v2d(:,:) ELSE !--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION SELECT CASE(varname) CASE ('surfgeo') IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml)) IF(SIZE(phis)/=SIZE(champ)) THEN WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size' STOP END IF phis(:,:)=champ(:,:) END SELECT END IF END SUBROUTINE startget_phys2d ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE startget_dyn(varname, lon_in, lat_in, pls,workvar,& champ, val_exp, lon_in2, lat_in2, ibar) use assert_eq_m, only: assert_eq !------------------------------------------------------------------------------- ! Comment: ! This routine only works if the variable does not exist or is constant. !------------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: varname REAL, INTENT(IN) :: lon_in(:) ! dim(iml) REAL, INTENT(IN) :: lat_in(:) ! dim(jml) REAL, INTENT(IN) :: pls(:, :, :) ! dim(iml, jml, lml) REAL, INTENT(IN) :: workvar(:, :, :) ! dim(iml, jml, lml) REAL, INTENT(INOUT) :: champ(:, :, :) ! dim(iml, jml, lml) REAL, INTENT(IN) :: val_exp REAL, INTENT(IN) :: lon_in2(:) ! dim(iml) REAL, INTENT(IN) :: lat_in2(:) ! dim(jml2) LOGICAL, INTENT(IN) :: ibar !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" #include "dimensions.h" #include "comconst.h" #include "paramet.h" #include "comgeom2.h" INTEGER :: iml, jml INTEGER :: lml INTEGER :: jml2 REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL() CHARACTER(LEN=10) :: vname INTEGER :: il REAL :: xppn, xpps !------------------------------------------------------------------------------- NULLIFY(v3d) IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1), & & size(champ, 1), size(lon_in2)/), "startget_dyn iml") jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2), & & size(champ, 2), "startget_dyn jml") lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3), & & "startget_dyn lml") jml2 = size(lat_in2) !--- READING UNALLOCATED FILES IF(.NOT.ALLOCATED(psol_dyn)) & CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) !--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS SELECT CASE(varname) CASE('u'); vname='U' CASE('v'); vname='V' CASE('t','tpot'); vname='TEMP' CASE('q'); vname='R' CASE DEFAULT WRITE(lunout,*)'startget_dyn' WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)& //' from any data set'; STOP END SELECT CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2, & lon_in2, lat_in2, pls, champ,ibar ) !--- COMPUTING THE REQUIRED FILED SELECT CASE(varname) CASE('u') !--- Eastward wind DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO champ(iml,:,:)=champ(1,:,:) CASE('v') !--- Northward wind DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO champ(iml,:,:)=champ(1,:,:) CASE('tpot') !--- Temperature IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN champ=champ*cpp/workvar DO il=1,lml xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols champ(:,1 ,il) = xppn champ(:,jml,il) = xpps END DO ELSE WRITE(lunout,*)'Could not compute potential temperature as the' WRITE(lunout,*)'Exner function is missing or constant.'; STOP END IF CASE('q') !--- Relat. humidity IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN champ=0.01*champ*workvar WHERE(champ<0.) champ=1.0E-10 DO il=1,lml xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols champ(:,1 ,il) = xppn champ(:,jml,il) = xpps END DO ELSE WRITE(lunout,*)'Could not compute specific humidity as the' WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP END IF END SELECT END IF END SUBROUTINE startget_dyn ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu) ! !------------------------------------------------------------------------------- ! Arguments: INTEGER, INTENT(IN) :: iml, jml REAL, DIMENSION(iml), INTENT(IN) :: lon_in REAL, DIMENSION(jml), INTENT(IN) :: lat_in REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" CHARACTER(LEN=25) :: title CHARACTER(LEN=120) :: orofname LOGICAL :: check=.TRUE. REAL, DIMENSION(1) :: lev REAL :: date, dt INTEGER, DIMENSION(1) :: itau INTEGER :: fid, llm_tmp, ttm_tmp REAL, DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini !------------------------------------------------------------------------------- pi=2.0*ASIN(1.0); deg2rad=pi/180.0 orofname = 'Relief.nc'; title='RELIEF' IF(check) WRITE(lunout,*)'Reading the high resolution orography' CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& lev, ttm_tmp, itau, date, dt, fid) ALLOCATE(relief_hi(iml_rel,jml_rel)) CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) CALL flinclo(fid) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, & relief_hi, .FALSE.) DEALLOCATE(lon_ini,lat_ini) !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra& &g code' ALLOCATE(phis(iml,jml)) ! Geopotentiel au sol ALLOCATE(zstd(iml,jml)) ! Deviation standard de l'orographie sous-maille ALLOCATE(zsig(iml,jml)) ! Pente de l'orographie sous-maille ALLOCATE(zgam(iml,jml)) ! Anisotropie de l'orographie sous maille ALLOCATE(zthe(iml,jml)) ! Orientation axe +grande pente d'oro sous maille ALLOCATE(zpic(iml,jml)) ! Hauteur pics de la SSO ALLOCATE(zval(iml,jml)) ! Hauteur vallees de la SSO ALLOCATE(relief(iml,jml)) ! Orographie moyenne ALLOCATE(masque(iml,jml)) ! Masque terre ocean masque = -99999. IF(PRESENT(masque_lu)) masque=masque_lu CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, & lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque) phis = phis * 9.81 !--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! ) IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography' ALLOCATE(rugo (iml ,jml)) ALLOCATE(tmp_var(iml-1,jml)) CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, & lon_in, lat_in, tmp_var) rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:) DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad) RETURN END SUBROUTINE start_init_orog ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) ! !------------------------------------------------------------------------------- ! Arguments: INTEGER, INTENT(IN) :: iml, jml REAL, DIMENSION(iml), INTENT(IN) :: lon_in REAL, DIMENSION(jml), INTENT(IN) :: lat_in INTEGER, INTENT(IN) :: jml2 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 LOGICAL, INTENT(IN) :: ibar !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" CHARACTER(LEN=25) :: title CHARACTER(LEN=120) :: physfname LOGICAL :: check=.TRUE. REAL :: date, dt INTEGER, DIMENSION(1) :: itau INTEGER :: llm_tmp, ttm_tmp REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini !------------------------------------------------------------------------------- physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0 IF(check) WRITE(lunout,*)'Opening the surface analysis' CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys) ALLOCATE(lat_phys(iml_phys,jml_phys)) ALLOCATE(lon_phys(iml_phys,jml_phys)) ALLOCATE(levphys_ini(llm_tmp)) CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys, & lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys) DEALLOCATE(levphys_ini) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys)) lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) !--- SURFACE TEMPERATURE title='ST' ALLOCATE(tsol(iml,jml)) CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana) CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,& var_ana , ibar ) CALL interp_startvar(title, ibar, .TRUE., & iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & lon_in, lat_in, lon_in2, lat_in2, tsol) !--- SOIL MOISTURE title='CDSW' ALLOCATE(qsol(iml,jml)) CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana) CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,& var_ana, ibar ) CALL interp_startvar(title, ibar, .TRUE., & iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & lon_in, lat_in, lon_in2, lat_in2, qsol) CALL flinclo(fid_phys) DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) END SUBROUTINE start_init_phys ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) ! !------------------------------------------------------------------------------- ! Arguments: INTEGER, INTENT(IN) :: iml, jml REAL, DIMENSION(iml), INTENT(IN) :: lon_in REAL, DIMENSION(jml), INTENT(IN) :: lat_in INTEGER, INTENT(IN) :: jml2 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 LOGICAL, INTENT(IN) :: ibar !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" #include "dimensions.h" #include "paramet.h" #include "comgeom2.h" CHARACTER(LEN=25) :: title CHARACTER(LEN=120) :: physfname LOGICAL :: check=.TRUE. REAL :: date, dt INTEGER, DIMENSION(1) :: itau INTEGER :: i, j REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana, z REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini REAL, DIMENSION(:), ALLOCATABLE :: xppn, xpps !------------------------------------------------------------------------------- !--- KINETIC ENERGY physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0 IF(check) WRITE(lunout,*) 'Opening the surface analysis' CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn ALLOCATE(lat_dyn(iml_dyn,jml_dyn)) ALLOCATE(lon_dyn(iml_dyn,jml_dyn)) ALLOCATE(levdyn_ini(llm_dyn)) CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,& levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) !--- SURFACE GEOPOTENTIAL title='Z' ALLOCATE(z(iml,jml)) CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana) CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, & var_ana, ibar) CALL interp_startvar(title, ibar, .TRUE., & iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & lon_in, lat_in, lon_in2, lat_in2, z) !--- SURFACE PRESSURE title='SP' ALLOCATE(psol_dyn(iml,jml)) CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana) CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, & var_ana, ibar) CALL interp_startvar(title, ibar, .TRUE., & iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, & lon_in, lat_in, lon_in2, lat_in2, psol_dyn) DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) !--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE IF(.NOT.ALLOCATED(tsol)) THEN CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar) ELSE IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN WRITE(lunout,*)'start_init_dyn :' WRITE(lunout,*)'The temperature field we have does not have the right size' STOP END IF END IF IF(.NOT.ALLOCATED(phis)) THEN CALL start_init_orog(iml,jml,lon_in,lat_in) ELSE IF(SIZE(phis)/=SIZE(psol_dyn)) THEN WRITE(lunout,*)'start_init_dyn :' WRITE(lunout,*)'The orography field we have does not have the right size' STOP END IF END IF !--- PSOL IS COMPUTED IN PASCALS DO j = 1, jml DO i = 1, iml-1 psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j)) END DO psol_dyn(iml,j) = psol_dyn(1,j) END DO DEALLOCATE(z) ALLOCATE(xppn(iml-1),xpps(iml-1)) DO i = 1, iml-1 xppn(i) = aire( i,1) * psol_dyn( i,1) xpps(i) = aire( i,jml) * psol_dyn( i,jml) END DO DO i = 1, iml psol_dyn(i,1 ) = SUM(xppn)/apoln psol_dyn(i,jml) = SUM(xpps)/apols END DO DEALLOCATE(xppn,xpps) RETURN END SUBROUTINE start_init_dyn ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, & lon_in2, lat_in2, pls_in, var3d, ibar) use pchsp_95_m, only: pchsp_95 use pchfe_95_m, only: pchfe_95 ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: varname INTEGER, INTENT(IN) :: iml, jml, lml REAL, DIMENSION(iml), INTENT(IN) :: lon_in REAL, DIMENSION(jml), INTENT(IN) :: lat_in INTEGER, INTENT(IN) :: jml2 REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 REAL, DIMENSION(iml, jml, lml), INTENT(IN) :: pls_in REAL, DIMENSION(iml, jml, lml), INTENT(OUT) :: var3d LOGICAL, INTENT(IN) :: ibar !---------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" LOGICAL:: check=.TRUE., skip REAL chmin, chmax INTEGER ii, ij, il, ierr integer n_extrap ! number of extrapolated points REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder !--------------------------------------------------------------------------- IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D field.' IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, & ttm_dyn IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, & llm_dyn IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, & var_ana3d) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn)) CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, & levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar) DEALLOCATE(lon_ini, lat_ini) !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro DO il=1, llm_dyn CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, & lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, & lon_in2, lat_in2, var_tmp3d(:, :, il)) END DO DEALLOCATE(lon_rad, lat_rad) !--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND ax = lev_dyn(llm_dyn:1:-1) skip = .false. n_extrap = 0 DO ij=1, jml DO ii=1, iml-1 ay = var_tmp3d(ii, ij, llm_dyn:1:-1) yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.) CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & var3d(ii, ij, lml:1:-1), ierr) if (ierr < 0) stop 1 n_extrap = n_extrap + ierr END DO END DO if (n_extrap /= 0) then print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap end if var3d(iml, :, :) = var3d(1, :, :) DO il=1, lml CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax) WRITE(lunout, *)' '//TRIM(varname)//' min max l ', il, chmin, chmax END DO END SUBROUTINE start_inter_3d ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj, lon, lat, vari, & i1, j1, j2, lon1, lat1, lon2, lat2, varo) ! !------------------------------------------------------------------------------- USE inter_barxy_m, only: inter_barxy ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: vname LOGICAL, INTENT(IN) :: ibar, ibeg INTEGER, INTENT(IN) :: ii, jj REAL, DIMENSION(ii), INTENT(IN) :: lon REAL, DIMENSION(jj), INTENT(IN) :: lat REAL, DIMENSION(ii,jj), INTENT(IN) :: vari INTEGER, INTENT(IN) :: i1, j1, j2 REAL, DIMENSION(i1), INTENT(IN) :: lon1 REAL, DIMENSION(j1), INTENT(IN) :: lat1 REAL, DIMENSION(i1), INTENT(IN) :: lon2 REAL, DIMENSION(j2), INTENT(IN) :: lat2 REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo !------------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" REAL, DIMENSION(i1-1,j1) :: vtmp !------------------------------------------------------------------------------- IF(ibar) THEN IF(ibeg) THEN WRITE(lunout,*) & '---------------------------------------------------------------' WRITE(lunout,*) & '$$$ Utilisation de l interpolation barycentrique pour '//TRIM(vname)//' $$$' WRITE(lunout,*) & '---------------------------------------------------------------' END IF CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp) ELSE CALL grille_m (ii, jj, lon, lat, vari, i1-1, j1, lon1, lat1, vtmp) END IF CALL gr_int_dyn(vtmp, varo, i1-1, j1) END SUBROUTINE interp_startvar ! !------------------------------------------------------------------------------- #endif ! of #ifdef CPP_EARTH END MODULE startvar ! !*******************************************************************************