C====================================================================== PROGRAM newstart c======================================================================= c c c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick c ------ c Derniere modif : 12/03 c c c Objet: Create or modify the initial state for the LMD Mars GCM c ----- (fichiers NetCDF start et startfi) c c c======================================================================= implicit none #include "dimensions.h" #include "dimphys.h" #include "surfdat.h" #include "comsoil.h" #include "planete.h" #include "paramet.h" #include "comconst.h" #include "comvert.h" #include "comgeom2.h" #include "control.h" #include "logic.h" #include "description.h" #include "ener.h" #include "temps.h" #include "lmdstd.h" #include "comdissnew.h" #include "clesph0.h" #include "serre.h" #include "netcdf.inc" #include "advtrac.h" #include "tracer.h" c======================================================================= c Declarations c======================================================================= c Variables dimension du fichier "start_archive" c------------------------------------ CHARACTER relief*3 c Variables pour les lectures NetCDF des fichiers "start_archive" c-------------------------------------------------- INTEGER nid_dyn, nid_fi,nid,nvarid INTEGER length parameter (length = 100) INTEGER tab0 INTEGER NB_ETATMAX parameter (NB_ETATMAX = 100) REAL date REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec c Variable histoire c------------------ REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants REAL phis(iip1,jjp1) REAL q(iip1,jjp1,llm,nqmx) ! champs advectes c autre variables dynamique nouvelle grille c------------------------------------------ REAL pks(iip1,jjp1) REAL w(iip1,jjp1,llm+1) REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) ! REAL dv(ip1jm,llm),du(ip1jmp1,llm) ! REAL dh(ip1jmp1,llm),dp(ip1jmp1) REAL phi(iip1,jjp1,llm) integer klatdat,klongdat PARAMETER (klatdat=180,klongdat=360) c Physique sur grille scalaire c---------------------------- real zmeaS(iip1,jjp1),zstdS(iip1,jjp1) real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1) c variable physique c------------------ REAL tsurf(ngridmx) ! surface temperature REAL tsoil(ngridmx,nsoilmx) ! soil temperature REAL lati(ngridmx) ! REAL co2ice(ngridmx) ! CO2 ice layer REAL emis(ngridmx) ! surface emissivity real emisread ! added by RW REAL qsurf(ngridmx,nqmx) REAL q2(ngridmx,nlayermx+1) ! REAL rnaturfi(ngridmx) real alb(iip1,jjp1),albfi(ngridmx) ! albedos real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D) real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D) REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) INTEGER i,j,l,isoil,ig,idum ! real mugaz ! molar mass of the atmosphere. !TB15 : in comconst.h integer ierr c Variables on the new grid along scalar points c------------------------------------------------------ ! REAL p(iip1,jjp1) REAL t(iip1,jjp1,llm) REAL tset(iip1,jjp1,llm) real phisold_newgrid(iip1,jjp1) ! REAL tanh(ip1jmp1) REAL :: teta(iip1, jjp1, llm) REAL :: pk(iip1,jjp1,llm) REAL :: pkf(iip1,jjp1,llm) REAL :: ps(iip1, jjp1) REAL :: ppa(iip1,jjp1,llm) REAL :: masse(iip1,jjp1,llm) REAL :: xpn,xps,xppn(iim),xpps(iim) REAL :: p3d(iip1, jjp1, llm+1) REAL :: beta(iip1,jjp1,llm) ! REAL dteta(ip1jmp1,llm) c Variable de l'ancienne grille c------------------------------ real time real tab_cntrl(100) real tab_cntrl_bis(100) c variables diverses c------------------- real choix_1,choice,pp character*80 fichnom integer Lmodif,iq,thermo character modif*20 real z_reel(iip1,jjp1) ! real z_reel(ngridmx) real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove,tsurf_bb,geop real ptoto,pcap,patm,airetot,ptotn,patmn,psea character(len=50) :: surfacefile ! "surface.nc" (or equivalent file) ! real ssum character*1 yes logical :: flagtset=.false. , flagps0=.false. real val, val1, val2, val3, val4 ! to store temporary variables real :: iceith=2000 ! thermal inertia of subterranean ice integer iref,jref,iref1,jref1,iref2,jref2 INTEGER :: itau INTEGER :: nq,numvanle character(len=20) :: txt ! to store some text integer :: count ! MONS data: ! real :: MONS_Hdn(iip1,jjp1) ! Hdn: %WEH=Mass fraction of H2O ! real :: MONS_d21(iip1,jjp1) ! ice table "depth" (in kg/m2) ! ! coefficient to apply to convert d21 to 'true' depth (m) ! real :: MONS_coeff ! real :: MONS_coeffS ! coeff for southern hemisphere ! real :: MONS_coeffN ! coeff for northern hemisphere !! real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth c Special Pluto Map from Lellouch & Grundy c------------------------------------------ integer,parameter :: im_plu=59 integer,parameter :: jm_plu=30 real latv_plu(jm_plu),lonu_plu(im_plu+1) real map_pluto_dat(im_plu,jm_plu+1) real alb_pluto_dat(im_plu,jm_plu+1) real N2_pluto_dat(im_plu,jm_plu+1) real CH4_pluto_dat(im_plu,jm_plu+1) real CO_pluto_dat(im_plu,jm_plu+1) real ith_pluto_dat(im_plu,jm_plu+1) real alb_pluto_rein(im_plu+1,jm_plu+1) real N2_pluto_rein(im_plu+1,jm_plu+1) real CH4_pluto_rein(im_plu+1,jm_plu+1) real CO_pluto_rein(im_plu+1,jm_plu+1) real ith_pluto_rein(im_plu+1,jm_plu+1) real N2_pluto_gcm(iip1,jjp1) real CH4_pluto_gcm(iip1,jjp1) real CO_pluto_gcm(iip1,jjp1) real alb_pluto_gcm(iip1,jjp1),alb_pluto_fi(ngridmx) real ith_pluto_gcm(iip1,jjp1),ithf(ngridmx) c sortie visu pour les champs dynamiques c--------------------------------------- ! INTEGER :: visuid ! real :: time_step,t_ops,t_wrt ! CHARACTER*80 :: visu_file ! preff = 2. ! Pluto ! 610. ! for Mars, instead of 101325. (Earth) ! pa = 0.5 ! Pluto ! 20 ! for Mars, instead of 500 (Earth) c======================================================================= c Choice of the start file(s) to use c======================================================================= write(*,*) 'From which kind of files do you want to create new', . 'start and startfi files' write(*,*) ' 0 - from a file start_archive' write(*,*) ' 1 - from files start and startfi' c----------------------------------------------------------------------- c Open file(s) to modify (start or start_archive) c----------------------------------------------------------------------- DO read(*,*,iostat=ierr) choix_1 if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT ENDDO c Open start_archive c ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (choix_1.eq.0) then write(*,*) 'Creating start files from:' write(*,*) './start_archive.nc' write(*,*) fichnom = 'start_archive.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF tab0 = 50 Lmodif = 1 c OR open start and startfi files c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else write(*,*) 'Creating start files from:' write(*,*) './start.nc and ./startfi.nc' write(*,*) fichnom = 'start.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF fichnom = 'startfi.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening file:',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF tab0 = 0 Lmodif = 0 endif c----------------------------------------------------------------------- c Lecture du tableau des parametres du run (pour la dynamique) c----------------------------------------------------------------------- if (choix_1.eq.0) then write(*,*) 'reading 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 else if (choix_1.eq.1) then write(*,*) 'reading tab_cntrl START' c ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl) #else ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl) #endif c write(*,*) 'reading tab_cntrl STARTFI' c ierr = NF_INQ_VARID (nid_fi, "controle", nvarid) #ifdef NC_DOUBLE ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis) #else ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis) #endif c do i=1,50 tab_cntrl(i+50)=tab_cntrl_bis(i) enddo endif c TB15 : c do i=1,size(tab_cntrl) c write(*,*) i,tab_cntrl(i) c enddo c----------------------------------------------------------------------- c Initialisation des constantes dynamique c----------------------------------------------------------------------- kappa = tab_cntrl(9) etot0 = tab_cntrl(12) ptot0 = tab_cntrl(13) ztot0 = tab_cntrl(14) stot0 = tab_cntrl(15) ang0 = tab_cntrl(16) write(*,*) "Newstart: in nc: kappa,etot0,ptot0,ztot0,stot0,ang0" write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0 ! for vertical coordinate preff=tab_cntrl(18) ! reference surface pressure pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0 if (preff.eq.0) then preff=2. ! Pluto pa=0.5 ! Pluto endif write(*,*) "Newstart: preff=",preff," pa=",pa yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*) "Change the values of preff and pa? (y/n)" read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*)"New value of reference surface pressure preff?" write(*,*)" (for Mars, typically preff=610)" read(*,*) preff write(*,*)"New value of reference pressure pa for purely" write(*,*)"pressure levels (for hybrid coordinates)?" write(*,*)" (for Mars, typically pa=20)" read(*,*) pa endif c----------------------------------------------------------------------- c Lecture du tab_cntrl et initialisation des constantes physiques c - pour start: Lmodif = 0 => pas de modifications possibles c (modif dans le tabfi de readfi + loin) c - pour start_archive: Lmodif = 1 => modifications possibles c----------------------------------------------------------------------- if (choix_1.eq.0) then call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad, . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) else if (choix_1.eq.1) then call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad, . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) endif rad = p_rad omeg = p_omeg g = p_g cpp = p_cpp mugaz = p_mugaz daysec = p_daysec kappa= 8.314511E+0 *1000.E+0/(mugaz*cpp) c======================================================================= c INITIALISATIONS DIVERSES c======================================================================= ! Load tracer names: call iniadvtrac(nq,numvanle) ! tnom(:) now contains tracer names ! Initialize global tracer indexes (stored in tracer.h) call initracer() ! Load parameters from run.def file CALL defrun_new( 99, .TRUE. ) CALL iniconst CALL inigeom idum=-1 idum=0 c Initialisation coordonnees /aires c ------------------------------- ! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h" ! rlatu() and rlonv() are given in radians latfi(1)=rlatu(1) lonfi(1)=0. DO j=2,jjm DO i=1,iim latfi((j-2)*iim+1+i)=rlatu(j) lonfi((j-2)*iim+1+i)=rlonv(i) ENDDO ENDDO latfi(ngridmx)=rlatu(jjp1) lonfi(ngridmx)=0. CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) c======================================================================= c lecture topographie, albedo, inertie thermique, relief sous-maille c======================================================================= if (choix_1.ne.1) then ! pour ne pas avoir besoin du fichier ! surface.dat dans le cas des start c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla')) c write(*,*) c write(*,*) 'choix du relief (mola,pla)' c write(*,*) '(Topographie MGS MOLA, plat)' c read(*,fmt='(a3)') relief c relief="mola" c enddo write(*,*) '' write(*,*) 'Please choose:' write(*,*) '1- use of surface file netcdf' write(*,*) '2- define topography manually' write(*,*) '3- skip this part' read(*,*) choice if (choice.eq.1) then write(*,*) ' enter surface file: (e.g. "surface_mars.nc)' read(*,fmt='(a50)') surfacefile CALL datareadnc(relief,surfacefile,phis,alb,surfith, & zmeaS,zstdS,zsigS,zgamS,ztheS) else if (choice.eq.2) then ! specific case when not using a "surface.nc" file 302 write(*,*) 'Do you want a flat terrain? (y/n)' read(*,fmt='(a)') yes if (yes.eq.'y') then phis(:,:)=0 else !! Creation topography : !! Choice longitudes val1 and val2 write(*,*)'Choice range longitude (in deg.), and topography' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val1,val2 write(*,*)'Values are:',val1,val2 if (ierr.eq.0) then ! got a value ! do a sanity check if((val1.lt.-180.).or.(val1.gt.180.)) then write(*,*)'Longitude should be between -180 and 180 deg!' ierr=1 else if((val2.lt.-180.).or.(val2.gt.180.)) then write(*,*)'Longitude should be between -180 and 180 deg!' ierr=1 else if(val1.ge.val2) then write(*,*)'Longitude 1 should be lower than longitude 2 !' ierr=1 else ! find corresponding iref (nearest longitude) ! note: rlonv(:) contains increasing values of longitude ! starting from -PI to PI do i=1,iip1 ! get val1 write(*,*) 'val1:',rlonv(i)*180./pi,rlonv(i+1)*180./pi if ((rlonv(i)*180./pi.le.val1).and. & (rlonv(i+1)*180./pi.ge.val1)) then ! find which grid point is nearest to val: if (abs(rlonv(i)*180./pi-val1).le. & abs((rlonv(i+1)*180./pi-val1))) then iref1=i else iref1=i+1 endif write(*,*)'Will use nearest grid longitude 1:', & rlonv(iref1)*180./pi endif enddo ! of do j=1,jjp1 do i=1,iip1 ! get val2 if ((rlonv(i)*180./pi.le.val2).and. & (rlonv(i+1)*180./pi.ge.val2)) then ! find which grid point is nearest to val: if (abs(rlonv(i)*180./pi-val2).le. & abs((rlonv(i+1)*180./pi-val2))) then iref2=i else iref2=i+1 endif write(*,*)'Will use nearest grid longitude 2:', & rlonv(iref2)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif ! of if((val.lt.0.).or.(val.gt.90)) enddo ! of do while !! Choice latitudes val1 and val2 write(*,*)'Choice range latitudes (in deg.), and topography' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val1,val2 write(*,*)'Values are:',val1,val2 if (ierr.eq.0) then ! got a value ! do a sanity check if((val1.lt.-90.).or.(val1.gt.90.)) then write(*,*)'Latitude should be between -90 and 90 deg!' ierr=1 else if((val2.lt.-90.).or.(val2.gt.90.)) then write(*,*)'Latitude should be between -90 and 90 deg!' ierr=1 else if(val1.ge.val2) then write(*,*)'Latitude 1 should be higuer than latitude 2 !' ierr=1 else ! find corresponding jref (nearest latitude) ! note: rlatu(:) contains decreasing values of latitude ! starting from PI/2 to -PI/2 do j=1,jjp1 ! get val1 write(*,*) 'val1:',rlatu(j)*180./pi,rlatu(j+1)*180./pi if ((rlatu(j)*180./pi.ge.val1).and. & (rlatu(j+1)*180./pi.le.val1)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val1).le. & abs((rlatu(j+1)*180./pi-val1))) then jref1=j else jref1=j+1 endif write(*,*)'Will use nearest grid latitude:', & rlatu(jref1)*180./pi endif enddo ! of do j=1,jjp1 do j=1,jjp1 ! get val2 write(*,*) 'val2:',rlatu(j)*180./pi,rlatu(j+1)*180./pi if ((rlatu(j)*180./pi.ge.val2).and. & (rlatu(j+1)*180./pi.le.val2)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val2).le. & abs((rlatu(j+1)*180./pi-val2))) then jref2=j else jref2=j+1 endif write(*,*)'Will use nearest grid latitude:', & rlatu(jref2)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif ! enddo ! of do while write(*,*) 'Choice of topography (m) (eg: 1000 or -1000) ?' 301 read(*,*,iostat=ierr) val if(ierr.ne.0) goto 301 write(*,*) 'gravity g is : ',g geop=g*val c phis(:,:)=0 phis(iref1:iref2,jref2:jref1)=geop write(*,*) 'phis=',phis write(*,*) 'iip1=',iip1,'jjp1=',jjp1 write(*,*) 'iref1=',iref1,'iref2=',iref2 write(*,*) 'jref1=',jref1,'jref2=',jref2 write(*,*) 'Do you want another topography choice ?' read(*,fmt='(a)') yes if (yes.eq.'y') goto 302 endif ! 'yes flat terrain' zmeaS(:,:)=0 zstdS(:,:)=0 zsigS(:,:)=0 zgamS(:,:)=0 ztheS(:,:)=0 yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'Do you want to change soil albedo and IT (y/n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) "Enter value of albedo of the bare ground:" read(*,*) alb(1,1) alb(:,:)=alb(1,1) write(*,*) "Enter value of thermal inertia of soil:" read(*,*) surfith(1,1) surfith(:,:)=surfith(1,1) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) endif !yes for alb and IT changes else yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'Do you want to change soil albedo and IT (y/n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) "Enter value of albedo of the bare ground:" read(*,*) alb(1,1) alb(:,:)=alb(1,1) write(*,*) "Enter value of thermal inertia of soil:" read(*,*) surfith(1,1) surfith(:,:)=surfith(1,1) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) endif !yes for alb and IT changes endif ! of if (choice) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe) endif ! of if (choix_1.ne.1) c======================================================================= c Lecture des fichiers (start ou start_archive) c======================================================================= if (choix_1.eq.0) then write(*,*) 'Reading file START_ARCHIVE' CALL lect_start_archive(date,tsurf,tsoil,emis,q2, . t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf, & surfith,nid) write(*,*) "OK, read start_archive file" ! copy soil thermal inertia ithfi(:,:)=inertiedat(:,:) ierr= NF_CLOSE(nid) else if (choix_1.eq.1) then ! c'est l'appel a tabfi de phyeta0 qui ! permet de changer les valeurs du ! tab_cntrl Lmodif=1 tab0=0 Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0 write(*,*) 'Reading file START' fichnom = 'start.nc' CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse, . ps,phis,time) write(*,*) 'Reading file STARTFI' fichnom = 'startfi.nc' CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx, . day_ini,time, . tsurf,tsoil,emis,q2,qsurf) ! copy albedo and soil thermal inertia do i=1,ngridmx albfi(i) = albedodat(i) do j=1,nsoilmx ithfi(i,j) = inertiedat(i,j) enddo ! build a surfithfi(:) using 1st layer of ithfi(:), which might ! be neede later on if reinitializing soil thermal inertia surfithfi(i)=ithfi(i,1) enddo else CALL exit(1) endif dtvr = daysec/FLOAT(day_step) dtphys = dtvr * FLOAT(iphysiq) c======================================================================= c c======================================================================= do ! infinite loop on list of changes write(*,*) write(*,*) write(*,*) 'List of possible changes :' write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) write(*,*) 'flat : no topography ("aquaplanet")' write(*,*) 'bilball : uniform albedo and thermal inertia' write(*,*) 'coldspole : cold subsurface and high albedo at S.pole' write(*,*) 'qname : change tracer name' write(*,*) 'q=0 : ALL tracer =zero' write(*,*) 'q=x : give a specific uniform value to one tracer' write(*,*) 'ini_q : tracers initialisation for chemistry, water an $d ice ' write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and $ice ' write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on $ly ' write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45' write(*,*) 'nitrogencapn : N2 ice on permanent N polar cap ' write(*,*) 'nitrogencaps : N2 ice on permanent S polar cap ' write(*,*) 'oborealis : H2O ice across Vastitas Borealis' write(*,*) 'wetstart : start with a wet atmosphere' write(*,*) 'isotherm : Isothermal Temperatures, wind set to zero' write(*,*) 'initsoil : initialize soil Temperatures' write(*,*) 'coldstart : Start X K above the CO2 frost point and $set wind to zero (assumes 100% CO2)' write(*,*) 'co2ice=0 : remove CO2 polar cap' write(*,*) 'ptot : change total pressure' write(*,*) 'emis : change surface emissivity' write(*,*) 'therm_ini_s: Set soil thermal inertia to reference sur &face values' write(*,*) 'inert3d: give uniform value of thermal inertial' write(*,*) 'subsoilice_n: Put deep underground ice layer in northe &rn hemisphere' write(*,*) 'subsoilice_s: Put deep underground ice layer in southe &rn hemisphere' write(*,*) 'mons_ice: Put underground ice layer according to MONS- &derived data' write(*,*) 'plutomap: initialize surface like pluto from ...' write(*,*) 'subsoil_n2: Put deep underground n2 ice layer where & there is already n2 ice' write(*,*) 'initsurf: initial surface temperature' write(*,*) write(*,*) 'Change to perform ?' write(*,*) ' (enter keyword or return to end)' write(*,*) read(*,fmt='(a20)') modif if (modif(1:1) .eq. ' ') exit ! exit loop on changes write(*,*) write(*,*) trim(modif) , ' : ' c 'flat : no topography ("aquaplanet")' c ------------------------------------- if (modif(1:len_trim(modif)) .eq. 'flat') then ! set topo to zero CALL initial0(ip1jmp1,z_reel) CALL multscal(ip1jmp1,z_reel,g,phis) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) write(*,*) 'topography set to zero.' c Choice for surface pressure yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*) 'Do you wish to choose homogeneous surface', & 'pressure (y) or let newstart interpolate ', & ' the previous field (n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) 'New value for ps (Pa) ?' 201 read(*,*,iostat=ierr) patm if(ierr.ne.0) goto 201 write(*,*) write(*,*) ' new ps everywhere (Pa) = ', patm write(*,*) do j=1,jjp1 do i=1,iip1 ps(i,j)=patm enddo enddo end if c bilball : albedo, inertie thermique uniforme c -------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'bilball') then write(*,*) 'constante albedo and iner.therm:' write(*,*) 'New value for albedo (ex: 0.25) ?' 101 read(*,*,iostat=ierr) alb_bb if(ierr.ne.0) goto 101 write(*,*) ' uniform albedo (new value):',alb_bb write(*,*) 'New value for thermal inertia (eg: 247) ?' 150 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 150 write(*,*) 'uniform thermal inertia (new value):',ith_bb DO ig=1,ngridmx DO j=1,nsoilmx if(qsurf(ig,igcm_ch4_ice).gt.0)then ithfi(ig,j)=ith_bb else ithfi(ig,j)=400 endif ENDDO ENDDO DO j=1,jjp1 DO i=1,iip1 alb(i,j) = alb_bb ! albedo do isoil=1,nsoilmx ith(i,j,isoil) = ith_bb ! thermal inertia enddo END DO END DO CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi) CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) c TB15 : subsoil_n2 : choice of thermal inertia values for N2 only c ---------------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'subsoil_n2') then write(*,*) 'New value for subsoil thermal inertia ?' 102 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 102 write(*,*) 'thermal inertia (new value):',ith_bb write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(1),' - ',layer(nsoilmx),')' write(*,*)'write 0 for uniform value for all subsurf levels?' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val2 write(*,*)'val2:',val2,'ierr=',ierr if (ierr.eq.0) then ! got a value, but do a sanity check if(val2.gt.layer(nsoilmx)) then write(*,*)'Depth should be less than ',layer(nsoilmx) ierr=1 endif if(val2.lt.layer(1)) then if(val2.eq.0) then write(*,*)'Thermal inertia value set for all subsurface layers' ierr=0 else write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif endif enddo ! of do while ! find the reference index iref the depth corresponds to if(val2.eq.0) then iref=1 write(*,*)'Level selected is first level: ',layer(iref),' m' else do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil write(*,*)'Level selected : ',layer(isoil),' m' exit endif enddo endif DO ig=1,ngridmx DO j=iref,nsoilmx if((qsurf(ig,igcm_ch4_ice).lt.1.).and. & (qsurf(ig,igcm_co_ice).lt.1.))then ithfi(ig,j)=ith_bb endif ENDDO ENDDO CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) c TB15 Choice tsurf c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'initsurf') then write(*,*) 'New value for initial surface temperature ?' 105 read(*,*,iostat=ierr) tsurf_bb if(ierr.ne.0) goto 105 write(*,*) 'new surface temperature (K):',tsurf_bb DO ig=1,ngridmx tsurf(ig)=tsurf_bb ENDDO c coldspole : sous-sol de la calotte sud toujours froid c ----------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'coldspole') then write(*,*)'Sorry, not available' c ptot : Modification of the total pressure: ice + current atmosphere c ------------------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'ptot') then c calcul de la pression totale glace + atm actuelle patm=0. airetot=0. pcap=0. DO j=1,jjp1 DO i=1,iim ig=1+(j-2)*iim +i if(j.eq.1) ig=1 if(j.eq.jjp1) ig=ngridmx patm = patm + ps(i,j)*aire(i,j) airetot= airetot + aire(i,j) ENDDO ENDDO ptoto = pcap + patm print*,'Current total pressure at surface ', & ptoto/airetot print*,'new value?' read(*,*) ptotn ptotn=ptotn*airetot patmn=ptotn-pcap print*,'ptoto,patm,ptotn,patmn' print*,ptoto,patm,ptotn,patmn print*,'Mult. factor for pressure (atm only)', patmn/patm do j=1,jjp1 do i=1,iip1 ps(i,j)=ps(i,j)*patmn/patm enddo enddo c Correction pour la conservation des traceurs yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*) 'Do you wish to conserve tracer total mass (y)', & ' or tracer mixing ratio (n) ?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) 'OK : conservation of tracer total mass' DO iq =1, nqmx DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn ENDDO ENDDO ENDDO ENDDO else write(*,*) 'OK : conservation of tracer mixing ratio' end if c qname : change tracer name c -------------------------- else if (trim(modif).eq.'qname') then yes='y' do while (yes.eq.'y') write(*,*) 'Which tracer name do you want to change ?' do iq=1,nqmx write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq)) enddo write(*,'(a35,i3)') & '(enter tracer number; between 1 and ',nqmx write(*,*)' or any other value to quit this option)' read(*,*) iq if ((iq.ge.1).and.(iq.le.nqmx)) then write(*,*)'Change tracer name ',trim(tnom(iq)),' to ?' read(*,*) txt tnom(iq)=txt write(*,*)'Do you want to change another tracer name (y/n)?' read(*,'(a)') yes else ! inapropiate value of iq; quit this option yes='n' endif ! of if ((iq.ge.1).and.(iq.le.nqmx)) enddo ! of do while (yes.ne.'y') c q=0 : set tracers to zero c ------------------------- else if (modif(1:len_trim(modif)).eq.'q=0') then c mise a 0 des q (traceurs) write(*,*) 'Tracers set to 0 (1.E-30 in fact)' DO iq =1, nqmx DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=1.e-30 ENDDO ENDDO ENDDO ENDDO c set surface tracers to zero DO iq =1, nqmx DO ig=1,ngridmx qsurf(ig,iq)=0. ENDDO ENDDO c q=x : initialise tracer manually c -------------------------------- else if (modif(1:len_trim(modif)).eq.'q=x') then write(*,*) 'Which tracer do you want to modify ?' do iq=1,nqmx write(*,*)iq,' : ',trim(tnom(iq)),' : ',q(1,1,1,iq) enddo write(*,*) '(choose between 1 and ',nqmx,')' read(*,*) iq write(*,*)'mixing ratio of tracer ',trim(tnom(iq)), & ' ? (kg/kg)' read(*,*) val DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,iq)=val ENDDO ENDDO ENDDO yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'Do you want to change surface value (y/n)?' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) 'SURFACE value of tracer ',trim(tnom(iq)), & ' ? (kg/m2)' read(*,*) val DO ig=1,ngridmx qsurf(ig,iq)= val ENDDO endif c DO ig=1,ngridmx c if (ig.lt.ngridmx*11/32+1) then c qsurf(ig,iq)=300 c else if (ig.gt.ngridmx*26/32+1) then c qsurf(ig,iq)=val c else c qsurf(ig,iq)=0 c endif c ENDDO c endif c ini_q : Initialize tracers for chemistry c ----------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'ini_q') then write(*,*) 'Sorry, not available' c ini_q-H2O : as above exept for the water vapour tracer c ------------------------------------------------------ else if (modif(1:len_trim(modif)) .eq. 'ini_q-H2O') then write(*,*) 'Sorry, not available' c ini_q-iceH2O : as above exept for ice et H2O c ----------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'ini_q-iceH2O') then write(*,*) 'Sorry, not available' c wetstart : wet atmosphere with a north to south gradient c -------------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'wetstart') then write(*,*) "No water vapour tracer! Can't use this option" c noglacier : remove tropical water ice (to initialize high res sim) c -------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'noglacier') then write(*,*) "No water vapour tracer! Can't use this option" c watercapn : H20 ice on permanent northern cap c -------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'watercapn') then write(*,*) "No water vapour tracer! Can't use this option" c watercaps : H20 ice on permanent southern cap c ------------------------------------------------- else if (modif(1:len_trim(modif)) .eq. 'watercaps') then write(*,*) "No water vapour tracer! Can't use this option" c TB14 surface emissivity c -------------------- else if (modif(1:len_trim(modif)) .eq. 'emis') then write(*,*) 'Emissivity actual max/min values:' & ,minval(emis),maxval(emis) write(*,*) 'Values are in set in physiq !' write(*,*) 'Value emissivity N2' read(*,*) val write(*,*) 'Value emissivity ch4 ' read(*,*) val2 write(*,*) 'Value emissivity tholins ' read(*,*) val3 write(*,*) 'Value emissivity co ' read(*,*) val4 write(*,*) 'New emissivity values :' write(*,*) 'N2= ',val,' CH4= ',val2 write(*,*)' tholins= ',val3,'co= ',val4 DO ig=1,ngridmx if ((qsurf(ig,igcm_n2).eq.0).and. & (qsurf(ig,igcm_ch4_ice).eq.0).and. & (qsurf(ig,igcm_co_ice).eq.0)) then emis(ig)=val3 else if ((qsurf(ig,igcm_n2).eq.0).and. & (qsurf(ig,igcm_ch4_ice).eq.0)) then emis(ig)=val4 else if ((qsurf(ig,igcm_n2).eq.0).and. & (qsurf(ig,igcm_co_ice).eq.0)) then emis(ig)=val2 else emis(ig)=val endif ENDDO c isothermal temperatures and no winds c ------------------------------------------------ else if (modif(1:len_trim(modif)) .eq. 'isotherm') then write(*,*)'Sorry, not available' c initsoil: subsurface temperature c --------------------------------- else if (modif(1:len_trim(modif)) .eq. 'initsoil') then write(*,*)'Isothermal temperature of the subsurface' write(*,*) 'Value of this temperature ? :' 303 read(*,*,iostat=ierr) Tiso if(ierr.ne.0) goto 303 do l=1,nsoilmx do ig=1, ngridmx tsoil(ig,l) = Tiso end do end do c ----------------------------------------------------------------------- c coldstart : T set 1K above CO2 frost point and no winds c ------------------------------------------------ else if (modif(1:len_trim(modif)) .eq. 'coldstart') then write(*,*)'Sorry, not available : No CO2 tracer' c n2ice=0 : remove n2 polar ice caps' c ------------------------------------------------ else if (modif(1:len_trim(modif)) .eq. 'N2ice=0') then ! check that there is indeed a N2_ice tracer ... if (igcm_n2.ne.0) then do ig=1,ngridmx !co2ice(ig)=0 qsurf(ig,igcm_n2)=0 emis(ig)=emis(ngridmx/2) end do else write(*,*) "Can't remove N2 ice!! (no N2 tracer)" endif ! therm_ini_s: (re)-set soil thermal inertia to reference surface values ! ---------------------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'therm_ini_s') then ! write(*,*)"surfithfi(1):",surfithfi(1) do isoil=1,nsoilmx inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx) enddo write(*,*)'OK: Soil thermal inertia has been reset to referenc &e surface values' write(*,*)"inertiedat(1,1):",inertiedat(1,1) ithfi(:,:)=inertiedat(:,:) ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) ! Check: ! do i=1,iip1 ! do j=1,jjp1 ! do isoil=1,nsoilmx ! write(77,*) i,j,isoil," ",ith(i,j,isoil) ! enddo ! enddo ! enddo ! TB14 inert3d: set soil thermal inertia to specific uniform value ! ---------------------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'inert3d') then write(*,*) 'Actual value of surf Thermal inertia at ig=1: ', & inertiedat(1,1), ' SI' write(*,*) 'Value of Thermal inertia (SI) ? ' read(*,*) val do isoil=1,nsoilmx do ig=1,ngridmx inertiedat(ig,isoil)=val enddo enddo write(*,*)'OK: Soil TI set to ',val,' SI everywhere' ithfi(:,:)=inertiedat(:,:) ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) ! subsoilice_n: Put deep ice layer in northern hemisphere soil ! ------------------------------------------------------------ else if (modif(1:len_trim(modif)).eq.'subsoilice_n') then write(*,*)'From which latitude (in deg.), up to the north pole, &should we put subterranean ice?' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val if (ierr.eq.0) then ! got a value ! do a sanity check if((val.lt.0.).or.(val.gt.90)) then write(*,*)'Latitude should be between 0 and 90 deg. !!!' ierr=1 else ! find corresponding jref (nearest latitude) ! note: rlatu(:) contains decreasing values of latitude ! starting from PI/2 to -PI/2 do j=1,jjp1 if ((rlatu(j)*180./pi.ge.val).and. & (rlatu(j+1)*180./pi.le.val)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val).le. & abs((rlatu(j+1)*180./pi-val))) then jref=j else jref=j+1 endif write(*,*)'Will use nearest grid latitude which is:', & rlatu(jref)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif !of if (ierr.eq.0) enddo ! of do while ! Build layers() (as in soil_settings.F) val2=sqrt(mlayer(0)*mlayer(1)) val3=mlayer(1)/mlayer(0) do isoil=1,nsoilmx layer(isoil)=val2*(val3**(isoil-1)) enddo write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(nsoilmx),')' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val2 ! write(*,*)'val2:',val2,'ierr=',ierr if (ierr.eq.0) then ! got a value, but do a sanity check if(val2.gt.layer(nsoilmx)) then write(*,*)'Depth should be less than ',layer(nsoilmx) ierr=1 endif if(val2.lt.layer(1)) then write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif enddo ! of do while ! find the reference index iref the depth corresponds to ! if (val2.lt.layer(1)) then ! iref=1 ! else do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil exit endif enddo ! endif ! write(*,*)'iref:',iref,' jref:',jref ! write(*,*)'layer',layer ! write(*,*)'mlayer',mlayer ! thermal inertia of the ice: ierr=1 do while (ierr.ne.0) write(*,*)'What is the value of subterranean ice thermal inert &ia? (e.g.: 2000)' read(*,*,iostat=ierr)iceith enddo ! of do while ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) do j=1,jref ! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi do i=1,iip1 ! loop on longitudes ! Build "equivalent" thermal inertia for the mixed layer ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ & (((val2-layer(iref))/(ith(i,j,iref)**2))+ & ((layer(iref+1)-val2)/(iceith)**2))) ! Set thermal inertia of lower layers do isoil=iref+2,nsoilmx ith(i,j,isoil)=iceith ! ice enddo enddo ! of do i=1,iip1 enddo ! of do j=1,jjp1 CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) ! do i=1,nsoilmx ! write(*,*)'i:',i,'ithfi(1,i):',ithfi(1,i) ! enddo ! subsoilice_s: Put deep ice layer in southern hemisphere soil ! ------------------------------------------------------------ else if (modif(1:len_trim(modif)).eq.'subsoilice_s') then write(*,*)'From which latitude (in deg.), down to the south pol &e, should we put subterranean ice?' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val if (ierr.eq.0) then ! got a value ! do a sanity check if((val.gt.0.).or.(val.lt.-90)) then write(*,*)'Latitude should be between 0 and -90 deg. !!!' ierr=1 else ! find corresponding jref (nearest latitude) ! note: rlatu(:) contains decreasing values of latitude ! starting from PI/2 to -PI/2 do j=1,jjp1 if ((rlatu(j)*180./pi.ge.val).and. & (rlatu(j+1)*180./pi.le.val)) then ! find which grid point is nearest to val: if (abs(rlatu(j)*180./pi-val).le. & abs((rlatu(j+1)*180./pi-val))) then jref=j else jref=j+1 endif write(*,*)'Will use nearest grid latitude which is:', & rlatu(jref)*180./pi endif enddo ! of do j=1,jjp1 endif ! of if((val.lt.0.).or.(val.gt.90)) endif !of if (ierr.eq.0) enddo ! of do while ! Build layers() (as in soil_settings.F) val2=sqrt(mlayer(0)*mlayer(1)) val3=mlayer(1)/mlayer(0) do isoil=1,nsoilmx layer(isoil)=val2*(val3**(isoil-1)) enddo write(*,*)'At which depth (in m.) does the ice layer begin?' write(*,*)'(currently, the deepest soil layer extends down to:' & ,layer(nsoilmx),')' ierr=1 do while (ierr.ne.0) read(*,*,iostat=ierr) val2 ! write(*,*)'val2:',val2,'ierr=',ierr if (ierr.eq.0) then ! got a value, but do a sanity check if(val2.gt.layer(nsoilmx)) then write(*,*)'Depth should be less than ',layer(nsoilmx) ierr=1 endif if(val2.lt.layer(1)) then write(*,*)'Depth should be more than ',layer(1) ierr=1 endif endif enddo ! of do while ! find the reference index iref the depth corresponds to do isoil=1,nsoilmx-1 if((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1))) & then iref=isoil exit endif enddo ! thermal inertia of the ice: ierr=1 do while (ierr.ne.0) write(*,*)'What is the value of subterranean ice thermal inert &ia? (e.g.: 2000)' read(*,*,iostat=ierr)iceith enddo ! of do while ! recast ithfi() onto ith() call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) do j=jref,jjp1 ! write(*,*)'j:',j,'rlatu(j)*180./pi:',rlatu(j)*180./pi do i=1,iip1 ! loop on longitudes ! Build "equivalent" thermal inertia for the mixed layer ith(i,j,iref+1)=sqrt((layer(iref+1)-layer(iref))/ & (((val2-layer(iref))/(ith(i,j,iref)**2))+ & ((layer(iref+1)-val2)/(iceith)**2))) ! Set thermal inertia of lower layers do isoil=iref+2,nsoilmx ith(i,j,isoil)=iceith ! ice enddo enddo ! of do i=1,iip1 enddo ! of do j=jref,jjp1 CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) c 'mons_ice' : use MONS data to build subsurface ice table c -------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'mons_ice') then write(*,*) "Sorry, not available" c -------------------------------------------------------- c 'plutomap' : initialize pluto c -------------------------------------------------------- else if (modif(1:len_trim(modif)).eq.'plutomap') then write(*,*) 'Pluto_map.dat has to be in your simulation file !!' write(*,*) 'Initial value of CH4 ice (kg)' read(*,*,iostat=ierr) val write(*,*) 'Initial value of N2 ice (kg)' read(*,*,iostat=ierr) val2 write(*,*) 'Initial value of CO ice (kg)' read(*,*,iostat=ierr) val3 c Reading map codes (from Lellouch et al. 2011) c The map must be placed in your simulation file open(27,file='pluto_map.dat') do j=1,jm_plu+1 read(27,*) (map_pluto_dat(i,j),i=1,im_plu+1) do i=1,im_plu+1 N2_pluto_dat(i,j) =0. CH4_pluto_dat(i,j) =0. CO_pluto_dat(i,j) =0. if(map_pluto_dat(i,j).eq.1) N2_pluto_dat(i,j)=val2 if(map_pluto_dat(i,j).eq.2) CH4_pluto_dat(i,j)=val c if(map_pluto_dat(i,j).eq.3) -> nothing : Tholins if(map_pluto_dat(i,j).eq.4) CO_pluto_dat(i,j)=val3 c if(map_pluto_dat(i,j).eq.1) alb_pluto_dat(i,j)=0.657 c if(map_pluto_dat(i,j).eq.1) ith_pluto_dat(i,j)=100 end do end do close (27) c Option to remove N2 ice in the south as in Young/Olkin model yes=' ' do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'No nitrogen in the south? as in Olkin et al. 2014' read(*,fmt='(a)') yes end do if (yes.eq.'y') then write(*,*) 'Lower latitude of nitrogen (deg: -90 to 90)' read(*,*) val itau = nint(jm_plu*(90.-val)/180.)+1 write(*,*) "grid lower latitude=" write(*,*)90.-180.*(itau-1)/(jm_plu) do j=itau, jm_plu+1 do i=1,im_plu+1 if(map_pluto_dat(i,j).eq.1) N2_pluto_dat(i,j)=0. end do end do end if c Reindexation to shift the longitude grid like a LMD GCM grid... do j=1,jm_plu+1 N2_pluto_rein(1,j)=(N2_pluto_dat(1,j)+N2_pluto_dat(im_plu,j))/2 CH4_pluto_rein(1,j)=(CH4_pluto_dat(1,j)+CH4_pluto_dat(im_plu,j))/2 CO_pluto_rein(1,j)=(CO_pluto_dat(1,j)+CO_pluto_dat(im_plu,j))/2 c alb_pluto_rein(1,j)=(alb_pluto_dat(1,j)+alb_pluto_dat(im_plu,j))/2 c ith_pluto_rein(1,j)=(ith_pluto_dat(1,j)+ith_pluto_dat(im_plu,j))/2 do i=2,im_plu N2_pluto_rein(i,j)=(N2_pluto_dat(i-1,j)+N2_pluto_dat(i,j))/2 CH4_pluto_rein(i,j)=(CH4_pluto_dat(i-1,j)+CH4_pluto_dat(i,j))/2 CO_pluto_rein(i,j)=(CO_pluto_dat(i-1,j)+CO_pluto_dat(i,j))/2 c alb_pluto_rein(i,j)=(alb_pluto_dat(i-1,j)+alb_pluto_dat(i,j))/2 c ith_pluto_rein(i,j)=(ith_pluto_dat(i-1,j)+ith_pluto_dat(i,j))/2 end do N2_pluto_rein(im_plu+1,j) = N2_pluto_rein(1,j) CH4_pluto_rein(im_plu+1,j) = CH4_pluto_rein(1,j) CO_pluto_rein(im_plu+1,j) = CO_pluto_rein(1,j) c alb_pluto_rein(im_plu+1,j) = alb_pluto_rein(1,j) c ith_pluto_rein(im_plu+1,j)= ith_pluto_dat(1,j) end do c latitude and longitude in REindexed pluto map : latv_plu(1)=90. *pi/180. do j=2,jm_plu latv_plu(j)= (90-(j-1 -0.5)*180./(jm_plu-1))*pi/180. end do latv_plu(jm_plu+1)= -90. *pi/180. do i=1,im_plu+1 lonu_plu(i)=(-180+ (i-1)*360./(im_plu) ) *pi/180. enddo c write(*,*) 'TB15 im_plu = ',im_plu c write(*,*) 'TB15 jm_plu = ',jm_plu c write(*,*) 'TB15 imn = ',iim c write(*,*) 'TB15 jmn = ',jjm c write(*,*) 'TB15 lonu_plu = ',lonu_plu*180./pi c write(*,*) 'TB15 latv_plu = ',latv_plu*180./pi c write(*,*) 'TB15 rlonu = ',rlonu*180./pi c write(*,*) 'TB15 rlatv = ',rlatv*180./pi c interpolation in LMD GCM grid: c write(*,*) 'TB15 N2_pluto_rein = ',N2_pluto_rein call interp_horiz(N2_pluto_rein,N2_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(CH4_pluto_rein,CH4_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) call interp_horiz(CO_pluto_rein,CO_pluto_gcm, & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) c call interp_horiz(alb_pluto_rein,alb_pluto_gcm, c & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) c write(*,*) 'TB15 N2_pluto_gcm = ',N2_pluto_gcm c do i=1,iip1 c do j=1,jjp1 c N2_pluto_gcm(i,j) = 1000. c end do c end do c write(*,*) 'TB15 N2_pluto_gcm = ',N2_pluto_gcm c Option to remove N2 ice in the North c yes=' ' c do while ((yes.ne.'y').and.(yes.ne.'n')) c write(*,*)'No N2 ice at North pole ?' c read(*,fmt='(a)') yes c end do c if (yes.eq.'y') then c do i=1,iip1 c N2_pluto_gcm(i,1) = 0. c end do c end if cc call interp_horiz(ith_pluto_rein,ith_pluto_gcm, cc & im_plu,jm_plu,iim,jjm,1,lonu_plu,latv_plu,rlonu,rlatv) qsurf(:,:) =0. c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb_pluto_gcm,albfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & N2_pluto_gcm,qsurf(1,igcm_n2)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & CH4_pluto_gcm,qsurf(1,igcm_ch4_ice)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx, & CO_pluto_gcm,qsurf(1,igcm_co_ice)) c ---------------------------------------------------------------- c ---------------------------------------------------------------- else write(*,*) ' Unknown (misspelled?) option!!!' end if ! of if (modif(1:len_trim(modif)) .eq. '...') elseif ... enddo ! of do ! infinite loop on liste of changes 999 continue c======================================================================= c Correct pressure on the new grid (menu 0) c======================================================================= if ((choix_1.eq.0).and.(.not.flagps0)) then r = 8.314511E+0 *1000.E+0/mugaz do j=1,jjp1 do i=1,iip1 ps(i,j) = ps(i,j) * . exp((phisold_newgrid(i,j)-phis(i,j)) / . (t(i,j,1) * r)) end do end do c periodicity of surface ps in longitude do j=1,jjp1 ps(1,j) = ps(iip1,j) end do end if c======================================================================= c======================================================================= c======================================================================= c Initialisation de la physique / ecriture de newstartfi : c======================================================================= CALL inifilr CALL pression(ip1jmp1, ap, bp, ps, p3d) c----------------------------------------------------------------------- c Initialisation pks: c----------------------------------------------------------------------- CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf) ! Calcul de la temperature potentielle teta if (flagtset) then DO l=1,llm DO j=1,jjp1 DO i=1,iim teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l) ENDDO teta (iip1,j,l)= teta (1,j,l) ENDDO ENDDO else if (choix_1.eq.0) then DO l=1,llm DO j=1,jjp1 DO i=1,iim teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l) ENDDO teta (iip1,j,l)= teta (1,j,l) ENDDO ENDDO endif C Calcul intermediaire if (choix_1.eq.0) then CALL massdair( p3d, masse ) print *,' ALPHAX ',alphax DO l = 1, llm DO i = 1, iim xppn(i) = aire( i, 1 ) * masse( i , 1 , l ) xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l ) ENDDO xpn = SUM(xppn)/apoln xps = SUM(xpps)/apols DO i = 1, iip1 masse( i , 1 , l ) = xpn masse( i , jjp1 , l ) = xps ENDDO ENDDO endif phis(iip1,:) = phis(1,:) CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh, * tetagdiv, tetagrot , tetatemp ) itau=0 if (choix_1.eq.0) then day_ini=int(date) endif CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis , * phi,w, pbaru,pbarv,day_ini+time ) c CALL caldyn c $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis , c $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, day_ini ) CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx) CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps) C Ecriture etat initial physique call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx, . dtphys,float(day_ini), . time,tsurf,tsoil,emis,q2,qsurf, . airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe) write(*,*) '***************************************' write(*,*) 'FINAL VALUES IN RESTART: ' write(*,*) 'check consistency' write(*,*) 'to change values in restart use start2archive' write(*,*) '***************************************' write(*,*) 'radius= ',rad write(*,*) 'omeg= ',omeg write(*,*) 'gravity= ',g write(*,*) 'cpp= ', cpp write(*,*) 'mugaz= ',mugaz write(*,*) 'daysec=',daysec write(*,*) 'kappa=', 8.314511E+0 *1000.E+0/(mugaz*cpp) c======================================================================= c Formats c======================================================================= 1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema *rrage est differente de la valeur parametree iim =',i4//) 2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema *rrage est differente de la valeur parametree jjm =',i4//) 3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar *rage est differente de la valeur parametree llm =',i4//) end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c subroutine load_MONS_data(MONS_Hdn,MONS_d21)