PROGRAM nogcm IMPLICIT NONE c ...... Version du 01/09/15 .......... c avec coordonnees verticales hybrides c======================================================================= c c Auteur: F. Forget, T. Bertrand c ------- c c Objet: c ------ c c GCM LMD sans dynamique, pour modele saisonnier de surface c 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" 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 kp(ip1jmp1) ! TB15 = 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 p0 ! pression de reference REAL p00 ! globalaverage(kp) REAL qmean_ch4 ! mass mean mixing ratio vap ch4 REAL qmean_co ! mass mean mixing ratio vap co REAL pqch4(ip1jmp1) ! average methane mass index : ps*qch4 REAL pqco(ip1jmp1) ! average methane mass index : ps*q_co REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4 ! flag to set/remove calls to groupeun real globaverage2d !tau_n2, tau_ch4,tau_co 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. ) write(*,*) 'TB15: test0 r, cpp=',r,cpp 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' c IF(MOD(day_step,iphysiq).NE.0) c * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq' 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 c CALL inifilr c ...... P.Le Van ( modif le 29/04/97 ) ......... c CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh , c * tetagdiv, tetagrot , tetatemp ) 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 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(*,*) "TB16:",day_ini,nday,day_end 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 istdyn=day_step/4 ! stockage toutes les 6h=1jour/4 c istphy=istdyn/iphysiq c CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) c----------------------------------------------------------------------- c Debut de l integration temporelle: c ---------------------------------- c kp=exp(-z/H) Needed to define a reference pressure : c P0=pmean/globalaverage(kp) c thus P(i) = P0*exp(-z(i)/H) = P0*kp(i) c It is checked that globalaverage2d(Pi)=pmean DO ij=1,ip1jmp1 kp(ij) = exp(-phis(ij)/(r*38.)) ENDDO p00=globaverage2d(kp) ! mean pres at ref level c tau_n2 = 1. ! constante de rappel for pressure n2 (s) c tau_ch4 = 1.E7 ! constante de rappel for mix ratio qch4 (s) c tau_co = 1. !E5 ! constante de rappel for mix ratio qco (s) dt = dtvr dtphys = iphysiq * dtvr 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 c write(20,*) 'TB15 : ps1 nogcm= ',globaverage2d(ps) ! mean pressure 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) c write(21,*) 'TB15 : ps2 nogcm= ',globaverage2d(ps) ! mean pressure c saving pressure pplev : DO ij=1,ip1jmp1 oldps(ij)=ps(ij) ENDDO c increment q_ch4 with physical tendancy if (methane) then DO ij=1,ip1jmp1 q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)+ & dqfi(ij,1,igcm_ch4_gas)*dtphys ENDDO c write(31,*) 'TB15 : qch4 t1 = ', c & globaverage2d(q(:,1,igcm_ch4_gas)*ps(:)/g) endif c increment q_co with physical tendancy if (carbox) then DO ij=1,ip1jmp1 q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)+ & dqfi(ij,1,igcm_co_gas)*dtphys ENDDO ! write(41,*) 'TB15 : qco t1 = ', ! & globaverage2d(q(:,1,igcm_co_gas)*ps(:)/g) ENDIF c update pressure DO ij=1,ip1jmp1 ps(ij) = ps(ij) + dpfi(ij)*dtphys ENDDO psmean= globaverage2d(ps) ! mean pressure p0=psmean/p00 DO ij=1,ip1jmp1 ! Rappel newtonien vers psmean ps(ij)=ps(ij) +(p0*kp(ij) & -ps(ij))*(1-exp(-dtphys/tau_n2)) ENDDO write(72,*) 'TB15 : ps redistributed = ',psmean ! mean pressure write(72,*) ' ' c TB15 : test pressure negative DO ij=1,ip1jmp1 IF (ps(ij).le.0.) then write(*,*) 'Negative pressure :' write(*,*) 'ps = ',ps(ij),' ig = ',ij write(*,*) 'pmean = ',p0*kp(ij) write(*,*) 'tau_n2 = ',tau_n2 STOP ENDIF ENDDO if (methane) then ! Simulate redistribution by dynamics for qch4 DO ij=1,ip1jmp1 c update mmr q for new pressure c (both q and dqfi were calculated with old pressure in physiq) q(ij,1,igcm_ch4_gas)=q(ij,1,igcm_ch4_gas)*oldps(ij)/ & ps(ij) ! average methane mass index : ps * qch4 pqch4(ij)= ps(ij) * q(ij,1,igcm_ch4_gas) end do c write(32,*) 'TB15 : qch4 t2 = ', c & globaverage2d(q(:,1,igcm_ch4_gas)*ps(:)/g) ! weighted averaged CH4 mass mixing ratio: qmean_ch4= globaverage2d(pqch4) / globaverage2d(ps) DO ij=1,ip1jmp1 ! Rappel newtonien vers qmean 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)) ENDDO c write(33,*) 'TB15 : qch4 t3 = ', c & globaverage2d(q(:,1,igcm_ch4_gas)*ps(:)/g) endif if (carbox) then ! Simulate redistribution by dynamics for q_co DO ij=1,ip1jmp1 c update mmr q for new pressure c (both q and dqfi were calculated with old pressure in physiq) q(ij,1,igcm_co_gas)=q(ij,1,igcm_co_gas)*oldps(ij)/ & ps(ij) ! average methane mass index : ps * q_co pqco(ij)= ps(ij) * q(ij,1,igcm_co_gas) end do ! write(42,*) 'TB15 : qco t2 = ', ! & globaverage2d(q(:,1,igcm_co_gas)*ps(:)/g) ! weighted averaged CO mass mixing ratio: qmean_co= globaverage2d(pqco) / globaverage2d(ps) DO ij=1,ip1jmp1 ! Rappel newtonien vers qmean 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)) ENDDO ! write(43,*) 'TB15 : qco t3 = ', ! & globaverage2d(q(:,1,igcm_co_gas)*ps(:)/g) endif 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 day_end=int(day_ini+paleoyears*365.25/6.3872)+1 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 ENDIF 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 TBch4 : 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" REAL var(iip1,jjp1), globaverage2d integer i,j REAL airetot save airetot logical firstcall data firstcall/.true./ save firstcall if (firstcall) then firstcall=.false. c write(30,*) 'TB15 : aire = ',aire airetot =0. DO j=2,jjp1-1 ! lat DO i = 1, iim ! lon airetot = airetot + aire(i,j) END DO END DO airetot=airetot + aire(1,1)+aire(1,jjp1) c write(*,*) 'TB15 : nogcm totarea= ',airetot 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)*aire(1,1) globaverage2d = globaverage2d + var(1,jjp1)*aire(1,jjp1) globaverage2d = globaverage2d/airetot return end