SUBROUTINE lect_start_archive(ngrid,nlayer, & date,tsurf,tsoil,emis,q2, & t,ucov,vcov,ps,h,phisold_newgrid, & q,qsurf,surfith,nid, & rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice, & du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress) USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat USE tracer_h, ONLY: igcm_co2_ice USE infotrac, ONLY: tname, nqtot ! USE slab_ice_h, ONLY: noceanmx USE ocean_slab_mod, ONLY: nslay USE callkeys_mod, ONLY: ok_slab_ocean USE comvert_mod, ONLY: ap,bp,aps,bps,preff USE comconst_mod, ONLY: kappa,g,pi c======================================================================= c c Routine to load variables from the "start_archive.nc" file c c======================================================================= implicit none include "dimensions.h" include "paramet.h" include "comgeom2.h" include "netcdf.inc" c======================================================================= c Declarations c======================================================================= INTEGER,INTENT(IN) :: ngrid, nlayer c Old variables dimensions (from file) c------------------------------------ INTEGER imold,jmold,lmold,nsoilold,nqold c Variables pour les lectures des fichiers "ini" c-------------------------------------------------- ! INTEGER sizei, integer timelen,dimid ! INTEGER length ! parameter (length = 100) INTEGER tab0 INTEGER isoil,iq,iqmax CHARACTER*2 str2 REAL,INTENT(OUT) :: date INTEGER memo ! character (len=50) :: tmpname c Variable histoire c------------------ REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1) REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) c Physique sur grille scalaire c---------------------------- c variable physique c------------------ REAL,INTENT(OUT) :: tsurf(ngrid) ! surface temperature REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx) ! soil temperature REAL co2ice(ngrid) ! CO2 ice layer REAL,INTENT(OUT) :: emis(ngrid) REAL,INTENT(OUT) :: q2(ngrid,llm+1),qsurf(ngrid,nqtot) REAL,INTENT(OUT) :: tslab(ngrid,nslay) REAL ,INTENT(OUT) ::rnat(ngrid),pctsrf_sic(ngrid) REAL,INTENT(OUT) :: tsea_ice(ngrid),tice(ngrid),sea_ice(ngrid) c REAL phisfi(ngrid) REAL,INTENT(OUT):: du_nonoro_gwd(ngrid,llm) REAL,INTENT(OUT):: dv_nonoro_gwd(ngrid,llm) REAL,INTENT(OUT):: east_gwstress(ngrid,llm) REAL,INTENT(OUT):: west_gwstress(ngrid,llm) INTEGER i,j,l INTEGER,INTENT(IN) :: nid INTEGER :: nvarid c REAL year_day,periheli,aphelie,peri_day c REAL obliquit,z0,emin_turb,lmixmin c REAL emissiv,emisice(2),albedice(2) c REAL iceradius(2) , dtemisice(2) integer ierr integer, dimension(4) :: start,count c Variable nouvelle grille naturelle au point scalaire c------------------------------------------------------ real us(iip1,jjp1,llm),vs(iip1,jjp1,llm) REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1) REAL,INTENT(OUT) :: t(iip1,jjp1,llm) real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx) real inertiedatS(iip1,jjp1,nsoilmx) real co2iceS(iip1,jjp1) real emisS(iip1,jjp1) REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot) real tslabS(iip1,jjp1,nslay) real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1),ticeS(iip1,jjp1) real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1) real du_nonoro_gwdS(iip1,jjp1,llm),dv_nonoro_gwdS(iip1,jjp1,llm) real east_gwstressS(iip1,jjp1,llm),west_gwstressS(iip1,jjp1,llm) real ptotal, co2icetotal c Var intermediaires : vent naturel, mais pas coord scalaire c----------------------------------------------------------- real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) c Variable de l'ancienne grille c--------------------------------------------------------- real, dimension(:), allocatable :: timelist real, dimension(:), allocatable :: rlonuold, rlatvold real, dimension(:), allocatable :: rlonvold, rlatuold real, dimension(:), allocatable :: apsold,bpsold real, dimension(:), allocatable :: mlayerold real, dimension(:,:,:), allocatable :: uold,vold,told,q2old real, dimension(:,:,:), allocatable :: tsoilold,qsurfold real, dimension(:,:,:),allocatable :: tsoiloldnew ! tsoiloldnew: old soil values, but along new subterranean grid real, dimension(:,:,:), allocatable :: inertiedatold ! inertiedatoldnew: old inertia values, but along new subterranean grid real, dimension(:,:,:), allocatable :: inertiedatoldnew real, dimension(:,:), allocatable :: psold,phisold real, dimension(:,:), allocatable :: co2iceold real, dimension(:,:), allocatable :: tsurfold real, dimension(:,:), allocatable :: emisold real, dimension(:,:,:,:), allocatable :: qold real, dimension(:,:,:), allocatable :: tslabold real, dimension(:,:), allocatable :: rnatold,pctsrf_sicold real, dimension(:,:), allocatable :: tsea_iceold,ticeold real, dimension(:,:), allocatable :: sea_iceold real,allocatable :: du_nonoro_gwdold(:,:,:) real,allocatable :: dv_nonoro_gwdold(:,:,:) real,allocatable :: east_gwstressold(:,:,:) real,allocatable :: west_gwstressold(:,:,:) real tab_cntrl(100) real ptotalold, co2icetotalold logical :: olddepthdef=.false. ! flag ! olddepthdef=.true. if soil depths are in 'old' (unspecified) format logical :: depthinterpol=.false. ! flag ! depthinterpol=.true. if interpolation will be requiered logical :: therminertia_3D=.true. ! flag ! therminertia_3D=.true. if thermal inertia is 3D and read from datafile c Variable intermediaires iutilise pour l'extrapolation verticale c---------------------------------------------------------------- real, dimension(:,:,:), allocatable :: var,varp1 real, dimension(:), allocatable :: oldgrid, oldval real, dimension(:), allocatable :: newval real,intent(out) :: surfith(iip1,jjp1) ! surface thermal inertia ! surface thermal inertia at old horizontal grid resolution real, dimension(:,:), allocatable :: surfithold character(len=30) :: txt ! to store some text c======================================================================= ! 0. Preliminary stuff !----------------------------------------------------------------------- ! 1. Read data dimensions (i.e. size and length) !----------------------------------------------------------------------- ! 1.2 Read the various dimension lengths of data in file ierr= NF_INQ_DIMID(nid,"Time",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"temps",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,timelen) ierr= NF_INQ_DIMID(nid,"latitude",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"rlatu",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,jmold) jmold=jmold-1 ierr= NF_INQ_DIMID(nid,"longitude",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"rlonv",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,imold) imold=imold-1 ierr= NF_INQ_DIMID(nid,"altitude",dimid) if (ierr.ne.NF_NOERR) then ierr= NF_INQ_DIMID(nid,"sig_s",dimid) endif ierr= NF_INQ_DIMLEN(nid,dimid,lmold) nqold=0 do write(str2,'(i2.2)') nqold+1 ierr= NF_INQ_VARID(nid,'q'//str2,dimid) ! write(*,*) 'q'//str2 if (ierr.eq.NF_NOERR) then nqold=nqold+1 else exit endif enddo ! 1.2.1 find out the # of subsurface_layers nsoilold=0 !dummy initialisation ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid) if (ierr.eq.NF_NOERR) then ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold) if (ierr.ne.NF_NOERR) then write(*,*)'lec_start_archive: ', & 'Failed reading subsurface_layers length' endif else write(*,*)"lec_start_archive: did not find subsurface_layers" endif if (nsoilold.eq.0) then ! 'old' archive format; ! must use Tg//str2 fields to compute nsoilold write(*,*)"lec_start_archive: building nsoilold from Tg fields" do write(str2,'(i2.2)') nsoilold+1 ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid) if (ierr.eq.NF_NOERR) then nsoilold=nsoilold+1 else exit endif enddo endif if (nsoilold.ne.nsoilmx) then ! interpolation will be required depthinterpol=.true. endif ! 1.3 Report dimensions write(*,*) "Start_archive dimensions:" write(*,*) "longitude: ",imold write(*,*) "latitude: ",jmold write(*,*) "altitude: ",lmold write(*,*) "tracers: ",nqold write(*,*) "subsurface_layers: ",nsoilold if (depthinterpol) then write(*,*) " => Warning, nsoilmx= ",nsoilmx write(*,*) ' which implies that you want subterranean interpola &tion.' write(*,*) ' Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold endif write(*,*) "time lenght: ",timelen write(*,*) !----------------------------------------------------------------------- ! 2. Allocate arrays to store datasets !----------------------------------------------------------------------- allocate(timelist(timelen)) allocate(rlonuold(imold+1), rlatvold(jmold)) allocate(rlonvold(imold+1), rlatuold(jmold+1)) allocate (apsold(lmold),bpsold(lmold)) allocate(uold(imold+1,jmold+1,lmold)) allocate(vold(imold+1,jmold+1,lmold)) allocate(told(imold+1,jmold+1,lmold)) allocate(psold(imold+1,jmold+1)) allocate(phisold(imold+1,jmold+1)) allocate(qold(imold+1,jmold+1,lmold,nqtot)) allocate(co2iceold(imold+1,jmold+1)) allocate(tsurfold(imold+1,jmold+1)) allocate(emisold(imold+1,jmold+1)) allocate(q2old(imold+1,jmold+1,lmold+1)) ! allocate(tsoilold(imold+1,jmold+1,nsoilmx)) allocate(tsoilold(imold+1,jmold+1,nsoilold)) allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx)) allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx)) ! surface thermal inertia at old horizontal grid resolution allocate(surfithold(imold+1,jmold+1)) allocate(mlayerold(nsoilold)) allocate(qsurfold(imold+1,jmold+1,nqtot)) allocate(tslabold(imold+1,jmold+1,nslay)) allocate(rnatold(imold+1,jmold+1)) allocate(pctsrf_sicold(imold+1,jmold+1)) allocate(tsea_iceold(imold+1,jmold+1)) allocate(ticeold(imold+1,jmold+1)) allocate(sea_iceold(imold+1,jmold+1)) allocate(du_nonoro_gwdold(imold+1,jmold+1,lmold)) allocate(dv_nonoro_gwdold(imold+1,jmold+1,lmold)) allocate(east_gwstressold(imold+1,jmold+1,lmold)) allocate(west_gwstressold(imold+1,jmold+1,lmold)) allocate(var (imold+1,jmold+1,llm)) allocate(varp1 (imold+1,jmold+1,llm+1)) write(*,*) 'lect_start_archive: q2',ngrid,llm+1 write(*,*) 'lect_start_archive: q2S',iip1,jjp1,llm+1 write(*,*) 'lect_start_archive: q2old',imold+1,jmold+1,lmold+1 !----------------------------------------------------------------------- ! 3. Read time-independent data !----------------------------------------------------------------------- C----------------------------------------------------------------------- c 3.1. Lecture du tableau des parametres du run c (pour la lecture ulterieure de "ptotalold" et "co2icetotalold") c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "controle", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Lect_start_archive: champ not found" CALL abort ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) #else ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Lecture echoue pour " CALL abort ENDIF c tab0 = 50 c----------------------------------------------------------------------- c 3.2 Lecture des longitudes et latitudes c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "rlonv", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" 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*, "lect_start_archive: Failed loading " CALL abort ENDIF c ierr = NF_INQ_VARID (nid, "rlatu", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" 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*, "lect_start_archive: Failed loading " CALL abort ENDIF c ierr = NF_INQ_VARID (nid, "rlonu", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" 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*, "lect_start_archive: Failed loading " CALL abort ENDIF c ierr = NF_INQ_VARID (nid, "rlatv", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" 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*, "lect_start_archive: Failed loading " CALL abort ENDIF c c----------------------------------------------------------------------- c 3.3. Lecture des niveaux verticaux c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "aps", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" apsold=0 PRINT*, " set to 0" ELSE #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, apsold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " ENDIF ENDIF c ierr = NF_INQ_VARID (nid, "bps", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" PRINT*, "It must be an old start_archive, lets look for sig_s" ierr = NF_INQ_VARID (nid, "sig_s", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "Nothing to do..." CALL abort ENDIF ENDIF #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold) #else ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold) #endif IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Failed loading " CALL abort END IF c----------------------------------------------------------------------- c 3.4 Read Soil layers depths c----------------------------------------------------------------------- ierr=NF_INQ_VARID(nid,"soildepth",nvarid) if (ierr.ne.NF_NOERR) then write(*,*)'lect_start_archive: Could not find ' write(*,*)' => Assuming this is an archive in old format' olddepthdef=.true. depthinterpol=.true. ! this is how soil depth was defined in ye old days do isoil=1,nsoilold mlayerold(isoil)=sqrt(887.75/3.14)*((2.**(isoil-0.5))-1.) enddo else #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold) #else ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold) #endif if (ierr .NE. NF_NOERR) then PRINT*, "lect_start_archive: Failed reading " CALL abort endif endif !of if(ierr.ne.NF_NOERR) ! Read (or build) mlayer() if (depthinterpol) then ! Build (default) new soil depths (mlayer(:) is in comsoil.h), ! as in soil_settings.F write(*,*)' => Building default soil depths' do isoil=0,nsoilmx-1 mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) enddo write(*,*)' => mlayer: ',mlayer ! Also build (default) new soil interlayer depth layer(:) do isoil=1,nsoilmx layer(isoil)=sqrt(mlayer(0)*mlayer(1))* & ((mlayer(1)/mlayer(0))**(isoil-1)) enddo write(*,*)' => layer: ',layer else ! read mlayer() from file #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer) #else ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer) #endif if (ierr .NE. NF_NOERR) then PRINT*, "lect_start_archive: Failed reading " CALL abort endif endif ! of if (depthinterpol) c----------------------------------------------------------------------- c 3.5 Read Soil thermal inertia c----------------------------------------------------------------------- ierr=NF_INQ_VARID(nid,"inertiedat",nvarid) if (ierr.ne.NF_NOERR) then write(*,*)'lect_start_archive: Could not find ' write(*,*)' => Assuming this is an archive in old format' therminertia_3D=.false. write(*,*)' => Thermal inertia will be read from reference file' volcapa=1.e6 write(*,*)' and soil volumetric heat capacity is set to ', & volcapa else #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold) #else ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold) #endif if (ierr .NE. NF_NOERR) then PRINT*, "lect_start_archive: Failed reading " CALL abort endif endif c----------------------------------------------------------------------- c 3.6 Lecture geopotentiel au sol c----------------------------------------------------------------------- c ierr = NF_INQ_VARID (nid, "phisinit", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field not found" 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*, "lect_start_archive: Failed loading " CALL abort ENDIF C----------------------------------------------------------------------- c lecture de "ptotalold" et "co2icetotalold" c----------------------------------------------------------------------- ptotalold = tab_cntrl(tab0+49) co2icetotalold = tab_cntrl(tab0+50) c----------------------------------------------------------------------- c 4. Lecture du temps et choix c----------------------------------------------------------------------- c lecture du temps c ierr = NF_INQ_DIMID (nid, "Time", nvarid) IF (ierr .NE. NF_NOERR) THEN ierr = NF_INQ_DIMID (nid, "temps", nvarid) IF (ierr .NE. NF_NOERR) THEN PRINT*, "lect_start_archive: Field