PROGRAM gcm USE IOIPSL IMPLICIT NONE c ...... Version du 10/01/98 .......... c avec coordonnees verticales hybrides c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 ) c======================================================================= c c Auteur: P. Le Van /L. Fairhead/F.Hourdin c ------- c c Objet: c ------ c c GCM LMD nouvelle grille c c======================================================================= c c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv c et possibilite d'appeler une fonction f(y) a derivee tangente c hyperbolique a la place de la fonction a derivee sinusoidale. c ... Possibilite de choisir le shema de Van-leer pour l'advection de c q , en faisant iadv = 3 dans traceur (29/04/97) . c c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99) 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" INTEGER longcles PARAMETER ( longcles = 20 ) REAL clesphy0( longcles ) SAVE clesphy0 INTEGER*4 iday ! jour julien REAL time ! Heure de la journee en fraction d'1 jour REAL zdtvr INTEGER nbetatmoy, nbetatdem,nbetat c variables dynamiques REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants REAL teta(ip1jmp1,llm) ! temperature potentielle REAL q(ip1jmp1,llm,nqmx) ! champs advectes REAL ps(ip1jmp1) ! pression au sol 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 variables dynamiques intermediaire pour le transport REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse c variables dynamiques au pas -1 REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm) REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1) REAL massem1(ip1jmp1,llm) c tendances dynamiques REAL dv(ip1jm,llm),du(ip1jmp1,llm) REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1) c tendances de la dissipation REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm) REAL dhdis(ip1jmp1,llm) 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 REAL tppn(iim),tpps(iim),tpn,tps c INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq INTEGER itau,itaufinp1,iav EXTERNAL caldyn, traceur EXTERNAL dissip,geopot,iniconst,inifilr EXTERNAL integrd,SCOPY EXTERNAL iniav,writeav,writehist EXTERNAL inigeom EXTERNAL exner_hyb,addit EXTERNAL defrun_new, test_period REAL SSUM REAL time_0 , finvmaold(ip1jmp1,llm) LOGICAL lafin INTEGER ij,iq,l,numvanle,iapp_tracvl INTEGER fluxid, fluxvid,fluxdid integer histid, histvid, histaveid real time_step, t_wrt, t_ops REAL rdayvrai,rdaym_ini,rday_ecri LOGICAL first REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) LOGICAL offline ! Controle du stockage ds "fluxmass" PARAMETER (offline=.true.) character*80 dynhist_file, dynhistave_file character*20 modname character*80 abort_message C Calendrier LOGICAL true_calendar PARAMETER (true_calendar = .false.) c----------------------------------------------------------------------- c Initialisations: c ---------------- modname = 'gcm' descript = 'Run GCM LMDZ' lafin = .FALSE. dynhist_file = 'dyn_hist.nc' dynhistave_file = 'dyn_hist_ave.nc' if (true_calendar) then call ioconf_calendar('gregorian') else call ioconf_calendar('360d') endif c----------------------------------------------------------------------- c c .... Choix des shemas d'advection pour l'eau et les traceurs ... c ................................................................... c c iadv = 1 shema transport type "humidite specifique LMD" c iadv = 2 shema amont c iadv = 3 shema Van-leer c iadv = 4 schema Van-leer + humidite specifique c Modif F.Codron c c dans le tableau q(ij,l,iq) , iq = 1 pour l'eau vapeur c , iq = 2 pour l'eau liquide c et eventuellement , iq = 3, nqmx pour les autres traceurs c c iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq. c DO iq = 1, nqmx iadv( iq ) = 3 ENDDO c IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique' * ,' pour la vapeur d''eau' IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema amont',' pour la' * ,' vapeur d''eau ' IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour' * ,'la vapeur d''eau' IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema Van-Leer + humidite' * ,' specifique pour la vapeur d''eau' c IF( iadv(1).LE.0.OR.iadv(1).GT.4 ) THEN PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser * , car iadv(1) = ', iadv(1) CALL ABORT ENDIF DO iq = 2, nqmx IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique' * ,' pour le traceur no ', iq IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema amont',' pour le' * ,' traceur no ', iq IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour' * ,'le traceur no ', iq IF( iadv(iq).EQ.4 ) THEN PRINT *,' Le shema Van-Leer + humidite specifique ', * ' est uniquement pour la vapeur d eau .' PRINT *,' Corriger iadv( ',iq, ') et repasser ! ' CALL ABORT ENDIF IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 ) THEN PRINT *,' Erreur dans le choix de iadv . Corriger et repasser * . ' CALL ABORT ENDIF ENDDO c first = .TRUE. numvanle = nqmx + 1 DO iq = 1, nqmx IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN numvanle = iq first = .FALSE. ENDIF ENDDO c DO iq = 1, nqmx IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle ) THEN PRINT *,' Il y a discontinuite dans le choix du shema de ', * 'Van-leer pour les traceurs . Corriger et repasser . ' CALL ABORT ENDIF IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 ) THEN PRINT *,' Le choix de iadv est errone pour le traceur ', * iq STOP ENDIF ENDDO c CALL dynetat0("start.nc",nqmx,vcov,ucov, . teta,q,masse,ps,phis, time_0) c OPEN( 99,file ='run.def',status='old',form='formatted') CALL defrun_new( 99, .TRUE. , clesphy0 ) c CLOSE(99) 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' IF(MOD(day_step,iphysiq).NE.0) * 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 nbetatdem = nday / iecri nbetatmoy = nday / periodav + 1 dtvr = zdtvr CALL iniconst CALL inigeom CALL inifilr c c ...... P.Le Van ( modif le 29/04/97 ) ......... c CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , * tetagdiv, tetagrot , tetatemp ) c CALL pression ( ip1jmp1, ap, bp, ps, p ) CALL exner_hyb( ip1jmp1, ps, p,alpha,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 PRINT 300, itau,itaufin,day_ini,day_end CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx) ecripar = .TRUE. time_step = zdtvr t_ops = iecri * daysec t_wrt = iecri * daysec CALL inithist(dynhist_file,day_ini,anne_ini,time_step, . t_ops, t_wrt, nqmx, histid, histvid) t_ops = iperiod * time_step t_wrt = periodav * daysec CALL initdynav(dynhistave_file,day_ini,anne_ini,time_step, . t_ops, t_wrt, nqmx, histaveid) dtav = iperiod*dtvr/daysec c Quelques initialisations pour les traceurs call initial0(ijp1llm*nqmx,dq) istdyn=day_step/4 ! stockage toutes les 6h=1jour/4 istphy=istdyn/iphysiq c----------------------------------------------------------------------- c Debut de l'integration temporelle: c ---------------------------------- 1 CONTINUE c call nudge(itau,ucov,vcov,teta,masse,ps) c c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN CALL test_period ( ucov,vcov,teta,q,p,phis ) PRINT *,' ---- Test_period apres continue OK ! -----', itau c ENDIF c CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 ) CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 ) CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 ) CALL SCOPY( ijp1llm,masse, 1, massem1, 1 ) CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 ) forward = .TRUE. leapf = .FALSE. dt = dtvr c ... P.Le Van .26/04/94 .... CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 ) CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 ) 2 CONTINUE c----------------------------------------------------------------------- c date: c ----- c gestion des appels de la physique et des dissipations: c ------------------------------------------------------ c c ... P.Le Van ( 6/02/95 ) .... apphys = .FALSE. statcl = .FALSE. conser = .FALSE. apdiss = .FALSE. IF( purmats ) THEN IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE. IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE. IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward $ .AND. physic ) apphys = .TRUE. ELSE IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE. IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE. IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE. END IF c----------------------------------------------------------------------- c calcul des tendances dynamiques: c -------------------------------- CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) c c CALL caldyn $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis , $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini ) c----------------------------------------------------------------------- c calcul des tendances advection des traceurs (dont l'humidite) c ------------------------------------------------------------- IF( forward. OR . leapf ) THEN c DO 15 iq = 1, nqmx c IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq ) ELSE IF( iq.EQ. nqmx ) THEN c iapp_tracvl = 5 c cccc iapp_tracvl est la frequence en pas du groupement des flux cccc de masse pour Van-Leer dans la routine tracvl . c CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv, * p, masse, dq, iadv(1), teta, pk ) c c ... Modif F.Codron .... c ENDIF c 15 CONTINUE C IF (offline) THEN Cmaf stokage du flux de masse pour traceurs OFF-LINE CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, . time_step, itau,fluxid, fluxvid,fluxdid ) ENDIF c ENDIF c----------------------------------------------------------------------- c integrations dynamique et traceurs: c ---------------------------------- CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 , $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis , $ finvmaold ) c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd) c c----------------------------------------------------------------------- c calcul des tendances physiques: c ------------------------------- c ######## P.Le Van ( Modif le 6/02/95 ) ########### c IF( purmats ) THEN IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE. ELSE IF( itau+1. EQ. itaufin ) lafin = .TRUE. ENDIF c c IF( apphys ) THEN c c ....... Ajout P.Le Van ( 17/04/96 ) ........... c CALL pression ( ip1jmp1, ap, bp, ps, p ) CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf ) rdaym_ini = itau * dtvr / daysec rdayvrai = rdaym_ini + day_ini IF ( ecritphy.LT.1. ) THEN rday_ecri = rdaym_ini ELSE rday_ecri = NINT( rdayvrai ) ENDIF c CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time , $ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi , $ du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi ) c ajout des tendances physiques: c ------------------------------ CALL addfi( nqmx, dtphys, leapf, forward , $ ucov, vcov, teta , q ,ps , $ dufi, dvfi, dhfi , dqfi ,dpfi ) c ENDIF CALL pression ( ip1jmp1, ap, bp, ps, p ) CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) c----------------------------------------------------------------------- c c c dissipation horizontale et verticale des petites echelles: c ---------------------------------------------------------- IF(apdiss) THEN CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis) CALL addit( ijp1llm,ucov ,dudis,ucov ) CALL addit( ijmllm ,vcov ,dvdis,vcov ) CALL addit( ijp1llm,teta ,dhdis,teta ) c ....... P. Le Van ( ajout le 17/04/96 ) ........... c ... Calcul de la valeur moyenne, unique de h aux poles ..... c DO l = 1, llm DO ij = 1,iim tppn(ij) = aire( ij ) * teta( ij ,l) tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l) ENDDO tpn = SSUM(iim,tppn,1)/apoln tps = SSUM(iim,tpps,1)/apols DO ij = 1, iip1 teta( ij ,l) = tpn teta(ij+ip1jm,l) = tps ENDDO ENDDO DO ij = 1,iim tppn(ij) = aire( ij ) * ps ( ij ) tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm) ENDDO tpn = SSUM(iim,tppn,1)/apoln tps = SSUM(iim,tpps,1)/apols DO ij = 1, iip1 ps( ij ) = tpn ps(ij+ip1jm) = tps ENDDO END IF c ******************************************************************** c ******************************************************************** c .... fin de l'integration dynamique et physique pour le pas itau .. c ******************************************************************** c ******************************************************************** c preparation du pas d'integration suivant ...... IF ( .NOT.purmats ) THEN c ........................................................ c .............. schema matsuno + leapfrog .............. c ........................................................ IF(forward. OR. leapf) THEN itau= itau + 1 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 ENDIF IF( itau. EQ. itaufinp1 ) then abort_message = 'Simulation finished' call abort_gcm(modname,abort_message,0) ENDIF c----------------------------------------------------------------------- c ecriture du fichier histoire moyenne: c ------------------------------------- IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN IF(itau.EQ.itaufin) THEN iav=1 ELSE iav=0 ENDIF CALL writedynav(histaveid, nqmx, itau,vcov , , ucov,teta,pk,phi,q,masse,ps,phis) ENDIF c----------------------------------------------------------------------- c ecriture de la bande histoire: c ------------------------------ IF( MOD(itau,iecri*day_step).EQ.0) THEN nbetat = nbetatdem CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) CALL writehist( histid, histvid, nqmx, itau,vcov , , ucov,teta,phi,q,masse,ps,phis) ENDIF IF(itau.EQ.itaufin) THEN PRINT *,' Appel test_period avant redem ', itau 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 c----------------------------------------------------------------------- c gestion de l'integration temporelle: c ------------------------------------ IF( MOD(itau,iperiod).EQ.0 ) THEN GO TO 1 ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN IF( forward ) THEN c fin du pas forward et debut du pas backward forward = .FALSE. leapf = .FALSE. GO TO 2 ELSE c fin du pas backward et debut du premier pas leapfrog leapf = .TRUE. dt = 2.*dtvr GO TO 2 END IF ELSE c ...... pas leapfrog ..... leapf = .TRUE. dt = 2.*dtvr GO TO 2 END IF ELSE c ........................................................ c .............. schema matsuno ............... c ........................................................ IF( forward ) THEN itau = itau + 1 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 forward = .FALSE. IF( itau. EQ. itaufinp1 ) then abort_message = 'Simulation finished' call abort_gcm(modname,abort_message,0) ENDIF GO TO 2 ELSE IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN IF(itau.EQ.itaufin) THEN iav=1 ELSE iav=0 ENDIF CALL writedynav(histaveid, nqmx, itau,vcov , , ucov,teta,pk,phi,q,masse,ps,phis) ENDIF IF(MOD(itau,iecri*day_step).EQ.0) THEN nbetat = nbetatdem CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) CALL writehist( histid, histvid, nqmx, itau,vcov , , ucov,teta,phi,q,masse,ps,phis) ENDIF IF(itau.EQ.itaufin) . CALL dynredem1("restart.nc",0.0, . vcov,ucov,teta,q,nqmx,masse,ps) forward = .TRUE. GO TO 1 ENDIF END IF 300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, . 'c''est a dire du jour',i7,3x,'au jour',i7//) STOP END