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 "dimradmars.h" #include "yomaer.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" c======================================================================= c Declarations c======================================================================= c Variables dimension du fichier "start_archive" c------------------------------------ CHARACTER relief*3 c et autres: c---------- INTEGER lnblnk EXTERNAL lnblnk c Variables pour les lectures NetCDF des fichiers "start_archive" c-------------------------------------------------- INTEGER nid_dyn, nid_fi,nid,nvarid INTEGER size INTEGER length parameter (length = 100) INTEGER tab0 INTEGER NB_ETATMAX parameter (NB_ETATMAX = 100) REAL date REAL p_rad,p_omeg,p_g,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),tsoil(ngridmx,nsoilmx),co2ice(ngridmx) REAL emis(ngridmx),qsurf(ngridmx,nqmx) REAL q2(ngridmx,nlayermx+1) REAL rnaturfi(ngridmx) real alb(iip1,jjp1),albfi(ngridmx) real ith(iip1,jjp1),ithfi(ngridmx) REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) INTEGER i,j,l,idum,ig real mugaz ! masse molaire de l''atmosphere EXTERNAL iniconst,geopot,inigeom integer ierr, nbetat integer ismin external ismin c Variable nouvelle grille naturelle au point scalaire c------------------------------------------------------ REAL p(iip1,jjp1) REAL t(iip1,jjp1,llm) real phisold_newgrid(iip1,jjp1) REAL :: teta(iip1, jjp1, llm) REAL :: pk(iip1,jjp1,llm) REAL :: pkf(iip1,jjp1,llm) REAL :: ps(iip1, jjp1) 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 character*80 fichnom integer Lmodif,iq,thermo character modif*20 real z_reel(iip1,jjp1) real tsud,albsud,alb_bb,ith_bb,Tiso real ptoto,pcap,patm,airetot,ptotn,patmn real ssum character*1 yes logical :: flagiso=.false. , flagps0=.false. real val INTEGER :: itau c sortie visu pour les champs dynamiques c--------------------------------------- INTEGER :: visuid real :: time_step,t_ops,t_wrt CHARACTER*80 :: visu_file cpp = 744.499 ! au lieu de 1004.70885 (Terre) preff = 610. ! au lieu de 101325. (Terre) pa= 20 ! au lieu de 500 (Terre) c======================================================================= c Choix du fichier de demarrage a modifier 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 Ouverture des fichiers a modifier (start ou start_archive) c----------------------------------------------------------------------- DO read(*,*,iostat=ierr) choix_1 if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT ENDDO c Ouverture de start_archive c ~~~~~~~~~~~~~~~~~~~~~~~~~~ if (choix_1.eq.0) then write(*,*) 'Creating start files from:' write(*,*) './start_archive.dat and ./start_archive.dic' write(*,*) fichnom = 'start_archive.nc' ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) IF (ierr.NE.NF_NOERR) THEN write(6,*)' Problem opening: ',fichnom write(6,*)' ierr = ', ierr CALL ABORT ENDIF tab0 = 50 Lmodif = 1 c OU Ouverture de start 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: ',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: ',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 write(*,*) 'printing tab_cntrl', tab_cntrl do i=1,100 write(*,*) i,tab_cntrl(i) enddo endif 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(*,*) "kappa,etot0,ptot0,ztot0,stot0,ang0" write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0 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_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_mugaz,p_daysec,time) endif rad = p_rad omeg = p_omeg g = p_g mugaz = p_mugaz daysec = p_daysec write(*,*) 'aire',aire c======================================================================= c INITIALISATIONS DIVERSES c======================================================================= day_step=180 CALL defrun_new( 99, .TRUE. ) CALL iniconst CALL inigeom idum=-1 idum=0 c Initialisation coordonnees /aires c ------------------------------- 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 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,co2ice,teta,phisold_newgrid,q,qsurf,nid) 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,co2ice) do i=1,ngridmx albfi(i) = albedodat(i) ithfi(i) = inertiedat(i) enddo c save old topography do j=1,jjp1 do i=1,iip1 phisold_newgrid(i,j)=phis(i,j) enddo enddo else CALL exit(1) endif dtvr = daysec/FLOAT(day_step) dtphys = dtvr * FLOAT(iphysiq) 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 relief="mola" c enddo CALL datareadnc(relief,phis,alb,ith,zmeaS,zstdS,zsigS,zgamS, . ztheS) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) 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 c======================================================================= c c======================================================================= 888 continue write(*,*) write(*,*) write(*,*) 'Other 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(*,*) '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(*,*) 'watercapn : H20 ice on permanent N polar cap ' write(*,*) 'watercaps : H20 ice on permanent S polar cap ' write(*,*) 'wetstart : start with a wet atmosphere' write(*,*) 'isotherm : Isothermal Temperatures, wind set to zero' write(*,*) 'co2ice=0 : remove CO2 polar cap' write(*,*) 'ptot : change total pressure' do while(modif(1:1).ne.'hello') write(*,*) write(*,*) 'Changes to perform ?' write(*,*) ' (enter keyword or return to end)' write(*,*) read(*,fmt='(a20)') modif if (modif(1:1) .eq. ' ') goto 999 write(*,*) write(*,*) modif(1:lnblnk(modif)) , ' : ' c 'flat : no topography ("aquaplanet")' c ------------------------------------- if (modif(1:lnblnk(modif)) .eq. 'flat') then c 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.' write(*,*) 'WARNING : the subgrid topography parameters', & ' were not set to zero ! => set calllott to F' 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 flagps0=.true. 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:lnblnk(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(*,*) write(*,*) ' albedo uniform (new value):',alb_bb write(*,*) write(*,*) 'New value for iner.therm (ex: 247) ?' 102 read(*,*,iostat=ierr) ith_bb if(ierr.ne.0) goto 102 write(*,*) 'iner.therm uniform (new value):',ith_bb DO j=1,jjp1 DO i=1,iip1 alb(i,j) = alb_bb ith(i,j) = ith_bb END DO END DO CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) c coldspole : sous-sol de la calotte sud toujours froid c ----------------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'coldspole') then write(*,*)'nouvelle valeur de la temperature du', & 'sous sol de la calotte permanente sud ? (ex: 141 K)' 103 read(*,*,iostat=ierr) tsud if(ierr.ne.0) goto 103 write(*,*) write(*,*) ' nouvelle valeur de la temperature:',tsud c nouvelle temperature sous la calotte permanente do l=2,nsoilmx tsoil(ngridmx,l) = tsud end do write(*,*)'nouvelle valeur de l albedo', & 'de la calotte permanente sud ? (ex: 0.75 K)' 104 read(*,*,iostat=ierr) albsud if(ierr.ne.0) goto 104 write(*,*) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c Option 1: Seul l'albedo du pole est modifié : albfi(ngridmx)=albsud write(*,*) 'ig=',ngridmx,' albedo perennial cap ', & albfi(ngridmx) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~é c Option 2 A haute resolution : coordonnée de la vrai calotte ~ c DO j=1,jjp1 c DO i=1,iip1 c ig=1+(j-2)*iim +i c if(j.eq.1) ig=1 c if(j.eq.jjp1) ig=ngridmx c if ((rlatu(j)*180./pi.lt.-84.).and. c & (rlatu(j)*180./pi.gt.-91.).and. c & (rlonv(i)*180./pi.gt.-91.).and. c & (rlonv(i)*180./pi.lt.0.)) then cc albedo de la calotte permanente fixe a albsud c alb(i,j)=albsud c write(*,*) 'lat=',rlatu(j)*180./pi, c & ' lon=',rlonv(i)*180./pi cc fin de la condition sur les limites de la calotte permanente c end if c ENDDO c ENDDO c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) c ptot : Modification de la pression totale glace + atm actuelle c -------------------------------------------------------------- else if (modif(1:lnblnk(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) pcap = pcap + aire(i,j)*co2ice(ig)*g ENDDO ENDDO ptoto = pcap + patm print*,'Pression totale au sol actuelle (co2 ice + atm) ', & ptoto/airetot print*,'nouvelle valeur?' read(*,*) ptotn ptotn=ptotn*airetot patmn=ptotn-pcap print*,'ptoto,patm,ptotn,patmn' print*,ptoto,patm,ptotn,patmn print*,'Facteur mult de la pression (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 q=0 : traceurs a zero c --------------------- else if (modif(1:lnblnk(modif)) .eq. 'q=0') then c mise a 0 des q (traceurs) write(*,*) 'Traceurs mis a 0 (1.E-30 en fait)' 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 mise a 0 des qsurf (traceurs a la surface) 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:lnblnk(modif)) .eq. 'q=x') then write(*,*) 'Which tracer do you want to modify ?' write(*,*) '(choose between 1 and ',nqmx,')' read(*,*) iq write(*,*)'mixing ratio of tracer #',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 write(*,*) 'SURFACE value of tracer #', iq , ' ? (kg/m2)' read(*,*) val DO ig=1,ngridmx qsurf(ig,iq)=val ENDDO c ini_q : traceurs initialises pour la chimie c ----------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'ini_q') then c For more than 32 layers, possible to initiate thermosphere only thermo=0 yes=' ' if (llm.gt.32) then do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'', & 'initialisation for thermosphere only? (y/n)' read(*,fmt='(a)') yes if (yes.eq.'y') then thermo=1 else thermo=0 endif enddo endif call inichim_newstart(q,ps,sig,nqmx,latfi,lonfi,airefi, $ thermo) write(*,*) 'Especes chimiques initialisees' if (thermo.eq.0) then c mise a 0 des qsurf (traceurs a la surface) DO iq =1, nqmx DO ig=1,ngridmx qsurf(ig,iq)=0. ENDDO ENDDO endif c ini_q-H2O : idem sauf le dernier traceur (H2O) c ----------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'ini_q-H2O') then ! for more than 32 layers, possible to initiate thermosphere only thermo=0 yes=' ' if(llm.gt.32) then do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'', & 'initialisation for thermosphere only? (y/n)' read(*,fmt='(a)') yes if (yes.eq.'y') then thermo=1 else thermo=0 endif enddo endif call inichim_newstart(q,ps,sig,nqmx-1,latfi,lonfi,airefi, $ thermo) write(*,*) 'Esp. chim. initialisees sauf derniere (H2O)' if (thermo.eq.0) then c mise a 0 des qsurf (traceurs a la surface) DO iq =1, nqmx-1 DO ig=1,ngridmx qsurf(ig,iq)=0. ENDDO ENDDO endif c ini_q-iceH2O : idem sauf ice et H2O c ----------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'ini_q-iceH2O') then c For more than 32 layers, possible to initiate thermosphere only thermo=0 yes=' ' if(llm.gt.32) then do while ((yes.ne.'y').and.(yes.ne.'n')) write(*,*)'', & 'initialisation for thermosphere only? (y/n)' read(*,fmt='(a)') yes if (yes.eq.'y') then thermo=1 else thermo=0 endif enddo endif call inichim_newstart(q,ps,sig,nqmx-2,latfi,lonfi,airefi, $ thermo) write(*,*) 'Esp. chim. initialisees sauf ice et H2O' if (thermo.eq.0) then c mise a 0 des qsurf (traceurs a la surface) DO iq =1, nqmx-2 DO ig=1,ngridmx qsurf(ig,iq)=0. ENDDO ENDDO endif c wetstart : wet atmosphere with a north to south gradient c -------------------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'wetstart') then DO l=1,llm DO j=1,jjp1 DO i=1,iip1 q(i,j,l,nqmx)=150.e-6 * (rlatu(j)+pi/2.) / pi ENDDO ENDDO ENDDO write(*,*) 'Water mass mixing ratio at north pole=' * ,q(1,1,1,nqmx) write(*,*) '---------------------------south pole=' * ,q(1,jjp1,1,nqmx) c noglacier : remove tropical water ice (to initialize high res sim) c -------------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'noglacier') then do ig=1,ngridmx j=(ig-2)/iim +2 if(ig.eq.1) j=1 write(*,*) 'OK: remove surface ice for |lat|<45' if (abs(rlatu(j)*180./pi).lt.45.) then qsurf(ig,nqmx)=0. end if end do c watercapn : H20 ice sur la calotte permanente nord c -------------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'watercapn') then do ig=1,ngridmx j=(ig-2)/iim +2 if(ig.eq.1) j=1 if (rlatu(j)*180./pi.gt.80.) then qsurf(ig,nqmx)=1.e5 write(*,*) 'ig=',ig,' H2O ice mass (kg/m2)= ', & qsurf(ig,nqmx) write(*,*)' ==> Ice mesh South boundary (deg)= ', & rlatv(j)*180./pi end if enddo c watercaps : H20 ice sur la calotte permanente sud c ------------------------------------------------- else if (modif(1:lnblnk(modif)) .eq. 'watercaps') then do ig=1,ngridmx j=(ig-2)/iim +2 if(ig.eq.1) j=1 if (rlatu(j)*180./pi.lt.-80.) then qsurf(ig,nqmx)=1.e5 write(*,*) 'ig=',ig,' H2O ice mass (kg/m2)= ', & qsurf(ig,nqmx) write(*,*)' ==> Ice mesh North boundary (deg)= ', & rlatv(j-1)*180./pi end if enddo c isotherm : Temperatures isothermes et vents nuls c ------------------------------------------------ else if (modif(1:lnblnk(modif)) .eq. 'isotherm') then write(*,*)'Isothermal temperature of the atmosphere, & surface and subsurface' write(*,*) 'Value of THIS temperature ? :' 203 read(*,*,iostat=ierr) Tiso if(ierr.ne.0) goto 203 do ig=1, ngridmx tsurf(ig) = Tiso end do do l=2,nsoilmx do ig=1, ngridmx tsoil(ig,l) = Tiso end do end do flagiso=.true. call initial0(llm*ip1jmp1,ucov) call initial0(llm*ip1jm,vcov) call initial0(ngridmx*(llm+1),q2) c co2ice=0 : ellimination des calottes polaires de CO2' c ------------------------------------------------ else if (modif(1:lnblnk(modif)) .eq. 'co2ice=0') then do ig=1,ngridmx co2ice(ig)=0 emis(ig)=emis(ngridmx/2) end do else goto 888 end if end do 999 continue c======================================================================= c Correction de la pression pour nouvelle grille (menu 0) c======================================================================= if (flagps0.eqv..false.) then r = 1000.*8.31/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 periodicite de ps en 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 (flagiso) then DO l=1,llm DO j=1,jjp1 DO i=1,iim teta(i,j,l) = Tiso * 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 c if (choix_1.eq.0) then CALL massdair( p3d, masse ) c 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 c 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 C Ecriture etat initial physique C call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx, . dtphys,float(day_ini), . time,tsurf,tsoil,co2ice,emis,q2,qsurf, . airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe) 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