PROGRAM nogcm IMPLICIT NONE c ...... Version du 01/09/15 .......... c avec coordonnees verticales hybrides c======================================================================= c Auteur: F. Forget, T. Bertrand c ------- c Objet: c ------ c GCM LMD sans dynamique, pour modele saisonnier de surface c======================================================================= c c----------------------------------------------------------------------- c Declarations: c ------------- #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 "control.h" #include "ener.h" #include "netcdf.inc" #include "description.h" #include "serre.h" #include "tracstoke.h" #include "sponge.h" #include "advtrac.h" #include "callkeys.h" #include "tracer.h" #include "dimphys.h" #include "surfdat.h" !#include "comgeom2.h" INTEGER*4 iday ! jour julien REAL time ! Heure de la journee en fraction d''1 jour REAL zdtvr c variables dynamiques REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle REAL q(ip1jmp1,llm,nqmx) ! champs advectes REAL ps(ip1jmp1) ! pression au sol REAL kpd(ip1jmp1) ! = exp (-z/H) REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches REAL pks(ip1jmp1) ! exner au sol REAL pk(ip1jmp1,llm) ! exner au milieu des couches REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches REAL masse(ip1jmp1,llm) ! masse d''air REAL phis(ip1jmp1) ! geopotentiel au sol REAL phi(ip1jmp1,llm) ! geopotentiel REAL w(ip1jmp1,llm) ! vitesse verticale c tendances dynamiques REAL dv(ip1jm,llm),du(ip1jmp1,llm) REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx) c tendances physiques REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm) REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1) c variables pour le fichier histoire REAL dtav ! intervalle de temps elementaire INTEGER itau,itaufinp1 EXTERNAL traceur EXTERNAL iniconst EXTERNAL SCOPY EXTERNAL inigeom EXTERNAL exner_hyb,addit EXTERNAL defrun_new, test_period REAL time_0 LOGICAL lafin INTEGER ij REAL rdayvrai,rdaym_ini,rday_ecri REAL beta(ip1jmp1,llm) character*20 modname character*80 abort_message INTEGER nq, numvanle REAL psmean ! pression moyenne REAL pqch4mean ! moyenne globale ps*qch4 REAL pqcomean ! moyenne globale ps*qco !REAL pqhazemean REAL pqprechazemean REAL p0 ! pression de reference REAL p00d ! globalaverage(kpd) REAL qmean_ch4 ! mass mean mixing ratio vap ch4 REAL qmean_co ! mass mean mixing ratio vap co !REAL qmean_haze REAL qmean_prechaze REAL pqch4(ip1jmp1) ! average methane mass index : ps*qch4 REAL pqco(ip1jmp1) ! average methane mass index : ps*q_co REAL pqhaze(ip1jmp1) REAL pqprechaze(ip1jmp1) REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4 ! flag to set/remove calls to groupeun real globaverage2d logical firstphys data firstphys/.true./ save firstphys c----------------------------------------------------------------------- c Initialisations: c ---------------- modname = 'gcm' descript = 'Run GCM LMDZ' lafin = .FALSE. c----------------------------------------------------------------------- c Initialize tracers using iniadvtrac (Ehouarn, oct 2008) CALL iniadvtrac(nq,numvanle) CALL dynetat0("start.nc",nqmx,vcov,ucov, . teta,q,masse,ps,phis, time_0) CALL defrun_new( 99, .TRUE. ) c on recalcule eventuellement le pas de temps IF(MOD(day_step,iperiod).NE.0) * STOP'Il faut choisir un nb de pas par jour multiple de iperiod' zdtvr = daysec/FLOAT(day_step) IF(dtvr.NE.zdtvr) THEN PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr ENDIF c nombre d etats dans les fichiers demarrage et histoire dtvr = zdtvr CALL iniconst CALL inigeom CALL inifilr ! Needed for exner_hyb for filtrez c ...... P.Le Van ( modif le 29/04/97 ) ......... CALL pression ( ip1jmp1, ap, bp, ps, p ) call dump2d(iip1,jjp1,ps,'PRESSION SURFACE') CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf ) c numero de stockage pour les fichiers de redemarrage: c----------------------------------------------------------------------- c temps de depart et de fin: c -------------------------- itau = 0 iday = day_ini+itau/day_step time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0 IF(time.GT.1.) THEN time = time-1. iday = iday+1 ENDIF itaufin = nday*day_step itaufinp1 = itaufin +1 day_end = day_ini + nday write(*,*) "Nogcm: day ini, end, nday:",day_ini,day_end,nday PRINT 300, itau,itaufin,day_ini,day_end 300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, . 'c''est a dire du jour',i7,3x,'au jour',i7//) CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx) ecripar = .TRUE. dtav = iperiod*dtvr/daysec c Quelques initialisations pour les traceurs call initial0(ijp1llm*nqmx,dq) c TB22 : call once CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) c----------------------------------------------------------------------- c Debut de l integration temporelle: c ---------------------------------- c kpd=exp(-z/H) Needed to define a reference pressure : c P0=pmean/globalaverage(kpd) c thus P(i) = P0*exp(-z(i)/H) = P0*kpd(i) c It is checked that globalaverage2d(Pi)=pmean DO ij=1,ip1jmp1 kpd(ij) = exp(-phis(ij)/(r*38.)) ENDDO p00d=globaverage2d(kpd) ! mean pres at ref level dt = dtvr dtphys = iphysiq * dtvr c --------- c LOOP c --------- DO itau = 1,itaufin apphys = ( MOD(itau+1,iphysiq ).EQ.0.AND. physic) c --------------------------------- c calcul des tendances physiques: c --------------------------------- IF( itau+1.EQ.itaufin) lafin = .TRUE. c IF( apphys ) THEN rdaym_ini = (itau) * dtvr / daysec rdayvrai = rdaym_ini + day_ini IF ( ecritphy.LT.1. ) THEN rday_ecri = rdaym_ini ELSE rday_ecri = INT(rdaym_ini)+INT(day_ini) ENDIF CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time , $ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi , $ du,dv,dteta,dq,w,dufi,dvfi,dhfi,dqfi,dpfi,tracer) ! 1) saving pressure pplev : DO ij=1,ip1jmp1 oldps(ij)=ps(ij) ENDDO ! 2) increment q with physical tendancy (mixing ratio kg/kg) DO ij=1,ip1jmp1 if (methane) then q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)+ & dqfi(ij,1,igcm_ch4_gas)*dtphys endif if (carbox) then q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)+ & dqfi(ij,1,igcm_co_gas)*dtphys endif if (fasthaze) then q(ij,1,igcm_haze)=q(ij,1,igcm_haze)+ & dqfi(ij,1,igcm_haze)*dtphys q(ij,1,igcm_prec_haze)=q(ij,1,igcm_prec_haze)+ & dqfi(ij,1,igcm_prec_haze)*dtphys endif ENDDO ! 3) update pressure DO ij=1,ip1jmp1 ps(ij) = ps(ij) + dpfi(ij)*dtphys ENDDO ! 4) Rappel newtonien vers psmean psmean= globaverage2d(ps) ! mean pressure p0=psmean/p00d DO ij=1,ip1jmp1 ps(ij)=ps(ij) +(p0*kpd(ij) & -ps(ij))*(1-exp(-dtphys/tau_n2)) ENDDO write(72,*) 'psmean = ',psmean ! mean pressure redistributed ! 4bis) security check: test pressure negative DO ij=1,ip1jmp1 IF (ps(ij).le.0.) then write(*,*) 'Negative pressure :' write(*,*) 'ps = ',ps(ij),' ig = ',ij STOP ps(ij)=1.e-6*kpd(ij)/p00d ENDIF ENDDO ! 5) Simulate redistribution by dynamics for q ! intermediaire de calcul si methane or CO if ((methane).or.(carbox)) then psmean= globaverage2d(ps) endif !5a) update mmr q for new pressure (both q and dqfi were calculated with old pressure in physiq) DO ij=1,ip1jmp1 if (methane) then q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)*oldps(ij)/ & ps(ij) ! average methane mass index : ps * qch4 => 'pressure of ch4' pqch4(ij)= ps(ij) * q(ij,1,igcm_ch4_gas) endif if (carbox) then q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)*oldps(ij)/ & ps(ij) pqco(ij)= ps(ij) * q(ij,1,igcm_co_gas) endif if (fasthaze) then q(ij,1,igcm_haze)=q(ij,1,igcm_haze)*oldps(ij)/ & ps(ij) pqhaze(ij)= ps(ij) * q(ij,1,igcm_haze) q(ij,1,igcm_prec_haze)=q(ij,1,igcm_prec_haze)*oldps(ij)/ & ps(ij) pqprechaze(ij)= ps(ij) * q(ij,1,igcm_prec_haze) endif ENDDO !5b) Rappel newtonien vers q_mean if (methane) then pqch4mean=globaverage2d(pqch4) qmean_ch4= pqch4mean / psmean endif if (carbox) then pqcomean=globaverage2d(pqco) qmean_co= pqcomean / psmean endif if (fasthaze) then pqprechazemean=globaverage2d(pqprechaze) qmean_prechaze= pqprechazemean / psmean !pqhazemean=globaverage2d(pqhaze) !qmean_haze= pqhazemean / psmean endif DO ij=1,ip1jmp1 if (methane) then q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)+ & (qmean_ch4-q(ij,1,igcm_ch4_gas))* & (1.-exp(-dtphys/tau_ch4)) endif if (carbox) then q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)+ & (qmean_co-q(ij,1,igcm_co_gas))* & (1.-exp(-dtphys/tau_co)) endif ! CO if (fasthaze) then ! q(ij,1,igcm_haze)=q(ij,1,igcm_haze)+ & ! (qmean_haze-q(ij,1,igcm_haze))* & ! (1.-exp(-dtphys/tau_haze)) q(ij,1,igcm_prec_haze)=q(ij,1,igcm_prec_haze)+ & (qmean_prechaze-q(ij,1,igcm_prec_haze))* & (1.-exp(-dtphys/tau_prechaze)) endif ! fasthaze ENDDO ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc CALL pression ( ip1jmp1, ap, bp, ps, p ) CALL massdair ( p , masse ) CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf ) if(firstphys) then c Option Special "nogcm" c --------------- corrk = .false. calldifv = .false. calladj = .false. callconduct=.false. IF (paleo) then IF (triton.or.kbo) then day_end = day_ini + nday ELSE day_end=int(day_ini+paleoyears*365.25/6.3872)+1 ENDIF !TB22 : maybe not needed. Leads to neg values of phi-phis !CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) ! appel a dynredem0 apes passage dans physiq pour ! recuperer phisfi (=>phis) et paleoyears) CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx) ENDIF write(*,*) '**********************************************' write(*,*) 'WARNING: For nogcm, these options are forced:' write(*,*) '**********************************************' write(*,*) 'Option corrk = ', corrk write(*,*) 'Option calldifv = ', calldifv write(*,*) 'Option calladj = ', calladj write(*,*) 'Option callconduct = ', callconduct write(*,*) 'tau N2 - CH4 - CO = ', tau_n2,tau_ch4,tau_co write(*,*) '**********************************************' firstphys=.false. end if ! firstphys ENDIF !apphys c ******************************************************************** c .... fin de l'integration dynamique et physique pour le pas itau .. c ******************************************************************** c preparation du pas d'integration suivant ...... iday= day_ini+itau/day_step time= FLOAT(itau+1-(iday-day_ini)*day_step)/day_step+time_0 IF(time.GT.1.) THEN time = time-1. iday = iday+1 ENDIF c----------------------------------------------------------------------- IF(itau.EQ.itaufin) THEN PRINT *,' Appel test_period avant redem ', itau c TB : coupe test period car nous embete pour des pouilleme c sur le ch4_gas CALL test_period ( ucov,vcov,teta,q,p,phis ) CALL dynredem1("restart.nc",0.0, . vcov,ucov,teta,q,nqmx,masse,ps) CLOSE(99) ENDIF END DO ! fin de la boucle temporelle abort_message = 'Simulation finished' call abort_gcm(modname,abort_message,0) STOP END !********************************************************************************* FUNCTION globaverage2d(var) IMPLICIT NONE #include "dimensions.h" #include "paramet.h" #include "comgeom2.h" !ip1jmp1 called in paramet.h = iip1 x jjp1 REAL var(iip1,jjp1), globaverage2d integer i,j REAL airetot save airetot logical firstcall data firstcall/.true./ save firstcall if (firstcall) then firstcall=.false. airetot =0. DO j=1,jjp1 ! lat DO i = 1, iim ! lon airetot = airetot + aire(i,j) END DO END DO end if globaverage2d = 0. DO j=2,jjp1-1 DO i = 1, iim globaverage2d = globaverage2d + var(i,j)*aire(i,j) END DO END DO globaverage2d = globaverage2d + var(1,1)*sum(aire(1:iim,1)) globaverage2d = globaverage2d + var(1,jjp1)*sum(aire(1:iim,jjp1)) globaverage2d = globaverage2d/airetot return end !*********************************************************************************