Index: /trunk/DOC/chantiers/compilation.txt =================================================================== --- /trunk/DOC/chantiers/compilation.txt (revision 848) +++ /trunk/DOC/chantiers/compilation.txt (revision 849) @@ -21,2 +21,7 @@ Compils ok. +Ehouarn (22/11/2012): +- OK pour faire qu'on puisse compiler des exécutables situés dans phy* +- On est d'accord pour n'utiliser plus que "makelmdz" et "makelmdz_fcm" +- reste à faire que makelmdz_fcm soit capable de compiler des sources + Fortran et C (pour Titan) Index: /trunk/DOC/chantiers/meschantiers-Ehouarn.txt =================================================================== --- /trunk/DOC/chantiers/meschantiers-Ehouarn.txt (revision 848) +++ /trunk/DOC/chantiers/meschantiers-Ehouarn.txt (revision 849) @@ -16,5 +16,10 @@ "Terrestre" rev. 1678) avec la rev 841. Reste à faire l'exercice de validation des modifs en //... + Avec rev 841 et 843 Cp(T) en // OK + Reste toutefois à paufiner les choses au niveau de calfis_p et des + domaines d'utilisation de tpot2t_p et t2tpot_p (toutes les taches MPI + font l'intégralité du travail!) - Faire qu'on puisse compiler un exécutable qui serait dans la physique - avec makelmdz & makelmdz_fcm. + avec makelmdz & makelmdz_fcm. OK avec la rev 849. + Manque encore la possibilité de compiler des sources en C avec makelmdz_fcm Index: /trunk/LMDZ.COMMON/libf/dyn3d/calfis.F =================================================================== --- /trunk/LMDZ.COMMON/libf/dyn3d/calfis.F (revision 848) +++ /trunk/LMDZ.COMMON/libf/dyn3d/calfis.F (revision 849) @@ -462,5 +462,5 @@ c if (planet_type=="earth") then -#ifdef CPP_PHYS +#ifdef CPP_EARTH ! PVtheta calls tetalevel, which is in the (Earth) physics cIM calcul PV a teta=350, 380, 405K Index: unk/LMDZ.COMMON/libf/dyn3d/newstart-VT.F =================================================================== --- /trunk/LMDZ.COMMON/libf/dyn3d/newstart-VT.F (revision 848) +++ (revision ) @@ -1,1104 +1,0 @@ -C====================================================================== - PROGRAM newstart -c======================================================================= -c -c -c Auteur: S. Lebonnois, -c a partir des newstart/start_archive/lect_start_archive martiens -c -c Derniere modif : 02/09 (ecriture des q*) -c 01/12 (inclusion dans svn dyn3d) -c -c Objet: Modify the grid for the initial state (LMD GCM VENUS/TITAN) -c ----- (from file NetCDF start_archive.nc) -c -c -c======================================================================= - - use IOIPSL - USE filtreg_mod - USE startvar - USE control_mod - USE infotrac - - implicit none - -#include "dimensions.h" -#include "paramet.h" -#include "comconst.h" -#include "comdissnew.h" -#include "comvert.h" -#include "comgeom2.h" -#include "logic.h" -#include "temps.h" -#include "ener.h" -#include "description.h" -#include "serre.h" -#include "dimsoil.h" -#include "netcdf.inc" - -c----------------------------------------------------------------------- -c Declarations -c----------------------------------------------------------------------- - -c Variables pour fichier "ini" -c------------------------------------ - INTEGER imold,jmold,lmold,nqold,ip1jmp1old - INTEGER length - parameter (length = 100) - real tab_cntrl(2*length) - INTEGER isoil,iq,iqmax - CHARACTER*2 str2 - -c Variable histoire -c------------------ - REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants - REAL teta(iip1,jjp1,llm),ps(iip1,jjp1) - REAL phis(iip1,jjp1) ! geopotentiel au sol - REAL masse(ip1jmp1,llm) ! masse de l'atmosphere - REAL, ALLOCATABLE, DIMENSION(:,:,:,:):: q! champs advectes - REAL tab_cntrl_dyn(length) ! tableau des parametres de start - -c variable physique -c------------------ - integer ngridmx - parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm)) - REAL tab_cntrl_fi(length) ! tableau des parametres de startfi - real rlat(ngridmx),rlon(ngridmx) - REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx) - REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx) - real solsw(ngridmx),dlw(ngridmx) - REAL zmea(ngridmx), zstd(ngridmx) - REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx) - REAL zpic(ngridmx), zval(ngridmx) - real t_fi(ngridmx,llm) - -c Variable nouvelle grille naturelle au point scalaire -c------------------------------------------------------ - real us(iip1,jjp1,llm),vs(iip1,jjp1,llm) - REAL p3d(iip1,jjp1,llm+1) ! pression aux interfaces - REAL phisold_newgrid(iip1,jjp1) - REAL T(iip1,jjp1,llm) - real rlonS(iip1,jjp1),rlatS(iip1,jjp1) - real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx) - real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1) - real solswS(ip1jmp1),dlwS(ip1jmp1) - real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1) - real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1) - real zvalS(ip1jmp1) - - real ptotal - -c Var intermediaires : vent naturel, mais pas coord scalaire -c----------------------------------------------------------- - real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) - - REAL pks(iip1,jjp1) ! exner (f pour filtre) - REAL pk(iip1,jjp1,llm) - REAL pkf(iip1,jjp1,llm) - REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) - - -c Variable de l'ancienne grille -c--------------------------------------------------------- - - real, dimension(:), allocatable :: rlonuold, rlatvold - real, dimension(:), allocatable :: rlonvold, rlatuold - real, dimension(:), allocatable :: nivsigsold,nivsigold - real, dimension(:), allocatable :: apold,bpold - real, dimension(:), allocatable :: presnivsold - real, dimension(:,:,:), allocatable :: uold,vold,told - real, dimension(:,:,:,:), allocatable :: qold - real, dimension(:,:,:), allocatable :: tsoilold - real, dimension(:,:), allocatable :: psold,phisold - real, dimension(:,:), allocatable :: tsurfold - real, dimension(:,:), allocatable :: albeold,radsolold - real, dimension(:,:), allocatable :: sollwold,solswold - real, dimension(:,:), allocatable :: dlwold - real, dimension(:,:), allocatable :: zmeaold,zstdold,zsigold - real, dimension(:,:), allocatable :: zgamold,ztheold,zpicold - real, dimension(:,:), allocatable :: zvalold - - real ptotalold - -c Variable intermediaires iutilise pour l'extrapolation verticale -c---------------------------------------------------------------- - real, dimension(:,:,:), allocatable :: var,varp1 - -c divers local -c----------------- - - integer ierr,nid,nvarid - INTEGER ij, l,i,j - character*80 fichnom - integer, dimension(4) :: start,counter - REAL phisinverse(iip1,jjp1) ! geopotentiel au sol avant inversion - logical topoflag,albedoflag - real albedo - -c======================================================================= -c INITIALISATIONS DIVERSES -c======================================================================= - -c VENUS/TITAN - - iflag_trac = 1 -c----------------------------------------------------------------------- -c Initialisation des traceurs -c --------------------------- -c Choix du nombre de traceurs et du schema pour l'advection -c dans fichier traceur.def, par default ou via INCA - call infotrac_init - -c Allocation de la tableau q : champs advectes - allocate(q(iip1,jjp1,llm,nqtot)) - -c----------------------------------------------------------------------- -c Ouverture du fichier a modifier (start_archive.nc) -c----------------------------------------------------------------------- - - write(*,*) 'Creation d un etat initial a partir de' - write(*,*) './start_archive.nc' - write(*,*) - fichnom = 'start_archive.nc' - ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) - IF (ierr.NE.NF_NOERR) THEN - write(6,*)' Pb d''ouverture du fichier ',fichnom - write(6,*)' ierr = ', ierr - CALL ABORT - ENDIF - -c----------------------------------------------------------------------- -c Lecture du tableau des parametres du run (pour la dynamique) -c----------------------------------------------------------------------- - - write(*,*) 'lecture tab_cntrl START_ARCHIVE' -c - ierr = NF_INQ_VARID (nid, "controle", nvarid) -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) -#endif -c - write(*,*) 'Impression de tab_cntrl' - do i=1,200 - write(*,*) i,tab_cntrl(i) - enddo - -c----------------------------------------------------------------------- -c Initialisation des constantes -c----------------------------------------------------------------------- - - imold = tab_cntrl(1) - jmold = tab_cntrl(2) - lmold = tab_cntrl(3) - day_ref = tab_cntrl(4) - annee_ref = tab_cntrl(5) - rad = tab_cntrl(6) - omeg = tab_cntrl(7) - g = tab_cntrl(8) - cpp = tab_cntrl(9) - kappa = tab_cntrl(10) - daysec = tab_cntrl(11) - dtvr = tab_cntrl(12) - etot0 = tab_cntrl(13) - ptot0 = tab_cntrl(14) - ztot0 = tab_cntrl(15) - stot0 = tab_cntrl(16) - ang0 = tab_cntrl(17) - pa = tab_cntrl(18) - preff = tab_cntrl(19) -c - clon = tab_cntrl(20) - clat = tab_cntrl(21) - grossismx = tab_cntrl(22) - grossismy = tab_cntrl(23) - - IF ( tab_cntrl(24).EQ.1. ) THEN - fxyhypb = . TRUE . - dzoomx = tab_cntrl(25) - dzoomy = tab_cntrl(26) - taux = tab_cntrl(28) - tauy = tab_cntrl(29) - ELSE - fxyhypb = . FALSE . - ysinus = . FALSE . - IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE. - ENDIF - - ptotalold = tab_cntrl(2*length) - - write(*,*) "Old dimensions:" - write(*,*) "longitude: ",imold - write(*,*) "latitude: ",jmold - write(*,*) "altitude: ",lmold - write(*,*) - - ip1jmp1old = (imold+1-1/imold)*(jmold+1-1/jmold) - -c======================================================================= -c CHANGEMENT DE CONSTANTES CONTENUES DANS tab_cntrl -c======================================================================= -c changement de la resolution dans le fichier de controle - tab_cntrl(1) = REAL(iim) - tab_cntrl(2) = REAL(jjm) - tab_cntrl(3) = REAL(llm) - -c-------------------------------- -c We use a specific run.def to read new parameters that need to be changed -c-------------------------------- - -c Changement de cpp: - CALL getin('cpp',cpp) - kappa = (8.314511/43.44E-3)/cpp - tab_cntrl(9) = cpp - tab_cntrl(10) = kappa - -c CHANGING THE ZOOM PARAMETERS TO CHANGE THE GRID - CALL getin('clon',clon) - CALL getin('clat',clat) - tab_cntrl(20) = clon - tab_cntrl(21) = clat - CALL getin('grossismx',grossismx) - CALL getin('grossismy',grossismy) - tab_cntrl(22) = grossismx - tab_cntrl(23) = grossismy - CALL getin('fxyhypb',fxyhypb) - IF ( fxyhypb ) THEN - CALL getin('dzoomx',dzoomx) - CALL getin('dzoomy',dzoomy) - tab_cntrl(25) = dzoomx - tab_cntrl(26) = dzoomy - CALL getin('taux',taux) - CALL getin('tauy',tauy) - tab_cntrl(28) = taux - tab_cntrl(29) = tauy - ELSE - CALL getin('ysinus',ysinus) - tab_cntrl(27) = ysinus - ENDIF - -c changes must be done BEFORE these lines... - DO l=1,length - tab_cntrl_dyn(l)= tab_cntrl(l) - tab_cntrl_fi(l) = tab_cntrl(l+length) - ENDDO - -c----------------------------------------------------------------------- -c Autres initialisations -c----------------------------------------------------------------------- - - CALL iniconst - CALL inigeom - call inifilr - -c----------------------------------------------------------------------- -c Allocations des anciennes variables -c----------------------------------------------------------------------- - - allocate(rlonuold(imold+1), rlatvold(jmold)) - allocate(rlonvold(imold+1), rlatuold(jmold+1)) - allocate(nivsigsold(lmold+1),nivsigold(lmold)) - allocate(apold(lmold),bpold(lmold)) - allocate(presnivsold(lmold)) - allocate(uold(imold+1,jmold+1,lmold)) - allocate(vold(imold+1,jmold+1,lmold)) - allocate(told(imold+1,jmold+1,lmold)) - allocate(qold(imold+1,jmold+1,lmold,nqtot)) - allocate(psold(imold+1,jmold+1)) - allocate(phisold(imold+1,jmold+1)) - allocate(tsurfold(imold+1,jmold+1)) - allocate(tsoilold(imold+1,jmold+1,nsoilmx)) - allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1)) - allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1)) - allocate(dlwold(imold+1,jmold+1)) - allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1)) - allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1)) - allocate(ztheold(imold+1,jmold+1),zpicold(imold+1,jmold+1)) - allocate(zvalold(imold+1,jmold+1)) - - allocate(var (imold+1,jmold+1,llm)) - allocate(varp1 (imold+1,jmold+1,llm+1)) - - print*,"Initialisations OK" - -c----------------------------------------------------------------------- -c Lecture des longitudes et latitudes -c----------------------------------------------------------------------- -c - ierr = NF_INQ_VARID (nid, "rlonv", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "rlatu", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "rlonu", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "rlatv", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - CALL abort - ENDIF -c - - print*,"Lecture lon/lat OK" - -c----------------------------------------------------------------------- -c Lecture des niveaux verticaux -c----------------------------------------------------------------------- -c - ierr = NF_INQ_VARID (nid, "nivsigs", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigsold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigsold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "nivsig", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "ap", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ELSE -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, apold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - ENDIF - ENDIF -c - ierr = NF_INQ_VARID (nid, "bp", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, bpold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - CALL abort - END IF - - ierr = NF_INQ_VARID (nid, "presnivs", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, presnivsold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, presnivsold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " - CALL abort - ENDIF - -c----------------------------------------------------------------------- -c Lecture geopotentiel au sol -c----------------------------------------------------------------------- -c - ierr = NF_INQ_VARID (nid, "phisinit", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Le champ est absent de start.nc" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, phisold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "new_grid: Lecture echouee pour " - CALL abort - ENDIF - - print*,"Lecture niveaux et geopot OK" - -c----------------------------------------------------------------------- -c Lecture des champs 2D () -c----------------------------------------------------------------------- - - start=(/1,1,1,0/) - counter=(/imold+1,jmold+1,1,0/) - - ierr = NF_INQ_VARID (nid, "ps", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,psold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,psold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "tsurf", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,tsurfold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,tsurfold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - do isoil=1,nsoilmx - write(str2,'(i2.2)') isoil -c - ierr = NF_INQ_VARID (nid, "Tsoil"//str2, nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: - . Le champ <","Tsoil"//str2,"> est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter, - . tsoilold(1,1,isoil)) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter, - . tsoilold(1,1,isoil)) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: - . Lecture echouee pour <","Tsoil"//str2,">" - CALL abort - ENDIF - end do -c - ierr = NF_INQ_VARID (nid, "albe", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,albeold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,albeold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "radsol", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,radsolold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,radsolold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "sollw", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "solsw", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,solswold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,solswold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "dlw", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,dlwold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,dlwold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zmea", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - zmeaold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zmeaold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zmeaold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zstd", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - zstdold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zstdold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zstdold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zsig", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - zsigold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zsigold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zsigold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zgam", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - zgamold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zgamold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zgamold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zthe", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - ztheold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,ztheold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,ztheold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zpic", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - zpicold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zpicold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zpicold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid, "zval", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - PRINT*, " Il est donc initialise a zero" - zvalold=0. - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zvalold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zvalold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - - print*,"Lecture champs 2D OK" - -c----------------------------------------------------------------------- -c Lecture des champs 3D () -c----------------------------------------------------------------------- - - start=(/1,1,1,1/) - counter=(/imold+1,jmold+1,lmold,1/) -c - ierr = NF_INQ_VARID (nid,"temp", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, counter, told) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid, start, counter, told) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid,"u", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,uold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,uold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - ierr = NF_INQ_VARID (nid,"v", nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,vold) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,vold) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "lect_start_archive: Lecture echouee pour " - CALL abort - ENDIF -c - -c TNAME: IL EST LU A PARTIR DE traceur.def (mettre l'ancien si -c changement du nombre de traceurs) - - IF(iflag_trac.eq.1) THEN - DO iq=1,nqtot - ierr = NF_INQ_VARID (nid, tname(iq), nvarid) - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent" - PRINT*, " Il est donc initialise a zero" - qold(:,:,:,iq)=0. - ELSE -#ifdef NC_DOUBLE - ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,counter,qold(1,1,1,iq)) -#else - ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,qold(1,1,1,iq)) -#endif - IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour "//tname(iq) - CALL abort - ENDIF - ENDIF - ENDDO - ENDIF - - - print*,"Lecture champs 3D OK" - -c======================================================================= -c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables -c======================================================================= -c Interpolation horizontale puis passage dans la grille physique pour -c les variables physique -c Interpolation verticale puis horizontale pour chaque variable 3D -c======================================================================= - -c----------------------------------------------------------------------- -c Variable 2d : -c----------------------------------------------------------------------- - -c Relief -! topoflag = F: we keep the existing topography -! = T: we read it from the Relief.nc file -! topoflag need to be in the specific run.def for newstart - topoflag = . FALSE . - CALL getin('topoflag',topoflag) - - IF ( topoflag ) then - print*,"Topography (phis) read in file Relief.nc" - print*,"For Venus, map directly inverted in Relief.nc" - CALL startget_phys2d('surfgeo',iip1,jjp1,rlonv,rlatu,phis,0.0, - . jjm ,rlonu,rlatv,.true.) -c needed ? - phis(iip1,:) = phis(1,:) - - CALL startget_phys1d('zmea',iip1,jjp1,rlonv,rlatu,ngridmx,zmea, - . 0.0,jjm,rlonu,rlatv,.true.) - CALL startget_phys1d('zstd',iip1,jjp1,rlonv,rlatu,ngridmx,zstd, - . 0.0,jjm,rlonu,rlatv,.true.) - CALL startget_phys1d('zsig',iip1,jjp1,rlonv,rlatu,ngridmx,zsig, - . 0.0,jjm,rlonu,rlatv,.true.) - CALL startget_phys1d('zgam',iip1,jjp1,rlonv,rlatu,ngridmx,zgam, - . 0.0,jjm,rlonu,rlatv,.true.) - CALL startget_phys1d('zthe',iip1,jjp1,rlonv,rlatu,ngridmx,zthe, - . 0.0,jjm,rlonu,rlatv,.true.) - CALL startget_phys1d('zpic',iip1,jjp1,rlonv,rlatu,ngridmx,zpic, - . 0.0,jjm,rlonu,rlatv,.true.) - CALL startget_phys1d('zval',iip1,jjp1,rlonv,rlatu,ngridmx,zval, - . 0.0,jjm,rlonu,rlatv,.true.) - - ELSE - print*,'Using existing topography' - call interp_horiz (phisold,phis,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - - call interp_horiz (zmeaold,zmeaS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,zmeaS,zmea) - call interp_horiz (zstdold,zstdS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,zstdS,zstd) - call interp_horiz (zsigold,zsigS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,zsigS,zsig) - call interp_horiz (zgamold,zgamS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,zgamS,zgam) - call interp_horiz (ztheold,ztheS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,ztheS,zthe) - call interp_horiz (zpicold,zpicS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,zpicS,zpic) - call interp_horiz (zvalold,zvalS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,zvalS,zval) - ENDIF - - print*,"New surface geopotential OK" - -c Lat et lon pour physique - do i=1,iip1 - rlatS(i,:)=rlatu(:)*180./pi - enddo - call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlatS,rlat) - - do j=2,jjm - rlonS(:,j)=rlonv(:)*180./pi - enddo - rlonS(:,1)=0. - rlonS(:,jjp1)=0. - call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlonS,rlon) - -c Temperature de surface - call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf) -c write(44,*) 'tsurf', tsurf - -c Temperature du sous-sol - call interp_horiz(tsoilold,tsoilS, - & imold,jmold,iim,jjm,nsoilmx, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil) -c write(45,*) 'tsoil',tsoil - -! CHANGING ALBEDO: may be done through run.def - albedoflag = . FALSE . - CALL getin('albedoflag',albedoflag) - - IF ( albedoflag ) then - print*,"Albedo is reinitialized to the albedo value in run.def" - print*,"We may want to consider a map later on..." - albedo=0.1 - CALL getin('albedo',albedo) - albe=albedo - ELSE - call interp_horiz (albeold,albeS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,albeS,albe) - ENDIF - - call interp_horiz (radsolold,radsolS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,radsolS,radsol) - - call interp_horiz (sollwold,sollwS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwS,sollw) - - call interp_horiz (solswold,solswS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw) - - call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw) - - print*,"Nouvelles var physiques OK" - -c----------------------------------------------------------------------- -c Traitement special de la pression au sol : -c----------------------------------------------------------------------- - -c Extrapolation la pression dans la nouvelle grille - call interp_horiz(psold,ps,imold,jmold,iim,jjm,1, - & rlonuold,rlatvold,rlonu,rlatv) - -c On assure la conservation de la masse de l'atmosphere -c-------------------------------------------------------------- - - ptotal = 0. - DO j=1,jjp1 - DO i=1,iim - ptotal=ptotal+ps(i,j)*aire(i,j)/g - END DO - END DO - - write(*,*) - write(*,*)'Ancienne grille: masse de l atm :',ptotalold - write(*,*)'Nouvelle grille: masse de l atm :',ptotal - write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold - write(*,*) - - - DO j=1,jjp1 - DO i=1,iip1 - ps(i,j)=ps(i,j) * ptotalold/ptotal - END DO - END DO - -c la pression de surface et les temperatures ne sont pas reequilibrees en fonction -c de la nouvelle topographie... -c Si l'ajustement inevitable du debut pose des problemes, voir le newstart martien. - - print*,"Nouvelle ps OK" - print*,"If changes done on topography, beware !" - print*,"Some time may be needed for adjustments at the beginning" - print*,"so if unstable, relax the timestep and/or dissipation." - -c----------------------------------------------------------------------- -c Variable 3d : -c----------------------------------------------------------------------- - - CALL pression(ip1jmp1, ap, bp, ps, p3d) - if (disvert_type==1) then - CALL exner_hyb( ip1jmp1, ps, p3d,alpha,beta,pks, pk, pkf ) - else ! we assume that we are in the disvert_type==2 case - CALL exner_milieu( ip1jmp1, ps, p3d, beta, pks, pk, pkf ) - endif - -c temperatures atmospheriques - -c ATTENTION: peut servir, mais bon... -c modif: profil uniforme -c do l=1,lmold -c do j=1,jmold+1 -c do i=1,imold+1 -c told(i,j,l)=told(1,jmold/2,l) -c enddo -c enddo -c enddo - - write (*,*) 'told ', told (1,jmold+1,1) ! INFO - call interp_vert - & (told,var,lmold,llm,apold,bpold,ap,bp, - & psold,(imold+1)*(jmold+1)) - write (*,*) 'var ', var (1,jmold+1,1) ! INFO - call interp_horiz(var,t,imold,jmold,iim,jjm,llm, - & rlonuold,rlatvold,rlonu,rlatv) - write (*,*) 'T ', T(1,jjp1,1) ! INFO - -c passage grille physique pour restartphy.nc - call gr_dyn_fi (llm,iip1,jjp1,ngridmx,T,t_fi) - -! ADAPTATION GCM POUR CP(T) -c passage en temperature potentielle - call t2tpot(ip1jmp1*llm,T,teta,pk) -c on assure la periodicite - teta(iip1,:,:) = teta(1,:,:) - -c calcul des champ de vent; passage en vent covariant - write (*,*) 'uold ', uold (1,2,1) ! INFO - call interp_vert - & (uold,var,lmold,llm,apold,bpold,ap,bp, - & psold,(imold+1)*(jmold+1)) - write (*,*) 'var ', var (1,2,1) ! INFO - call interp_horiz(var,us,imold,jmold,iim,jjm,llm, - & rlonuold,rlatvold,rlonu,rlatv) - write (*,*) 'us ', us (1,2,1) ! INFO - - call interp_vert - & (vold,var,lmold,llm, - & apold,bpold,ap,bp,psold,(imold+1)*(jmold+1)) - call interp_horiz(var,vs,imold,jmold,iim,jjm,llm, - & rlonuold,rlatvold,rlonu,rlatv) - call scal_wind(us,vs,unat,vnat) - write (*,*) 'unat ', unat (1,2,1) ! INFO - do l=1,llm - do j = 1, jjp1 - do i=1,iip1 - ucov( i,j,l ) = unat( i,j,l ) * cu(i,j) -c ucov( i,j,l ) = 0 - end do - end do - end do - write (*,*) 'ucov ', ucov (1,2,1) ! INFO -c write(48,*) 'ucov',ucov - do l=1,llm - do j = 1, jjm - do i=1,iim - vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j) -c vcov( i,j,l ) = 0 - end do - vcov( iip1,j,l ) = vcov( 1,j,l ) - end do - end do -c write(49,*) 'ucov',vcov - -c masse: on la recalcule (ps a été ajustée pour conserver la masse totale) - call massdair(p3d,masse) - -c traceurs 3D - do iq = 1, nqtot - call interp_vert(qold(1,1,1,iq),var,lmold,llm, - & apold,bpold,ap,bp,psold,(imold+1)*(jmold+1)) - call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm, - & rlonuold,rlatvold,rlonu,rlatv) - enddo - - print*,"Nouvelles var dynamiques OK" - -c======================================================================= -c Ecriture des restart.nc et restartphy.nc : -c======================================================================= - - call writerestart('restart.nc',tab_cntrl_dyn, - . phis,vcov,ucov,teta,masse,q,ps) - - print*,"restart OK" - - call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm, - . rlat,rlon,tsurf,tsoil,albe,solsw, sollw,dlw, - . radsol, - . zmea,zstd,zsig,zgam,zthe,zpic,zval, - . t_fi) - - print*,"restartphy OK" - - end Index: unk/LMDZ.COMMON/libf/dyn3d/start2archive-VT.F =================================================================== --- /trunk/LMDZ.COMMON/libf/dyn3d/start2archive-VT.F (revision 848) +++ (revision ) @@ -1,311 +1,0 @@ -c======================================================================= - PROGRAM start2archive -c======================================================================= -c -c -c Date: 01/1997 -c ---- -c -c Version Venus: 09/2007 -c Titan: 02/2009 -c -c Objet: Passage des fichiers netcdf d etat initial "start" et -c ----- "startphy" a un fichier netcdf unique "start_archive" -c -c======================================================================= - - USE filtreg_mod - USE infotrac - USE control_mod - - implicit none - -#include "dimensions.h" -#include "paramet.h" -#include "comconst.h" -#include "comdissnew.h" -#include "comvert.h" -#include "comgeom.h" -#include "logic.h" -#include "temps.h" -#include "ener.h" -#include "description.h" -#include "dimsoil.h" -#include "netcdf.inc" - -c----------------------------------------------------------------------- -c Declarations -c----------------------------------------------------------------------- - -c variables dynamiques du GCM -c ----------------------------- - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants - REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes - REAL pks(ip1jmp1) ! exner (f pour filtre) - REAL pk(ip1jmp1,llm) - REAL pkf(ip1jmp1,llm) - REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) - REAL phis(ip1jmp1) ! geopotentiel au sol - REAL masse(ip1jmp1,llm) ! masse de l'atmosphere - REAL ps(ip1jmp1) ! pression au sol - REAL p3d(iip1, jjp1, llm+1) ! pression aux interfaces - -c Variable Physiques (grille physique) -c ------------------------------------ - integer ngridmx,nlayermx - parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm)) - parameter (nlayermx=llm) - - real rlat(ngridmx),rlon(ngridmx) - REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx) - REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx) - real solsw(ngridmx),dlw(ngridmx) - REAL zmea(ngridmx), zstd(ngridmx) - REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx) - REAL zpic(ngridmx), zval(ngridmx) - - INTEGER start,length - PARAMETER (length = 100) - REAL tab_cntrl_fi(length) ! tableau des parametres de startfi - REAL tab_cntrl_dyn(length) ! tableau des parametres de start - INTEGER*4 day_ini_fi - -c Variable naturelle / grille scalaire -c ------------------------------------ - REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm) - REAL tsurfS(ip1jmp1),tsoilS(ip1jmp1,nsoilmx) - real rlatS(ip1jmp1),rlonS(ip1jmp1) - real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1) - real solswS(ip1jmp1),dlwS(ip1jmp1) - real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1) - real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1) - real zvalS(ip1jmp1) - -c Variables intermediaires : vent naturel, mais pas coord scalaire -c---------------------------------------------------------------- - REAL vn(ip1jm,llm),un(ip1jmp1,llm) - -c Autres variables -c ----------------- - REAL ptotal - - CHARACTER*2 str2 - - INTEGER ij, l,i,j,isoil,iq - character*80 fichnom - integer :: ierr - -c Netcdf -c------- - integer varid,dimid - INTEGER nid - -c----------------------------------------------------------------------- -c Initialisations -c----------------------------------------------------------------------- - -c VENUS/TITAN - - iflag_trac = 1 -c----------------------------------------------------------------------- -c Initialisation des traceurs -c --------------------------- -c Choix du nombre de traceurs et du schema pour l advection -c dans fichier traceur.def, par default ou via INCA - call infotrac_init - -c Allocation de la tableau q : champs advectes - allocate(q(ip1jmp1,llm,nqtot)) - -c======================================================================= -c Lecture des donnees -c======================================================================= - - fichnom = 'start.nc' - CALL readstart(fichnom,nqtot,vcov,ucov,teta,q,masse, - . ps,phis,tab_cntrl_dyn) - - fichnom = 'startphy.nc' - CALL readstartphy(fichnom, - . rlat,rlon,tsurf,tsoil, - . albe, solsw, sollw, - . dlw,radsol, - . zmea,zstd,zsig,zgam,zthe,zpic,zval, - . tab_cntrl_fi) - -c----------------------------------------------------------------------- -c Initialisations -c----------------------------------------------------------------------- - - call iniconst - call inigeom - call inifilr - CALL pression(ip1jmp1, ap, bp, ps, p3d) - if (disvert_type==1) then - CALL exner_hyb( ip1jmp1, ps, p3d,alpha,beta,pks, pk, pkf ) - else ! we assume that we are in the disvert_type==2 case - CALL exner_milieu( ip1jmp1, ps, p3d, beta, pks, pk, pkf ) - endif - -c======================================================================= -c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire -c======================================================================= -c Les variables modeles dependent de la resolution. Il faut donc -c eliminer les facteurs responsables de cette dependance -c (pour utiliser newstart) -c======================================================================= - -c----------------------------------------------------------------------- -c Vent (depend de la resolution horizontale) -c----------------------------------------------------------------------- -c -c ucov --> un et vcov --> vn -c un --> us et vn --> vs -c -c----------------------------------------------------------------------- - - call covnat(llm,ucov, vcov, un, vn) - call wind_scal(un,vn,us,vs) - -c----------------------------------------------------------------------- -c Temperature (depend de la resolution verticale => de "sigma.def") -c----------------------------------------------------------------------- -c -c h --> T -c -c----------------------------------------------------------------------- -! ADAPTATION GCM POUR CP(T) - - call tpot2t(ip1jmp1*llm,teta,T,pk) - -c----------------------------------------------------------------------- -c Variable physique -c----------------------------------------------------------------------- -c -c tsurf --> tsurfS -c et autres... -c -c----------------------------------------------------------------------- - - call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS) - call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,rlat,rlatS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,rlon,rlonS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,albe,albeS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,radsol,radsolS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,sollw,sollwS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,solsw,solswS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,dlw,dlwS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zmea,zmeaS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zstd,zstdS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zsig,zsigS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zgam,zgamS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zthe,ztheS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zpic,zpicS) - call gr_fi_dyn(1,ngridmx,iip1,jjp1,zval,zvalS) - -c======================================================================= -c Info pour controler -c======================================================================= - - ptotal = 0. - DO j=1,jjp1 - DO i=1,iim - ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g - ENDDO - ENDDO - write(*,*)'Ancienne grille : masse de l''atm :',ptotal - -c----------------------------------------------------------------------- -c Passage de "ptotal" par tab_cntrl_fi -c----------------------------------------------------------------------- - - tab_cntrl_fi(length) = ptotal - -c======================================================================= -c Ecriture dans le fichier "start_archive" -c======================================================================= - -c----------------------------------------------------------------------- -c Ouverture de "start_archive" -c----------------------------------------------------------------------- - - ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid) - -c----------------------------------------------------------------------- -c si "start_archive" n'existe pas: -c 1_ ouverture -c 2_ creation de l'entete dynamique ("ini_archive") -c----------------------------------------------------------------------- -c ini_archive: -c On met dans l'entete le tab_cntrl_dyn (1 a length) -c On y ajoute les valeurs du tab_cntrl_fi (length+1 a 2*length) -c----------------------------------------------------------------------- - - if (ierr.ne.NF_NOERR) then - write(*,*)'OK, Could not open file "start_archive.nc"' - write(*,*)'So let s create a new "start_archive"' - ierr = NF_CREATE('start_archive.nc', NF_CLOBBER, nid) - call ini_archive(nid,phis,tab_cntrl_dyn,tab_cntrl_fi) - else - write(*,*)'Attention, start_archive.nc existe déjà...' - call abort - endif - -c----------------------------------------------------------------------- -c Ecriture des champs -c----------------------------------------------------------------------- - - call write_archive(nid,'u','Vent zonal','m.s-1',3,us) - call write_archive(nid,'v','Vent merid','m.s-1',3,vs) - call write_archive(nid,'temp','temperature','K',3,T) -c----------------------------------------------------------------------- -c Ecriture du champs q ( q[1,nqtot] ) -c----------------------------------------------------------------------- - do iq=1,nqtot - write(str2,'(i2.2)') iq - call write_archive(nid,tname(iq),'tracer','kg/kg', - . 3,q(1,1,iq)) - end do -c----------------------------------------------------------------------- - call write_archive(nid,'masse','Masse','kg',3,masse) - call write_archive(nid,'ps','Psurf','Pa',2,ps) - call write_archive(nid,'tsurf','surf T','K',2,tsurfS) -c----------------------------------------------------------------------- -c Ecriture du champs tsoil ( Tsoil[1,nsoilmx] ) -c----------------------------------------------------------------------- -c "tsoil" Temperature au sol definie dans nsoilmx couches dans le sol -c Les nsoilmx couches sont lues comme nsoilmx champs -c nommees Tsoil[1,nsoilmx] - do isoil=1,nsoilmx - write(str2,'(i2.2)') isoil - call write_archive(nid,'Tsoil'//str2,'Ground Temperature ', - . 'K',2,tsoilS(1,isoil)) - enddo -c----------------------------------------------------------------------- - call write_archive(nid,'rlat','Latitude','rad',2,rlatS) - call write_archive(nid,'rlon','Longitude','rad',2,rlonS) - call write_archive(nid,'albe','Albedo','',2,albeS) - call write_archive(nid,'radsol', - . 'Net flux at surface','W m-2',2,radsolS) - call write_archive(nid,'sollw', - . 'LW flux at surface','W m-2',2,sollwS) - call write_archive(nid,'solsw', - . 'SW flux at surface','W m-2',2,solswS) - call write_archive(nid,'dlw','LW derive','?',2,dlwS) - call write_archive(nid,'zmea','param oro sous-maille','m',2,zmeaS) - call write_archive(nid,'zstd','param oro sous-maille','m',2,zstdS) - call write_archive(nid,'zsig','param oro sous-maille','m',2,zsigS) - call write_archive(nid,'zgam','param oro sous-maille','m',2,zgamS) - call write_archive(nid,'zthe','param oro sous-maille','m',2,ztheS) - call write_archive(nid,'zpic','param oro sous-maille','m',2,zpicS) - call write_archive(nid,'zval','param oro sous-maille','m',2,zvalS) - - ierr=NF_CLOSE(nid) - -c----------------------------------------------------------------------- -c Fin -c----------------------------------------------------------------------- - - end Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F =================================================================== --- /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F (revision 848) +++ /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F (revision 849) @@ -512,5 +512,5 @@ IF (is_sequential.and.(planet_type=="earth")) THEN -#ifdef CPP_PHYS +#ifdef CPP_EARTH ! PVtheta calls tetalevel, which is in the physics cIM calcul PV a teta=350, 380, 405K Index: /trunk/LMDZ.COMMON/makelmdz =================================================================== --- /trunk/LMDZ.COMMON/makelmdz (revision 848) +++ /trunk/LMDZ.COMMON/makelmdz (revision 849) @@ -544,8 +544,32 @@ ################################################################# -source_code=${code}.F -if [[ -r $LMDGCM/libf/dyn${dimc}d${FLAG_PARA}/${code}.F90 ]] -then - source_code=${code}.F90 +# find code suffix and directory where code is located +if [[ -r $LMDGCM/libf/dyn${dimc}d${FLAG_PARA}/${code}.F ]] +then + source_code=${code}.F + code_dir=dyn${dimc}d${FLAG_PARA} +else + if [[ -r $LMDGCM/libf/dyn${dimc}d${FLAG_PARA}/${code}.F90 ]] + then + source_code=${code}.F90 + code_dir=dyn${dimc}d${FLAG_PARA} + else + if [[ -r $LMDGCM/libf/phy${physique}/${code}.F ]] + then + source_code=${code}.F + code_dir=phy${physique} + else + # last possibility: + if [[ -r $LMDGCM/libf/phy${physique}/${code}.F ]] + then + source_code=${code}.F90 + code_dir=phy${physique} + else + echo "Error: cannot find ${code}.F[90]" + echo " neither in dyn${dimc}d${FLAG_PARA} nor in phy${physique}" + exit + fi + fi + fi fi @@ -625,4 +649,5 @@ MOD_SUFFIX="mod" \ AR=$arcommand \ +DIRMAIN=$code_dir \ SOURCE=$source_code \ PROG=$code @@ -654,4 +679,5 @@ MOD_SUFFIX="mod" \ AR=$arcommand \ +DIRMAIN=$code_dir \ SOURCE=$source_code \ PROG=$code Index: /trunk/LMDZ.TITAN/libf/phytitan/physiq.F =================================================================== --- /trunk/LMDZ.TITAN/libf/phytitan/physiq.F (revision 848) +++ /trunk/LMDZ.TITAN/libf/phytitan/physiq.F (revision 849) @@ -56,5 +56,5 @@ c====================================================================== USE ioipsl - USE histcom +! USE histcom ! not needed; histcom is included in ioipsl USE infotrac USE control_mod @@ -279,4 +279,5 @@ REAL zday REAL zls,zlsdeg + REAL zlsm1 ! for write_hist* (NB: it is not properly initialized!!) c INTEGER i, k, iq, ig, j, ll, l Index: /trunk/LMDZ.VENUS/libf/phyvenus/newstart.F =================================================================== --- /trunk/LMDZ.VENUS/libf/phyvenus/newstart.F (revision 849) +++ /trunk/LMDZ.VENUS/libf/phyvenus/newstart.F (revision 849) @@ -0,0 +1,1104 @@ +C====================================================================== + PROGRAM newstart +c======================================================================= +c +c +c Auteur: S. Lebonnois, +c a partir des newstart/start_archive/lect_start_archive martiens +c +c Derniere modif : 02/09 (ecriture des q*) +c 01/12 (inclusion dans svn dyn3d) +c +c Objet: Modify the grid for the initial state (LMD GCM VENUS/TITAN) +c ----- (from file NetCDF start_archive.nc) +c +c +c======================================================================= + + use IOIPSL + USE filtreg_mod + USE startvar + USE control_mod + USE infotrac + + implicit none + +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comdissnew.h" +#include "comvert.h" +#include "comgeom2.h" +#include "logic.h" +#include "temps.h" +#include "ener.h" +#include "description.h" +#include "serre.h" +#include "dimsoil.h" +#include "netcdf.inc" + +c----------------------------------------------------------------------- +c Declarations +c----------------------------------------------------------------------- + +c Variables pour fichier "ini" +c------------------------------------ + INTEGER imold,jmold,lmold,nqold,ip1jmp1old + INTEGER length + parameter (length = 100) + real tab_cntrl(2*length) + INTEGER isoil,iq,iqmax + CHARACTER*2 str2 + +c Variable histoire +c------------------ + REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants + REAL teta(iip1,jjp1,llm),ps(iip1,jjp1) + REAL phis(iip1,jjp1) ! geopotentiel au sol + REAL masse(ip1jmp1,llm) ! masse de l'atmosphere + REAL, ALLOCATABLE, DIMENSION(:,:,:,:):: q! champs advectes + REAL tab_cntrl_dyn(length) ! tableau des parametres de start + +c variable physique +c------------------ + integer ngridmx + parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm)) + REAL tab_cntrl_fi(length) ! tableau des parametres de startfi + real rlat(ngridmx),rlon(ngridmx) + REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx) + REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx) + real solsw(ngridmx),dlw(ngridmx) + REAL zmea(ngridmx), zstd(ngridmx) + REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx) + REAL zpic(ngridmx), zval(ngridmx) + real t_fi(ngridmx,llm) + +c Variable nouvelle grille naturelle au point scalaire +c------------------------------------------------------ + real us(iip1,jjp1,llm),vs(iip1,jjp1,llm) + REAL p3d(iip1,jjp1,llm+1) ! pression aux interfaces + REAL phisold_newgrid(iip1,jjp1) + REAL T(iip1,jjp1,llm) + real rlonS(iip1,jjp1),rlatS(iip1,jjp1) + real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx) + real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1) + real solswS(ip1jmp1),dlwS(ip1jmp1) + real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1) + real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1) + real zvalS(ip1jmp1) + + real ptotal + +c Var intermediaires : vent naturel, mais pas coord scalaire +c----------------------------------------------------------- + real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) + + REAL pks(iip1,jjp1) ! exner (f pour filtre) + REAL pk(iip1,jjp1,llm) + REAL pkf(iip1,jjp1,llm) + REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) + + +c Variable de l'ancienne grille +c--------------------------------------------------------- + + real, dimension(:), allocatable :: rlonuold, rlatvold + real, dimension(:), allocatable :: rlonvold, rlatuold + real, dimension(:), allocatable :: nivsigsold,nivsigold + real, dimension(:), allocatable :: apold,bpold + real, dimension(:), allocatable :: presnivsold + real, dimension(:,:,:), allocatable :: uold,vold,told + real, dimension(:,:,:,:), allocatable :: qold + real, dimension(:,:,:), allocatable :: tsoilold + real, dimension(:,:), allocatable :: psold,phisold + real, dimension(:,:), allocatable :: tsurfold + real, dimension(:,:), allocatable :: albeold,radsolold + real, dimension(:,:), allocatable :: sollwold,solswold + real, dimension(:,:), allocatable :: dlwold + real, dimension(:,:), allocatable :: zmeaold,zstdold,zsigold + real, dimension(:,:), allocatable :: zgamold,ztheold,zpicold + real, dimension(:,:), allocatable :: zvalold + + real ptotalold + +c Variable intermediaires iutilise pour l'extrapolation verticale +c---------------------------------------------------------------- + real, dimension(:,:,:), allocatable :: var,varp1 + +c divers local +c----------------- + + integer ierr,nid,nvarid + INTEGER ij, l,i,j + character*80 fichnom + integer, dimension(4) :: start,counter + REAL phisinverse(iip1,jjp1) ! geopotentiel au sol avant inversion + logical topoflag,albedoflag + real albedo + +c======================================================================= +c INITIALISATIONS DIVERSES +c======================================================================= + +c VENUS/TITAN + + iflag_trac = 1 +c----------------------------------------------------------------------- +c Initialisation des traceurs +c --------------------------- +c Choix du nombre de traceurs et du schema pour l'advection +c dans fichier traceur.def, par default ou via INCA + call infotrac_init + +c Allocation de la tableau q : champs advectes + allocate(q(iip1,jjp1,llm,nqtot)) + +c----------------------------------------------------------------------- +c Ouverture du fichier a modifier (start_archive.nc) +c----------------------------------------------------------------------- + + write(*,*) 'Creation d un etat initial a partir de' + write(*,*) './start_archive.nc' + write(*,*) + fichnom = 'start_archive.nc' + ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) + IF (ierr.NE.NF_NOERR) THEN + write(6,*)' Pb d''ouverture du fichier ',fichnom + write(6,*)' ierr = ', ierr + CALL ABORT + ENDIF + +c----------------------------------------------------------------------- +c Lecture du tableau des parametres du run (pour la dynamique) +c----------------------------------------------------------------------- + + write(*,*) 'lecture tab_cntrl START_ARCHIVE' +c + ierr = NF_INQ_VARID (nid, "controle", nvarid) +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) +#endif +c + write(*,*) 'Impression de tab_cntrl' + do i=1,200 + write(*,*) i,tab_cntrl(i) + enddo + +c----------------------------------------------------------------------- +c Initialisation des constantes +c----------------------------------------------------------------------- + + imold = tab_cntrl(1) + jmold = tab_cntrl(2) + lmold = tab_cntrl(3) + day_ref = tab_cntrl(4) + annee_ref = tab_cntrl(5) + rad = tab_cntrl(6) + omeg = tab_cntrl(7) + g = tab_cntrl(8) + cpp = tab_cntrl(9) + kappa = tab_cntrl(10) + daysec = tab_cntrl(11) + dtvr = tab_cntrl(12) + etot0 = tab_cntrl(13) + ptot0 = tab_cntrl(14) + ztot0 = tab_cntrl(15) + stot0 = tab_cntrl(16) + ang0 = tab_cntrl(17) + pa = tab_cntrl(18) + preff = tab_cntrl(19) +c + clon = tab_cntrl(20) + clat = tab_cntrl(21) + grossismx = tab_cntrl(22) + grossismy = tab_cntrl(23) + + IF ( tab_cntrl(24).EQ.1. ) THEN + fxyhypb = . TRUE . + dzoomx = tab_cntrl(25) + dzoomy = tab_cntrl(26) + taux = tab_cntrl(28) + tauy = tab_cntrl(29) + ELSE + fxyhypb = . FALSE . + ysinus = . FALSE . + IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE. + ENDIF + + ptotalold = tab_cntrl(2*length) + + write(*,*) "Old dimensions:" + write(*,*) "longitude: ",imold + write(*,*) "latitude: ",jmold + write(*,*) "altitude: ",lmold + write(*,*) + + ip1jmp1old = (imold+1-1/imold)*(jmold+1-1/jmold) + +c======================================================================= +c CHANGEMENT DE CONSTANTES CONTENUES DANS tab_cntrl +c======================================================================= +c changement de la resolution dans le fichier de controle + tab_cntrl(1) = REAL(iim) + tab_cntrl(2) = REAL(jjm) + tab_cntrl(3) = REAL(llm) + +c-------------------------------- +c We use a specific run.def to read new parameters that need to be changed +c-------------------------------- + +c Changement de cpp: + CALL getin('cpp',cpp) + kappa = (8.314511/43.44E-3)/cpp + tab_cntrl(9) = cpp + tab_cntrl(10) = kappa + +c CHANGING THE ZOOM PARAMETERS TO CHANGE THE GRID + CALL getin('clon',clon) + CALL getin('clat',clat) + tab_cntrl(20) = clon + tab_cntrl(21) = clat + CALL getin('grossismx',grossismx) + CALL getin('grossismy',grossismy) + tab_cntrl(22) = grossismx + tab_cntrl(23) = grossismy + CALL getin('fxyhypb',fxyhypb) + IF ( fxyhypb ) THEN + CALL getin('dzoomx',dzoomx) + CALL getin('dzoomy',dzoomy) + tab_cntrl(25) = dzoomx + tab_cntrl(26) = dzoomy + CALL getin('taux',taux) + CALL getin('tauy',tauy) + tab_cntrl(28) = taux + tab_cntrl(29) = tauy + ELSE + CALL getin('ysinus',ysinus) + tab_cntrl(27) = ysinus + ENDIF + +c changes must be done BEFORE these lines... + DO l=1,length + tab_cntrl_dyn(l)= tab_cntrl(l) + tab_cntrl_fi(l) = tab_cntrl(l+length) + ENDDO + +c----------------------------------------------------------------------- +c Autres initialisations +c----------------------------------------------------------------------- + + CALL iniconst + CALL inigeom + call inifilr + +c----------------------------------------------------------------------- +c Allocations des anciennes variables +c----------------------------------------------------------------------- + + allocate(rlonuold(imold+1), rlatvold(jmold)) + allocate(rlonvold(imold+1), rlatuold(jmold+1)) + allocate(nivsigsold(lmold+1),nivsigold(lmold)) + allocate(apold(lmold),bpold(lmold)) + allocate(presnivsold(lmold)) + allocate(uold(imold+1,jmold+1,lmold)) + allocate(vold(imold+1,jmold+1,lmold)) + allocate(told(imold+1,jmold+1,lmold)) + allocate(qold(imold+1,jmold+1,lmold,nqtot)) + allocate(psold(imold+1,jmold+1)) + allocate(phisold(imold+1,jmold+1)) + allocate(tsurfold(imold+1,jmold+1)) + allocate(tsoilold(imold+1,jmold+1,nsoilmx)) + allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1)) + allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1)) + allocate(dlwold(imold+1,jmold+1)) + allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1)) + allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1)) + allocate(ztheold(imold+1,jmold+1),zpicold(imold+1,jmold+1)) + allocate(zvalold(imold+1,jmold+1)) + + allocate(var (imold+1,jmold+1,llm)) + allocate(varp1 (imold+1,jmold+1,llm+1)) + + print*,"Initialisations OK" + +c----------------------------------------------------------------------- +c Lecture des longitudes et latitudes +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "rlonv", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "rlatu", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "rlonu", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "rlatv", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + CALL abort + ENDIF +c + + print*,"Lecture lon/lat OK" + +c----------------------------------------------------------------------- +c Lecture des niveaux verticaux +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "nivsigs", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigsold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigsold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Lecture echouee pour " + CALL abort + ENDIF + + ierr = NF_INQ_VARID (nid, "nivsig", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, nivsigold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, nivsigold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Lecture echouee pour " + CALL abort + ENDIF + + ierr = NF_INQ_VARID (nid, "ap", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ELSE +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, apold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + ENDIF + ENDIF +c + ierr = NF_INQ_VARID (nid, "bp", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, bpold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + CALL abort + END IF + + ierr = NF_INQ_VARID (nid, "presnivs", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, presnivsold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, presnivsold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Lecture echouee pour " + CALL abort + ENDIF + +c----------------------------------------------------------------------- +c Lecture geopotentiel au sol +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "phisinit", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Le champ est absent de start.nc" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, phisold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "new_grid: Lecture echouee pour " + CALL abort + ENDIF + + print*,"Lecture niveaux et geopot OK" + +c----------------------------------------------------------------------- +c Lecture des champs 2D () +c----------------------------------------------------------------------- + + start=(/1,1,1,0/) + counter=(/imold+1,jmold+1,1,0/) + + ierr = NF_INQ_VARID (nid, "ps", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,psold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,psold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "tsurf", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,tsurfold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,tsurfold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + do isoil=1,nsoilmx + write(str2,'(i2.2)') isoil +c + ierr = NF_INQ_VARID (nid, "Tsoil"//str2, nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: + . Le champ <","Tsoil"//str2,"> est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter, + . tsoilold(1,1,isoil)) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter, + . tsoilold(1,1,isoil)) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: + . Lecture echouee pour <","Tsoil"//str2,">" + CALL abort + ENDIF + end do +c + ierr = NF_INQ_VARID (nid, "albe", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,albeold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,albeold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "radsol", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,radsolold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,radsolold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "sollw", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "solsw", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,solswold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,solswold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "dlw", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,dlwold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,dlwold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zmea", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + zmeaold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zmeaold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zmeaold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zstd", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + zstdold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zstdold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zstdold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zsig", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + zsigold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zsigold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zsigold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zgam", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + zgamold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zgamold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zgamold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zthe", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + ztheold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,ztheold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,ztheold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zpic", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + zpicold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zpicold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zpicold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "zval", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + PRINT*, " Il est donc initialise a zero" + zvalold=0. + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,zvalold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,zvalold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + + print*,"Lecture champs 2D OK" + +c----------------------------------------------------------------------- +c Lecture des champs 3D () +c----------------------------------------------------------------------- + + start=(/1,1,1,1/) + counter=(/imold+1,jmold+1,lmold,1/) +c + ierr = NF_INQ_VARID (nid,"temp", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, counter, told) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid, start, counter, told) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid,"u", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,uold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,uold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid,"v", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Le champ est absent" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,vold) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,vold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Lecture echouee pour " + CALL abort + ENDIF +c + +c TNAME: IL EST LU A PARTIR DE traceur.def (mettre l'ancien si +c changement du nombre de traceurs) + + IF(iflag_trac.eq.1) THEN + DO iq=1,nqtot + ierr = NF_INQ_VARID (nid, tname(iq), nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent" + PRINT*, " Il est donc initialise a zero" + qold(:,:,:,iq)=0. + ELSE +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,counter,qold(1,1,1,iq)) +#else + ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,qold(1,1,1,iq)) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "dynetat0: Lecture echouee pour "//tname(iq) + CALL abort + ENDIF + ENDIF + ENDDO + ENDIF + + + print*,"Lecture champs 3D OK" + +c======================================================================= +c INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables +c======================================================================= +c Interpolation horizontale puis passage dans la grille physique pour +c les variables physique +c Interpolation verticale puis horizontale pour chaque variable 3D +c======================================================================= + +c----------------------------------------------------------------------- +c Variable 2d : +c----------------------------------------------------------------------- + +c Relief +! topoflag = F: we keep the existing topography +! = T: we read it from the Relief.nc file +! topoflag need to be in the specific run.def for newstart + topoflag = . FALSE . + CALL getin('topoflag',topoflag) + + IF ( topoflag ) then + print*,"Topography (phis) read in file Relief.nc" + print*,"For Venus, map directly inverted in Relief.nc" + CALL startget_phys2d('surfgeo',iip1,jjp1,rlonv,rlatu,phis,0.0, + . jjm ,rlonu,rlatv,.true.) +c needed ? + phis(iip1,:) = phis(1,:) + + CALL startget_phys1d('zmea',iip1,jjp1,rlonv,rlatu,ngridmx,zmea, + . 0.0,jjm,rlonu,rlatv,.true.) + CALL startget_phys1d('zstd',iip1,jjp1,rlonv,rlatu,ngridmx,zstd, + . 0.0,jjm,rlonu,rlatv,.true.) + CALL startget_phys1d('zsig',iip1,jjp1,rlonv,rlatu,ngridmx,zsig, + . 0.0,jjm,rlonu,rlatv,.true.) + CALL startget_phys1d('zgam',iip1,jjp1,rlonv,rlatu,ngridmx,zgam, + . 0.0,jjm,rlonu,rlatv,.true.) + CALL startget_phys1d('zthe',iip1,jjp1,rlonv,rlatu,ngridmx,zthe, + . 0.0,jjm,rlonu,rlatv,.true.) + CALL startget_phys1d('zpic',iip1,jjp1,rlonv,rlatu,ngridmx,zpic, + . 0.0,jjm,rlonu,rlatv,.true.) + CALL startget_phys1d('zval',iip1,jjp1,rlonv,rlatu,ngridmx,zval, + . 0.0,jjm,rlonu,rlatv,.true.) + + ELSE + print*,'Using existing topography' + call interp_horiz (phisold,phis,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + + call interp_horiz (zmeaold,zmeaS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,zmeaS,zmea) + call interp_horiz (zstdold,zstdS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,zstdS,zstd) + call interp_horiz (zsigold,zsigS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,zsigS,zsig) + call interp_horiz (zgamold,zgamS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,zgamS,zgam) + call interp_horiz (ztheold,ztheS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,ztheS,zthe) + call interp_horiz (zpicold,zpicS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,zpicS,zpic) + call interp_horiz (zvalold,zvalS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,zvalS,zval) + ENDIF + + print*,"New surface geopotential OK" + +c Lat et lon pour physique + do i=1,iip1 + rlatS(i,:)=rlatu(:)*180./pi + enddo + call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlatS,rlat) + + do j=2,jjm + rlonS(:,j)=rlonv(:)*180./pi + enddo + rlonS(:,1)=0. + rlonS(:,jjp1)=0. + call gr_dyn_fi (1,iip1,jjp1,ngridmx,rlonS,rlon) + +c Temperature de surface + call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf) +c write(44,*) 'tsurf', tsurf + +c Temperature du sous-sol + call interp_horiz(tsoilold,tsoilS, + & imold,jmold,iim,jjm,nsoilmx, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil) +c write(45,*) 'tsoil',tsoil + +! CHANGING ALBEDO: may be done through run.def + albedoflag = . FALSE . + CALL getin('albedoflag',albedoflag) + + IF ( albedoflag ) then + print*,"Albedo is reinitialized to the albedo value in run.def" + print*,"We may want to consider a map later on..." + albedo=0.1 + CALL getin('albedo',albedo) + albe=albedo + ELSE + call interp_horiz (albeold,albeS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,albeS,albe) + ENDIF + + call interp_horiz (radsolold,radsolS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,radsolS,radsol) + + call interp_horiz (sollwold,sollwS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwS,sollw) + + call interp_horiz (solswold,solswS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw) + + call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw) + + print*,"Nouvelles var physiques OK" + +c----------------------------------------------------------------------- +c Traitement special de la pression au sol : +c----------------------------------------------------------------------- + +c Extrapolation la pression dans la nouvelle grille + call interp_horiz(psold,ps,imold,jmold,iim,jjm,1, + & rlonuold,rlatvold,rlonu,rlatv) + +c On assure la conservation de la masse de l'atmosphere +c-------------------------------------------------------------- + + ptotal = 0. + DO j=1,jjp1 + DO i=1,iim + ptotal=ptotal+ps(i,j)*aire(i,j)/g + END DO + END DO + + write(*,*) + write(*,*)'Ancienne grille: masse de l atm :',ptotalold + write(*,*)'Nouvelle grille: masse de l atm :',ptotal + write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold + write(*,*) + + + DO j=1,jjp1 + DO i=1,iip1 + ps(i,j)=ps(i,j) * ptotalold/ptotal + END DO + END DO + +c la pression de surface et les temperatures ne sont pas reequilibrees en fonction +c de la nouvelle topographie... +c Si l'ajustement inevitable du debut pose des problemes, voir le newstart martien. + + print*,"Nouvelle ps OK" + print*,"If changes done on topography, beware !" + print*,"Some time may be needed for adjustments at the beginning" + print*,"so if unstable, relax the timestep and/or dissipation." + +c----------------------------------------------------------------------- +c Variable 3d : +c----------------------------------------------------------------------- + + CALL pression(ip1jmp1, ap, bp, ps, p3d) + if (disvert_type==1) then + CALL exner_hyb( ip1jmp1, ps, p3d,alpha,beta,pks, pk, pkf ) + else ! we assume that we are in the disvert_type==2 case + CALL exner_milieu( ip1jmp1, ps, p3d, beta, pks, pk, pkf ) + endif + +c temperatures atmospheriques + +c ATTENTION: peut servir, mais bon... +c modif: profil uniforme +c do l=1,lmold +c do j=1,jmold+1 +c do i=1,imold+1 +c told(i,j,l)=told(1,jmold/2,l) +c enddo +c enddo +c enddo + + write (*,*) 'told ', told (1,jmold+1,1) ! INFO + call interp_vert + & (told,var,lmold,llm,apold,bpold,ap,bp, + & psold,(imold+1)*(jmold+1)) + write (*,*) 'var ', var (1,jmold+1,1) ! INFO + call interp_horiz(var,t,imold,jmold,iim,jjm,llm, + & rlonuold,rlatvold,rlonu,rlatv) + write (*,*) 'T ', T(1,jjp1,1) ! INFO + +c passage grille physique pour restartphy.nc + call gr_dyn_fi (llm,iip1,jjp1,ngridmx,T,t_fi) + +! ADAPTATION GCM POUR CP(T) +c passage en temperature potentielle + call t2tpot(ip1jmp1*llm,T,teta,pk) +c on assure la periodicite + teta(iip1,:,:) = teta(1,:,:) + +c calcul des champ de vent; passage en vent covariant + write (*,*) 'uold ', uold (1,2,1) ! INFO + call interp_vert + & (uold,var,lmold,llm,apold,bpold,ap,bp, + & psold,(imold+1)*(jmold+1)) + write (*,*) 'var ', var (1,2,1) ! INFO + call interp_horiz(var,us,imold,jmold,iim,jjm,llm, + & rlonuold,rlatvold,rlonu,rlatv) + write (*,*) 'us ', us (1,2,1) ! INFO + + call interp_vert + & (vold,var,lmold,llm, + & apold,bpold,ap,bp,psold,(imold+1)*(jmold+1)) + call interp_horiz(var,vs,imold,jmold,iim,jjm,llm, + & rlonuold,rlatvold,rlonu,rlatv) + call scal_wind(us,vs,unat,vnat) + write (*,*) 'unat ', unat (1,2,1) ! INFO + do l=1,llm + do j = 1, jjp1 + do i=1,iip1 + ucov( i,j,l ) = unat( i,j,l ) * cu(i,j) +c ucov( i,j,l ) = 0 + end do + end do + end do + write (*,*) 'ucov ', ucov (1,2,1) ! INFO +c write(48,*) 'ucov',ucov + do l=1,llm + do j = 1, jjm + do i=1,iim + vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j) +c vcov( i,j,l ) = 0 + end do + vcov( iip1,j,l ) = vcov( 1,j,l ) + end do + end do +c write(49,*) 'ucov',vcov + +c masse: on la recalcule (ps a été ajustée pour conserver la masse totale) + call massdair(p3d,masse) + +c traceurs 3D + do iq = 1, nqtot + call interp_vert(qold(1,1,1,iq),var,lmold,llm, + & apold,bpold,ap,bp,psold,(imold+1)*(jmold+1)) + call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm, + & rlonuold,rlatvold,rlonu,rlatv) + enddo + + print*,"Nouvelles var dynamiques OK" + +c======================================================================= +c Ecriture des restart.nc et restartphy.nc : +c======================================================================= + + call writerestart('restart.nc',tab_cntrl_dyn, + . phis,vcov,ucov,teta,masse,q,ps) + + print*,"restart OK" + + call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm, + . rlat,rlon,tsurf,tsoil,albe,solsw, sollw,dlw, + . radsol, + . zmea,zstd,zsig,zgam,zthe,zpic,zval, + . t_fi) + + print*,"restartphy OK" + + end Index: /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F =================================================================== --- /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F (revision 848) +++ /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F (revision 849) @@ -55,5 +55,5 @@ c====================================================================== USE ioipsl - USE histcom +! USE histcom ! not needed; histcom is included in ioipsl USE infotrac USE control_mod Index: /trunk/LMDZ.VENUS/libf/phyvenus/start2archive.F =================================================================== --- /trunk/LMDZ.VENUS/libf/phyvenus/start2archive.F (revision 849) +++ /trunk/LMDZ.VENUS/libf/phyvenus/start2archive.F (revision 849) @@ -0,0 +1,311 @@ +c======================================================================= + PROGRAM start2archive +c======================================================================= +c +c +c Date: 01/1997 +c ---- +c +c Version Venus: 09/2007 +c Titan: 02/2009 +c +c Objet: Passage des fichiers netcdf d etat initial "start" et +c ----- "startphy" a un fichier netcdf unique "start_archive" +c +c======================================================================= + + USE filtreg_mod + USE infotrac + USE control_mod + + implicit none + +#include "dimensions.h" +#include "paramet.h" +#include "comconst.h" +#include "comdissnew.h" +#include "comvert.h" +#include "comgeom.h" +#include "logic.h" +#include "temps.h" +#include "ener.h" +#include "description.h" +#include "dimsoil.h" +#include "netcdf.inc" + +c----------------------------------------------------------------------- +c Declarations +c----------------------------------------------------------------------- + +c variables dynamiques du GCM +c ----------------------------- + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants + REAL teta(ip1jmp1,llm) ! temperature potentielle + REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes + REAL pks(ip1jmp1) ! exner (f pour filtre) + REAL pk(ip1jmp1,llm) + REAL pkf(ip1jmp1,llm) + REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) + REAL phis(ip1jmp1) ! geopotentiel au sol + REAL masse(ip1jmp1,llm) ! masse de l'atmosphere + REAL ps(ip1jmp1) ! pression au sol + REAL p3d(iip1, jjp1, llm+1) ! pression aux interfaces + +c Variable Physiques (grille physique) +c ------------------------------------ + integer ngridmx,nlayermx + parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm)) + parameter (nlayermx=llm) + + real rlat(ngridmx),rlon(ngridmx) + REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx) + REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx) + real solsw(ngridmx),dlw(ngridmx) + REAL zmea(ngridmx), zstd(ngridmx) + REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx) + REAL zpic(ngridmx), zval(ngridmx) + + INTEGER start,length + PARAMETER (length = 100) + REAL tab_cntrl_fi(length) ! tableau des parametres de startfi + REAL tab_cntrl_dyn(length) ! tableau des parametres de start + INTEGER*4 day_ini_fi + +c Variable naturelle / grille scalaire +c ------------------------------------ + REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm) + REAL tsurfS(ip1jmp1),tsoilS(ip1jmp1,nsoilmx) + real rlatS(ip1jmp1),rlonS(ip1jmp1) + real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1) + real solswS(ip1jmp1),dlwS(ip1jmp1) + real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1) + real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1) + real zvalS(ip1jmp1) + +c Variables intermediaires : vent naturel, mais pas coord scalaire +c---------------------------------------------------------------- + REAL vn(ip1jm,llm),un(ip1jmp1,llm) + +c Autres variables +c ----------------- + REAL ptotal + + CHARACTER*2 str2 + + INTEGER ij, l,i,j,isoil,iq + character*80 fichnom + integer :: ierr + +c Netcdf +c------- + integer varid,dimid + INTEGER nid + +c----------------------------------------------------------------------- +c Initialisations +c----------------------------------------------------------------------- + +c VENUS/TITAN + + iflag_trac = 1 +c----------------------------------------------------------------------- +c Initialisation des traceurs +c --------------------------- +c Choix du nombre de traceurs et du schema pour l advection +c dans fichier traceur.def, par default ou via INCA + call infotrac_init + +c Allocation de la tableau q : champs advectes + allocate(q(ip1jmp1,llm,nqtot)) + +c======================================================================= +c Lecture des donnees +c======================================================================= + + fichnom = 'start.nc' + CALL readstart(fichnom,nqtot,vcov,ucov,teta,q,masse, + . ps,phis,tab_cntrl_dyn) + + fichnom = 'startphy.nc' + CALL readstartphy(fichnom, + . rlat,rlon,tsurf,tsoil, + . albe, solsw, sollw, + . dlw,radsol, + . zmea,zstd,zsig,zgam,zthe,zpic,zval, + . tab_cntrl_fi) + +c----------------------------------------------------------------------- +c Initialisations +c----------------------------------------------------------------------- + + call iniconst + call inigeom + call inifilr + CALL pression(ip1jmp1, ap, bp, ps, p3d) + if (disvert_type==1) then + CALL exner_hyb( ip1jmp1, ps, p3d,alpha,beta,pks, pk, pkf ) + else ! we assume that we are in the disvert_type==2 case + CALL exner_milieu( ip1jmp1, ps, p3d, beta, pks, pk, pkf ) + endif + +c======================================================================= +c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire +c======================================================================= +c Les variables modeles dependent de la resolution. Il faut donc +c eliminer les facteurs responsables de cette dependance +c (pour utiliser newstart) +c======================================================================= + +c----------------------------------------------------------------------- +c Vent (depend de la resolution horizontale) +c----------------------------------------------------------------------- +c +c ucov --> un et vcov --> vn +c un --> us et vn --> vs +c +c----------------------------------------------------------------------- + + call covnat(llm,ucov, vcov, un, vn) + call wind_scal(un,vn,us,vs) + +c----------------------------------------------------------------------- +c Temperature (depend de la resolution verticale => de "sigma.def") +c----------------------------------------------------------------------- +c +c h --> T +c +c----------------------------------------------------------------------- +! ADAPTATION GCM POUR CP(T) + + call tpot2t(ip1jmp1*llm,teta,T,pk) + +c----------------------------------------------------------------------- +c Variable physique +c----------------------------------------------------------------------- +c +c tsurf --> tsurfS +c et autres... +c +c----------------------------------------------------------------------- + + call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS) + call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,rlat,rlatS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,rlon,rlonS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,albe,albeS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,radsol,radsolS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,sollw,sollwS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,solsw,solswS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,dlw,dlwS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zmea,zmeaS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zstd,zstdS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zsig,zsigS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zgam,zgamS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zthe,ztheS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zpic,zpicS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,zval,zvalS) + +c======================================================================= +c Info pour controler +c======================================================================= + + ptotal = 0. + DO j=1,jjp1 + DO i=1,iim + ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g + ENDDO + ENDDO + write(*,*)'Ancienne grille : masse de l''atm :',ptotal + +c----------------------------------------------------------------------- +c Passage de "ptotal" par tab_cntrl_fi +c----------------------------------------------------------------------- + + tab_cntrl_fi(length) = ptotal + +c======================================================================= +c Ecriture dans le fichier "start_archive" +c======================================================================= + +c----------------------------------------------------------------------- +c Ouverture de "start_archive" +c----------------------------------------------------------------------- + + ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid) + +c----------------------------------------------------------------------- +c si "start_archive" n'existe pas: +c 1_ ouverture +c 2_ creation de l'entete dynamique ("ini_archive") +c----------------------------------------------------------------------- +c ini_archive: +c On met dans l'entete le tab_cntrl_dyn (1 a length) +c On y ajoute les valeurs du tab_cntrl_fi (length+1 a 2*length) +c----------------------------------------------------------------------- + + if (ierr.ne.NF_NOERR) then + write(*,*)'OK, Could not open file "start_archive.nc"' + write(*,*)'So let s create a new "start_archive"' + ierr = NF_CREATE('start_archive.nc', NF_CLOBBER, nid) + call ini_archive(nid,phis,tab_cntrl_dyn,tab_cntrl_fi) + else + write(*,*)'Attention, start_archive.nc existe déjà...' + call abort + endif + +c----------------------------------------------------------------------- +c Ecriture des champs +c----------------------------------------------------------------------- + + call write_archive(nid,'u','Vent zonal','m.s-1',3,us) + call write_archive(nid,'v','Vent merid','m.s-1',3,vs) + call write_archive(nid,'temp','temperature','K',3,T) +c----------------------------------------------------------------------- +c Ecriture du champs q ( q[1,nqtot] ) +c----------------------------------------------------------------------- + do iq=1,nqtot + write(str2,'(i2.2)') iq + call write_archive(nid,tname(iq),'tracer','kg/kg', + . 3,q(1,1,iq)) + end do +c----------------------------------------------------------------------- + call write_archive(nid,'masse','Masse','kg',3,masse) + call write_archive(nid,'ps','Psurf','Pa',2,ps) + call write_archive(nid,'tsurf','surf T','K',2,tsurfS) +c----------------------------------------------------------------------- +c Ecriture du champs tsoil ( Tsoil[1,nsoilmx] ) +c----------------------------------------------------------------------- +c "tsoil" Temperature au sol definie dans nsoilmx couches dans le sol +c Les nsoilmx couches sont lues comme nsoilmx champs +c nommees Tsoil[1,nsoilmx] + do isoil=1,nsoilmx + write(str2,'(i2.2)') isoil + call write_archive(nid,'Tsoil'//str2,'Ground Temperature ', + . 'K',2,tsoilS(1,isoil)) + enddo +c----------------------------------------------------------------------- + call write_archive(nid,'rlat','Latitude','rad',2,rlatS) + call write_archive(nid,'rlon','Longitude','rad',2,rlonS) + call write_archive(nid,'albe','Albedo','',2,albeS) + call write_archive(nid,'radsol', + . 'Net flux at surface','W m-2',2,radsolS) + call write_archive(nid,'sollw', + . 'LW flux at surface','W m-2',2,sollwS) + call write_archive(nid,'solsw', + . 'SW flux at surface','W m-2',2,solswS) + call write_archive(nid,'dlw','LW derive','?',2,dlwS) + call write_archive(nid,'zmea','param oro sous-maille','m',2,zmeaS) + call write_archive(nid,'zstd','param oro sous-maille','m',2,zstdS) + call write_archive(nid,'zsig','param oro sous-maille','m',2,zsigS) + call write_archive(nid,'zgam','param oro sous-maille','m',2,zgamS) + call write_archive(nid,'zthe','param oro sous-maille','m',2,ztheS) + call write_archive(nid,'zpic','param oro sous-maille','m',2,zpicS) + call write_archive(nid,'zval','param oro sous-maille','m',2,zvalS) + + ierr=NF_CLOSE(nid) + +c----------------------------------------------------------------------- +c Fin +c----------------------------------------------------------------------- + + end Index: /trunk/LMDZ.VENUS/libf/phyvenus/testphys1d.F =================================================================== --- /trunk/LMDZ.VENUS/libf/phyvenus/testphys1d.F (revision 848) +++ /trunk/LMDZ.VENUS/libf/phyvenus/testphys1d.F (revision 849) @@ -1,4 +1,6 @@ PROGRAM testphys1d + + use control_mod IMPLICIT NONE @@ -27,5 +29,5 @@ #include "dimsoil.h" #include "comcstfi.h" -#include "control.h" +!#include "control.h" #include "comvert.h" #include "netcdf.inc" @@ -226,5 +228,5 @@ c Vertical Coordinates (hybrids) c """""""""""""""""""" - CALL disvert + CALL disvert_noterre c Calcul au milieu des couches : Vient de la version Mars @@ -474,5 +476,5 @@ c*********************************************************************** -#include "../dyn3d/disvert.F" +#include "../dyn3d/disvert_noterre.F" #include "../dyn3d/abort_gcm.F" #include "../dyn3d/dump2d.F"